Trigonometric argument reduction (reduction modulo 2π)

Question

I am trying to create a calculator for sin and cosine that technically only operates on the range of 0-pi/2. Right now this might seem silly but later on it will be used so I can employ Taylor series.

I have a mostly working implementation however I have a serious issue when theta is in the form of x * (pi/2) where x is an arbitrary integer. It would appear that on these values that sometimes they are pushed into a nearby quadrants they dont belong to. There are also some occasional outright errors that I cannot explain.

How can I shore this up making it more efficient and correct?

Here is the code for doing this.

#define T_PI (2.0 * M_PI)
#define H_PI (0.5 * M_PI)
void sincos(float theta, float* cosine, float* cosine) {
  int mode;
  prepareForRange(&theta, cosine, sine);
  Assert(!(f < 0.0 || f > H_PI));
  *cosine = cos(theta);
  *sine = sin(theta);
  range_output(mode, cosine, sine);
}
void prepareForRange(float* theta, int* mode, float *cosine, float* sine) {
  if (*theta < 0.0) *theta += ceil(-*theta / T_PI) * T_PI;
  *mode = (int)floor(*theta / H_PI) % 4 + 1;
  *theta = fmodf(*theta, H_PI);
}
void range_output(int mode, float *cos, float *sin) {
  float temp;
  switch (mode) {
    case 1:
      break;
    case 2:
      temp = *cos;
      *cos = -*sin;
      *sin = temp;
      break;
    case 3:
      *cos = -*cos;
      *sin = -*sin;
      break;
    case 4:
      temp = *cos;
      *cos = *sin;
      *sin = -temp;
      break;
    default:
      break;
    }
}

When you post a question like this, you should include at least one input case for which it does not produce the desired result, along with the results you get from the code and the results you expect. The code is likely suffering from floating-point rounding issues. Floating-point does not represent real numbers exact. In particular, `M_PI` cannot be exactly π, and working with multiples of `M_PI` will multiply the deviation. See [the answers to this question](https://stackoverflow.com/questions/11576202/modulo-2pi-using-sse-sse2) for a lead on a proper algorithm. — Eric Postpischil, Apr 29 '18 at 00:30
My method is to use a similar reduction but involving **rational numbers**: `k = prepareForRange(int x2)`. Then your function sincos(x) would be of the form `sincos(k*pi)` where `k` is the trigonometrically reduced fraction. — Ṃųỻịgǻňạcểơửṩ, Apr 29 '18 at 01:46

score 1 · Answer 1 · answered Apr 29 '18 at 19:58

1

You are running into a long-standing and often-unrecognized problem area called range reduction. The basic problem is that a float PI constant is only defined 8 digit accuracy, so by the time you are trying to compute (x - n * PI) for n ~ 10^4, you have lost half the digits of accuracy in your results, and worse as n gets larger. There is no simple software solution to this problem. To really solve it in my own numerical library, I had to effectively implement arbitrary-precision floating point arithmetic and store a 320-digit (1078-bit) PI constant. Some implementations of libc effectively do this for you, but not all, so you can't safely assume it.

answered Apr 29 '18 at 19:58

David Wright

765
6
15

`a float PI constant is only defined 8 digit accuracy` where did you get that? The OP is declaring a `double` M_PI, which is ~15-17 digits of precision, and it's only cast down to `float` in variables and arguments. Anyway a single-precision `float` has only ~7 digits if precision, not 8. And note that [accuracy and precision are different](https://stackoverflow.com/q/8270789/995714) – phuclv Apr 30 '18 at 00:45
@LưuVĩnhPhúc: I grant you all your nit-picks. Of course, none of them are central to the issue, which is that his PI and the arithmetic operations using it have nowhere near the required precision to do range reduction. – David Wright Apr 30 '18 at 00:58

Trigonometric argument reduction (reduction modulo 2π)

1 Answers1

Linked