3

I have been reading C99 Rationale, where I found this mysterious excerpt (emphasis is mine):

7.12.4 Trigonometric functions

Implementation note: trignometric argument reduction should be performed by a method that causes no catastrophic discontinuities in the error of the computed result. In particular, methods based solely on naive application of a calculation like

x - (2*pi) * (int)(x/(2*pi))

are ill-advised.

What exactly is wrong with this formula of reduction? It seems, that it looks fine, according to property of periodity, with the interval of 2*pi.

Community
  • 1
  • 1
Grzegorz Szpetkowski
  • 36,988
  • 6
  • 90
  • 137
  • 1
    For one thing, it gets progressively worse as *x* increases, because `2 * pi` is not equal to 2π for any representable floating point value `pi`. – John Bollinger Mar 16 '17 at 20:36

1 Answers1

2

π is an irrational number unable to be represented exactly by a finite floating-point value - which are all rational.

Various implementation support a constant like M_PI which is nearly, but not exactly π. So the following introduces error. Of course it is a problem if (x/(2*pi) exceed the int range.

double pi = M_PI;
double x;  // radians
double y;  // reduced radians.
y  = x - (2*pi) * (int)(x/(2*pi))

If this error is important to code is application specific. The typical issue is with tan(x) where x is near π*(n +1/2) and a slight variation on x will generate + or - infinity/DBL_MAX.

Some platforms supply functions for π reduction.

A good reference for this problem is ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit K. C. Ng and the members of the FP group of SunPro

To reduce in degrees:

A range reduction for degrees is fmod(x,360.0) which can be expected to reduce x to the range -360.0 < x < +360.0 exactly. Better to use remquo: sind() example

Community
  • 1
  • 1
chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
  • 1
    @Grzegorz Szpetkowski Note that `y = fmod(x, 2*M_PI);` is certainly better than `y = x - (2*pi) * (int)(x/(2*pi))`. If that is good enough, is application specific. – chux - Reinstate Monica Mar 16 '17 at 20:51
  • Hmm, I noticed `x - (2*pi) * (int)(x/(2*pi))` also introduces large cancellation for values _near_ multiples of 2*π. The inexactness of `M_PI`, multiplication (with its round-off error) and that subtractive cancellation of nearly the same significant bits can leave the remaining bits meaningless. With a trig function with a discontinuity near multiples of 2*π ([cotan()](https://www.google.com/search?client=opera&q=cotangent+grpah&sourceid=opera&ie=UTF-8&oe=UTF-8#q=cotangent+graph&*)), this is problematic. – chux - Reinstate Monica Mar 16 '17 at 21:06
  • The decision to use sin(x) as its fundamental primitive rather than sin(2πx) and 1-cos(2πx) means that there's no good way for programs to precisely compute the above functions for precisely-representable values of x. I wonder how many non-contrived situations arise where code which is interested in accuracy should ever evaluate trig functions for radian-measured angles greater than π/2? Even if one wants to compute e.g. sin(x/1000000) [i.e. a function with a period of 2000000π] accuracy would demand that one perform argument reduction on x before doing the division (which would of course... – supercat Mar 17 '17 at 17:33
  • ...imply that it was done before doing the call]. I suspect that more often than not, code which computes sin(x) has just scaled x by some multiple of π, and that what code really wants to compute is the mathematical value of sin(x * π / roundedPi). I know that radians are useful in calculus, but in most other applications precisely-representable fractions of circles are more useful. – supercat Mar 17 '17 at 17:36
  • @supercat Sounds like a `sin2pi(revolutions)`, etc. would be a good companion to `sin(radians)` for the situations you describe. Much like `log(x)` has companions: `log10(x), log2()`. The [referenced article](https://www.csee.umbc.edu/~phatak/645/supl/Ng-ArgReduction.pdf) details the way to counter "no good way for programs to precisely compute the above functions (like `sin(x)`)". – chux - Reinstate Monica Mar 17 '17 at 17:56
  • @supercat "but in most other applications precisely-representable fractions of circles are more useful." cuts to the core of good function design. There are many goals: correctness, speed, code/memory footprint, broadness of application. `sin()` is a good compromise, but as you have well identified, comes up short with a number of practical applications. Thanks for your insights. – chux - Reinstate Monica Mar 17 '17 at 18:02
  • I don't see how that algorithm would be helpful to someone who wants to e.g. compute sin(2πx) within 3/4ulp for values of x near 1/8. Multiplication of values near 1/8 by 2π would cause an error of ±1/2ulp in the argument, which would affect the result by about 0.35lsb in addition to the ±1/2 lsb introduced by rounding the result, leaving a net uncertainty of about 0.85lsb. – supercat Mar 17 '17 at 18:22
  • Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/138351/discussion-between-chux-and-supercat). – chux - Reinstate Monica Mar 17 '17 at 18:27