Sin and Cos Give Unexpected Results For Well-Known Angles

Sin and Cos give unexpected results for well-known angles

C/C++ provides sin(a), cos(a), tan(a), etc. functions that require a parameter with radian units rather than degrees. double DegreesToRadians(d) performs a conversion that is close but an approximate as the conversion results are rounded. Also machine M_PI is close, but not the same value as the the mathematical irrational π.

OP's code with 180 passed to DegreesToRadians(d) and then to sin()/cos() gives results that differ than expected due to rounding, finite precision of double() and possible a weak value for PI.

An improvement is to perform argument reduction in degrees before calling the trig function. The below reduces the angle first to a -45° to 45° range and then calls sin(). This will insure that large values of N in sind(90.0*N) --> -1.0, 0.0, 1.0. . Note: sind(360.0*N +/- 30.0) may not exactly equal +/-0.5. Some additional considerations needed.

#include <math.h>
#include <stdio.h>

static double d2r(double d) {
  return (d / 180.0) * ((double) M_PI);
}

double sind(double x) {
  if (!isfinite(x)) {
    return sin(x);
  }
  if (x < 0.0) {
    return -sind(-x);
  }
  int quo;
  double x90 = remquo(fabs(x), 90.0, &quo);
  switch (quo % 4) {
    case 0:
      // Use * 1.0 to avoid -0.0
      return sin(d2r(x90)* 1.0);
    case 1:
      return cos(d2r(x90));
    case 2:
      return sin(d2r(-x90) * 1.0);
    case 3:
      return -cos(d2r(x90));
  }
  return 0.0;
}

int main(void) {
  int i;
  for (i = -360; i <= 360; i += 15) {
    printf("sin()  of %.1f degrees is  % .*e\n", 1.0 * i, DBL_DECIMAL_DIG - 1,
        sin(d2r(i)));
    printf("sind() of %.1f degrees is  % .*e\n", 1.0 * i, DBL_DECIMAL_DIG - 1,
        sind(i));
  }
  return 0;
}

Output

sin()  of -360.0 degrees is   2.4492935982947064e-16
sind() of -360.0 degrees is  -0.0000000000000000e+00  // Exact

sin()  of -345.0 degrees is   2.5881904510252068e-01  // 76-68 = 8 away
//                            2.5881904510252076e-01
sind() of -345.0 degrees is   2.5881904510252074e-01  // 76-74 = 2 away

sin()  of -330.0 degrees is   5.0000000000000044e-01  // 44 away
//  0.5                       5.0000000000000000e-01
sind() of -330.0 degrees is   4.9999999999999994e-01  //  6 away

sin()  of -315.0 degrees is   7.0710678118654768e-01  // 68-52 = 16 away
// square root 0.5 -->        7.0710678118654752e-01
sind() of -315.0 degrees is   7.0710678118654746e-01  // 52-46 = 6 away

sin()  of -300.0 degrees is   8.6602540378443860e-01
sind() of -300.0 degrees is   8.6602540378443871e-01
sin()  of -285.0 degrees is   9.6592582628906842e-01
sind() of -285.0 degrees is   9.6592582628906831e-01
sin()  of -270.0 degrees is   1.0000000000000000e+00  // Exact
sind() of -270.0 degrees is   1.0000000000000000e+00  // Exact
...

Sin cos function giving wrong results for degrees greater than 90

Functions cos and sin receive angles in radians instead of degrees.

Implementation of sine function in C not working

Each time your for loop progresses, n is increased by 2 and hence for DEPTH = 16, near the end of loop you have to calculate factorials of numbers as big as 30 and you are using unsigned int that can only store values as big as 2^32 = 4294967296 ~= 12! and this causes overflow in your factorial function which in turn gives you the wrong factorial.

Even if you used long double for it and I already stated in my comments that long double in MSCRT is mapped to double (Reference) you'd still see some anomalies probably at larger angles because although double can store values as big as 1.8E+308 but it loses its granularity at 2^53 = 9007199254740992 ~= 18! (i.e. 2^53 + 1 stored as a double is equal to 2^53). So once you go up in angles, the effect of this behavior becomes larger and larger to the point that it is noticeable in the 6 decimal precision that you are using with printf().

Although you are on the right track, you should use a bignum library like GMP or libcrypto. They can perform these calculations without the loss of precision.

BTW, since your are developing on Windows 7 that means you are either using x86 or x86-64. On these platforms, x87 is capable of performing extended precision (as per 754 standard) operations with 80 bits but I am unaware of compiler intrinsics that can give you that capability without resorting to assembly code.

I also would like to direct your attention to range reduction techniques. Although I still recommend using bignum libs, if you are good between 0 and 90 degrees (0 and 45 if I'm to be more strict), you can compute the sine() of all other angles just by simple trigonometric identities.

UPDATE:

Actually I'm gonna correct myself about using doubles in factorial calculations. After writing a simple program I verified that when I usedouble to store factorials, they are correct even if I go upper than 18. After giving it some thought I realized that in the case of factorials, the situation with double's granularity is a little bit more complex. I'll give you an example to make it clear:

19! = 19 * 18 * ... * 2 * 1

in this number 18, 16, 14, ... , 2 are all multiples of 2 and since a multiplication by 2 is equivalent to a shift to the left in binary representation, all lower bits in 19! are already 0 and hence when double's rounding kicks in for integers greater than 2^53, these factorials are unaffected. You can compute the number of least significant zeroes in the binary representation of 19! by counting the number of 2's which is 16. (for 20!, it is 18)

I'm gonna go up to 1.8e+308 and check if all the factorials are unaffected or not. I'll update you with the results.

UPDATE 2:

If we use doubles to hold factorials, they are affected by rounding from 23! onward. It can be easily shown, because 2^74 < 23! < 2^75 which means that at least 75 bits of precision is required to represent it, but since 23! has 19 least significant bits with the value of 0, so it needs 75 - 19 = 56 which is larger than 53 bits provided by double.

For 22!, it is 51 bits (you can calculate it yourself).

What discontinuities in the error are present in this trigonometric argument reduction?

π 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

Conversion of Degrees to Radians for use in trig functions in C

They probably will affect the accuracy by the precision to which you have pi - usually one unit in the last place. For cos and sin, |d f(x)/dx| is less than one so your value should be within the same error. For tan, |d f(x)/dx| is not bounded so a small input change can create a large change of output.

Whether such small changes will affect the operation of your code depends largely on whether your code assumes than the results are exact or not. If your code makes faulty assumptions about floating point values, then it will fail, if it allows some small tolerance on equality comparisons, then it wont.

Cosine Taylor approximation C language

OP's use of the Taylor's series suffers from numerical limitations with increasing values of x. The addition of large alternating sign terms accumulate too much error. Amend code to see these terms.

t = t*(-1)*x*x/(k*(k+1));
printf("%e\n", t);

Both sine and cosine calculation benefit with argument reduction to a
range of [-2π ... +2π]. The following is a good 1st step and will provide reduce the error for x over wider range.

x = fmod(x, 2*π);

Further reduction to the range [0 ... +π/4] can be employed using the usual trigonometric identities and remquo(). Degrees example

---

The trouble is that the above calculation relies on an approximate value of π. π, being an irrational number, cannot be represented exactly as a double. All finite double are rational numbers. So instead machine pi is used.

// example
#define M_PI 3.1415926535897932384626433832795
x = fmod(x, M_PI);

---

To achieve accurate calculations over the entire range of double x requires sophisticated extended precision techniques. Search for Argument reduction for huge arguments: Good to the last bit

---

Related Topics