Implementation of sine function in C not working

Question

I have attempted to implement the sine function in C, yet I am getting weird results. Here are the three functions I am using to calculate sine:

#define PI 3.14159265358979323846
#define DEPTH 16

double sine(long double);
long double pow(long double, unsigned int);
unsigned int fact(unsigned int);

double sine(long double x) {
    long double i_x = x *= PI/180;
    int n = 3, d = 0, sign = -1;

    // fails past 67 degrees
    for (; d < DEPTH; n += 2, d++, sign *= -1) {
        x += pow(i_x, n) / fact(n) * sign;
    }

    return x;
}

long double pow(long double base, unsigned int exp) {
    double answer = 1;
    while (exp) {
        answer *= base;
        exp--;
    }
    return answer;
}

unsigned int fact(unsigned int n) {
    unsigned int answer = 1;
    while (n > 1) {
        answer *= n--;
    }
    return answer;
}

To test it I have been testing it against the built in sine function as follows:

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

main() {
    for (int i = 0; i <= 180; i++) {
        printf("sin(%i) = %lf, %lf\n", i, sine(i), sin(i*3.14159265358979323846/180));
    }

    exit(EXIT_SUCCESS);
}

Up through 67 degrees, it calculates the same as the built in function. Though, as it increases past 67 it typically strays gets further and further from the actual value.

Here is an example output:

>> sin(100) = 0.987711, 0.984808
>> sin(101) = 0.986885, 0.981627
>> sin(102) = 0.987056, 0.978148
>> sin(103) = 0.988830, 0.974370
>> sin(104) = 0.993060, 0.970296
>> sin(105) = 1.000948, 0.965926
>> sin(106) = 1.014169, 0.961262
>> sin(107) = 1.035052, 0.956305
>> sin(108) = 1.066807, 0.951057
>> sin(109) = 1.113846, 0.945519
>> sin(110) = 1.182194, 0.939693
>> sin(111) = 1.280047, 0.933580
>> sin(112) = 1.418502, 0.927184
>> sin(113) = 1.612527, 0.920505
>> sin(114) = 1.882224, 0.913545
>> sin(115) = 2.254492, 0.906308
>> sin(116) = 2.765192, 0.898794
>> sin(117) = 3.461969, 0.891007
              ...
>> sin(180) = 8431648.192239, 0.000000

Does anybody know why this is happening? I am using Visual Studio 2017 on Windows 7, if that provides any useful information.

Answer 1

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).

Answer 2

There are multiple issues in your code:

• You redefine the standard function pow() with a different prototype. This may cause problems when you link the program as en executable. Use a different anme such as pow_int.

• You should define the pow_int and fact functions as static before the sine function. It may allow for better optimisation at compile time.

• Indeed fact is limited by the range of type unsigned int which is much less than the precision of type long double. Factorials beyond 12 have an incorrect value, causing a loss of precision.

• You could actually compute the terms incrementally, saving a lot of computations and avoiding potential loss of precision.

• The prototype for main() without arguments is int main(void)

• The computation of PI/180 is performed as double, which is less precise than long double. You should write the expression as x = x * PI / 180;

• DEPTH should be increased to improve the precision. At least 4 more terms bring a substantial improvement.

• You should apply a range reduction: taking advantage of the sine function symmetric and periodic nature, computation can be performed with fewer terms on x modulo 90 or even 45 degrees.

Here is a modified version:

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

#define PI_L   3.14159265358979323846264338327950288L
#define PI     3.14159265358979323846264338327950288
#define DEPTH  24

double sine(long double x) {
    long double res, term, x2, t1;
    int phase;

    x = remquol(x, 90, &phase);
    if (phase & 1)
        x = 90 - x;

    x = x * PI_L / 180; // convert x to radians
    x2 = x * x;         // pre-compute x^2

    // compute the sine series: x - x^3/3! + x^5/5! ...
    res = term = x;   // the first term is x
    for (int n = 1; n < DEPTH; n += 4) {
        // to reduce precision loss, compute 2 terms for each iteration
        t1 = term * x2 / ((n + 1) * (n + 2));
        term = t1 * x2 / ((n + 3) * (n + 4));
        // update the result with the difference of the terms
        res += term - t1;
    }
    if (phase & 2)
        res = -res;

    return (double)res;
}

int main(void) {
    printf("deg            sin                  sine         delta\n\n");
    for (int i = 0; i <= 360; i += 10) {
        double s1 = sin(i * PI / 180);
        double s2 = sine(i);
        printf("%3i  %20.17f  %20.17f  %g\n", i, s1, s2, s2 - s1);
    }
    return 0;
}

The output is:

deg            sin                  sine         delta

  0   0.00000000000000000   0.00000000000000000  0
 10   0.17364817766693033   0.17364817766693036  2.77556e-17
 20   0.34202014332566871   0.34202014332566871  0
 30   0.49999999999999994   0.50000000000000000  5.55112e-17
 40   0.64278760968653925   0.64278760968653936  1.11022e-16
 50   0.76604444311897801   0.76604444311897801  0
 60   0.86602540378443860   0.86602540378443860  0
 70   0.93969262078590832   0.93969262078590843  1.11022e-16
 80   0.98480775301220802   0.98480775301220802  0
 90   1.00000000000000000   1.00000000000000000  0
100   0.98480775301220802   0.98480775301220802  0
110   0.93969262078590843   0.93969262078590843  0
120   0.86602540378443882   0.86602540378443860  -2.22045e-16
130   0.76604444311897812   0.76604444311897801  -1.11022e-16
140   0.64278760968653947   0.64278760968653936  -1.11022e-16
150   0.49999999999999994   0.50000000000000000  5.55112e-17
160   0.34202014332566888   0.34202014332566871  -1.66533e-16
170   0.17364817766693025   0.17364817766693036  1.11022e-16
180   0.00000000000000012  -0.00000000000000000  -1.22465e-16
190  -0.17364817766693047  -0.17364817766693036  1.11022e-16
200  -0.34202014332566866  -0.34202014332566871  -5.55112e-17
210  -0.50000000000000011  -0.50000000000000000  1.11022e-16
220  -0.64278760968653925  -0.64278760968653936  -1.11022e-16
230  -0.76604444311897790  -0.76604444311897801  -1.11022e-16
240  -0.86602540378443837  -0.86602540378443860  -2.22045e-16
250  -0.93969262078590821  -0.93969262078590843  -2.22045e-16
260  -0.98480775301220802  -0.98480775301220802  0
270  -1.00000000000000000  -1.00000000000000000  0
280  -0.98480775301220813  -0.98480775301220802  1.11022e-16
290  -0.93969262078590854  -0.93969262078590843  1.11022e-16
300  -0.86602540378443860  -0.86602540378443860  0
310  -0.76604444311897812  -0.76604444311897801  1.11022e-16
320  -0.64278760968653958  -0.64278760968653936  2.22045e-16
330  -0.50000000000000044  -0.50000000000000000  4.44089e-16
340  -0.34202014332566855  -0.34202014332566871  -1.66533e-16
350  -0.17364817766693127  -0.17364817766693036  9.15934e-16
360  -0.00000000000000024   0.00000000000000000  2.44929e-16

As can be seen above, the sine() function seems more precise than the standard sin function on my system: sin(180 * M_PI / 128) should be 0 precisely. Similarly, sin(150 * M_PI / 128) should be 0.5 exactly.

Answer 3

Your way of polynomial series evaluation is numerically unstable. Try horner's method which is more stable than power calculations.