Sine function implementation in C does not work

Question

I tried to implement the sine function in C, but I get weird results. Here are three functions that I use to calculate the 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 this, I tested 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); } 

To the level of 67 degrees, it calculates the same as the built-in function. Although, as it rises to 67, it usually suffers from actual value further and further.

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 anyone know why this is happening? I am using Visual Studio 2017 on Windows 7 if this provides any useful information.

Answer 1

Every time your for loop progresses, n increases by 2 and therefore by DEPTH = 16 , near the end of the loop you should calculate factorials with numbers greater than 30 , and you use unsigned int , which can only store such values like 2^32 = 4294967296 ~= 12! , and this causes an 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 displayed on double ( Link ), you will still see some anomalies, probably from large angles, because although double can store values larger than 1.8E+308 , but it loses its grain at 2^53 = 9007199254740992 ~= 18! (i.e. 2^53 + 1 stored as double is equal to 2^53 ). Therefore, as soon as you rise at an angle, the effect of this behavior becomes more and more to such an extent that it is noticeable in 6 decimal precision, which you use with printf() .

While you're on the right track, you should use a bignum library like GMP or libcrypto . They can perform these calculations without loss of accuracy.

By the way, since you are running Windows 7, this means that you are using x86 or x86-64. On these platforms, x87 is capable of performing extended precision (in accordance with 754 standard) operations with 80 bits, but I do not know about the built-in compiler scripts that can provide you with this feature without resorting to assembler.

I would also like to draw your attention to range reduction methods. Although I still recommend using bignum libs if you are good between 0 and 90 degrees ( 0 and 45 if I will be more strict), you can calculate sine() all other angles simply using simple trigonometric identities .

UPDATE:

In fact, I'm going to correct myself about using double in factorial calculations. After writing a simple program, I confirmed that when I use double to store factorials, they are true even if I believe above 18 . Thinking, I realized that in the case of factorials, the situation with double granularity is a little more complicated. I will give you an example to make it clear:

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

including 18, 16, 14, ... , 2 all multiples of 2 , and since multiplying by 2 equivalent to a left shift in binary representation, all the least significant bits are 19! already 0 and, therefore, when double rounding works for integers greater than 2^53 , these factorials are not affected. You can calculate the number of least significant zeros in binary representation 19! counting the number 2 , which is 16 . (for 20! it's 18 )

I am going to go up to 1.8E+308 and check if all factorials are affected or not. I will tell you about the results.

UPDATE 2:

If we use double to store factorials, they are affected by rounding from 23! forward. This is easy to show because 2^74 < 23! < 2^75 2^74 < 23! < 2^75 means that it requires at least 75 bits of precision to represent it, but since 23! has 19 least significant bits with a value of 0 , so it needs 75 - 19 = 56 , which is more than 53 bits provided by double .

For 22! this is bit 51 (you can calculate it yourself).

Answer 2

There are several problems in the code:

• You will redefine the standard pow() function with another prototype. This can cause problems when you link the program as executable. Use another anme, for example pow_int .

• You must define the pow_int and fact functions as static before the sine function. This may provide better compilation time optimization.

• Indeed, fact limited to a range of type unsigned int , which is much less accurate than long double . Factorials beyond 12 have the wrong value, which leads to a loss of accuracy.

• You could compute terms gradually, saving a lot of computation and avoiding the potential loss of accuracy.

• The prototype for main() with no arguments is int main(void)

• The PI/180 calculation is performed as a double , which is less accurate than a long double . You must write the expression as x = x * PI / 180;

• DEPTH should be increased to increase accuracy. At least four more conditions significantly improve.

• You must apply a range reduction: using a sinusoidal function of a symmetric and periodic nature, the calculation can be performed with fewer x terms 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; } 

Output:

 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 you can see above, the sine() function seems more accurate than the standard sin function in my system: sin(180 * M_PI / 128) must be 0 sure. Similarly, sin(150 * M_PI / 128) should be 0.5 sure.

Answer 3

Your method for estimating polynomial series is numerically unstable . Try the horner method , which is more stable than power calculations.

Answer 4

Your problem is here:

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

You use d < DEPTH in error, when it should be n < DEPTH , d not related to your calculations in the loop. The following should work - although I am not compiled for testing.

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

Note: a DEPTH of 12 (for example, expanding the Taylor series with members 1, 3, 5, ... 11 ) is enough for a 3e-10 error - 3 ten billion in 60-degrees . (although the error increases with increasing angle between 0-360 , a DEPTH 20 will contain an error less than 1.0e-8 in the entire range.)

Enabling compiler warnings might catch on an unused d in sine .

Here is a sample code with changes (note: Gnu provides the M_PI constant for PI):

 #include <stdio.h> #include <stdint.h> #include <math.h> #define DEPTH 16 /* n factorial */ uint64_t nfact (int n) { if (n <= 0) return 1; uint64_t s = n; while (--n) s *= n; return s; } /* y ^ x */ double powerd (const double y, const int x) { if (!x) return 1; double r = y; for (int i = 1; i < x; i++) r *= y; return r; } double sine (double deg) { double rad = deg * M_PI / 180.0, x = rad; int sign = -1; for (int n = 3; n < DEPTH; n += 2, sign *= -1) x += sign * powerd (rad, n) / nfact (n); return x; } int main (void) { printf (" deg sin sine\n\n"); for (int i = 0; i < 180; i++) printf ("%3d %11.8f %11.8f\n", i, sin (i * M_PI / 180.0), sine (i)); return 0; } 

Usage / Output Example

 $ ./bin/sine deg sin sine 0 0.00000000 0.00000000 1 0.01745241 0.01745241 2 0.03489950 0.03489950 3 0.05233596 0.05233596 4 0.06975647 0.06975647 5 0.08715574 0.08715574 6 0.10452846 0.10452846 7 0.12186934 0.12186934 8 0.13917310 0.13917310 9 0.15643447 0.15643447 10 0.17364818 0.17364818 11 0.19080900 0.19080900 12 0.20791169 0.20791169 13 0.22495105 0.22495105 14 0.24192190 0.24192190 15 0.25881905 0.25881905 16 0.27563736 0.27563736 17 0.29237170 0.29237170 18 0.30901699 0.30901699 19 0.32556815 0.32556815 20 0.34202014 0.34202014 21 0.35836795 0.35836795 22 0.37460659 0.37460659 23 0.39073113 0.39073113 24 0.40673664 0.40673664 25 0.42261826 0.42261826 26 0.43837115 0.43837115 27 0.45399050 0.45399050 28 0.46947156 0.46947156 29 0.48480962 0.48480962 30 0.50000000 0.50000000 31 0.51503807 0.51503807 32 0.52991926 0.52991926 33 0.54463904 0.54463904 34 0.55919290 0.55919290 35 0.57357644 0.57357644 36 0.58778525 0.58778525 37 0.60181502 0.60181502 38 0.61566148 0.61566148 39 0.62932039 0.62932039 40 0.64278761 0.64278761 41 0.65605903 0.65605903 42 0.66913061 0.66913061 43 0.68199836 0.68199836 44 0.69465837 0.69465837 45 0.70710678 0.70710678 46 0.71933980 0.71933980 47 0.73135370 0.73135370 48 0.74314483 0.74314483 49 0.75470958 0.75470958 50 0.76604444 0.76604444 51 0.77714596 0.77714596 52 0.78801075 0.78801075 53 0.79863551 0.79863551 54 0.80901699 0.80901699 55 0.81915204 0.81915204 56 0.82903757 0.82903757 57 0.83867057 0.83867057 58 0.84804810 0.84804810 59 0.85716730 0.85716730 60 0.86602540 0.86602540 61 0.87461971 0.87461971 62 0.88294759 0.88294759 63 0.89100652 0.89100652 64 0.89879405 0.89879405 65 0.90630779 0.90630779 66 0.91354546 0.91354546 67 0.92050485 0.92050485 68 0.92718385 0.92718385 69 0.93358043 0.93358043 70 0.93969262 0.93969262 71 0.94551858 0.94551858 72 0.95105652 0.95105652 73 0.95630476 0.95630476 74 0.96126170 0.96126170 75 0.96592583 0.96592583 76 0.97029573 0.97029573 77 0.97437006 0.97437006 78 0.97814760 0.97814760 79 0.98162718 0.98162718 80 0.98480775 0.98480775 81 0.98768834 0.98768834 82 0.99026807 0.99026807 83 0.99254615 0.99254615 84 0.99452190 0.99452190 85 0.99619470 0.99619470 86 0.99756405 0.99756405 87 0.99862953 0.99862953 88 0.99939083 0.99939083 89 0.99984770 0.99984770 90 1.00000000 1.00000000 91 0.99984770 0.99984770 92 0.99939083 0.99939083 93 0.99862953 0.99862953 94 0.99756405 0.99756405 95 0.99619470 0.99619470 96 0.99452190 0.99452190 97 0.99254615 0.99254615 98 0.99026807 0.99026807 99 0.98768834 0.98768834 100 0.98480775 0.98480775 101 0.98162718 0.98162718 102 0.97814760 0.97814760 103 0.97437006 0.97437006 104 0.97029573 0.97029573 105 0.96592583 0.96592583 106 0.96126170 0.96126170 107 0.95630476 0.95630476 108 0.95105652 0.95105652 109 0.94551858 0.94551858 110 0.93969262 0.93969262 111 0.93358043 0.93358043 112 0.92718385 0.92718385 113 0.92050485 0.92050485 114 0.91354546 0.91354546 115 0.90630779 0.90630779 116 0.89879405 0.89879405 117 0.89100652 0.89100652 118 0.88294759 0.88294759 119 0.87461971 0.87461971 120 0.86602540 0.86602540 121 0.85716730 0.85716730 122 0.84804810 0.84804810 123 0.83867057 0.83867057 124 0.82903757 0.82903757 125 0.81915204 0.81915204 126 0.80901699 0.80901699 127 0.79863551 0.79863551 128 0.78801075 0.78801075 129 0.77714596 0.77714596 130 0.76604444 0.76604444 131 0.75470958 0.75470958 132 0.74314483 0.74314482 133 0.73135370 0.73135370 134 0.71933980 0.71933980 135 0.70710678 0.70710678 136 0.69465837 0.69465836 137 0.68199836 0.68199835 138 0.66913061 0.66913060 139 0.65605903 0.65605902 140 0.64278761 0.64278760 141 0.62932039 0.62932038 142 0.61566148 0.61566146 143 0.60181502 0.60181501 144 0.58778525 0.58778523 145 0.57357644 0.57357642 146 0.55919290 0.55919288 147 0.54463904 0.54463901 148 0.52991926 0.52991924 149 0.51503807 0.51503804 150 0.50000000 0.49999996 151 0.48480962 0.48480958 152 0.46947156 0.46947152 153 0.45399050 0.45399045 154 0.43837115 0.43837109 155 0.42261826 0.42261820 156 0.40673664 0.40673657 157 0.39073113 0.39073105 158 0.37460659 0.37460651 159 0.35836795 0.35836786 160 0.34202014 0.34202004 161 0.32556815 0.32556804 162 0.30901699 0.30901686 163 0.29237170 0.29237156 164 0.27563736 0.27563720 165 0.25881905 0.25881887 166 0.24192190 0.24192170 167 0.22495105 0.22495084 168 0.20791169 0.20791145 169 0.19080900 0.19080873 170 0.17364818 0.17364788 171 0.15643447 0.15643414 172 0.13917310 0.13917274 173 0.12186934 0.12186895 174 0.10452846 0.10452803 175 0.08715574 0.08715526 176 0.06975647 0.06975595 177 0.05233596 0.05233537 178 0.03489950 0.03489886 179 0.01745241 0.01745170 

DEPTH-based error checking

In response to the comment regarding error calculation, you examine the error associated with the Taylor-Series extensions for sin and cos databases by the number of members, changing the DEPTH and setting the maximum error EMAX 1.0e-8 using something similar to the following for range 0-360 (or 0-2PI ),

 #define DEPTH 20 #define EMAX 1.0e-8 ... /* sine as above */ ... /* cos with taylor series expansion to n = DEPTH */ long double cose (const long double deg) { long double rad = deg * M_PI / 180.0, x = 1.0; int sign = -1; for (int n = 2; n < DEPTH; n += 2, sign *= -1) x += sign * powerd (rad, n) / nfact (n); return x; } int main (void) { for (int i = 0; i < 180; i++) { long double sinlibc = sin (i * M_PI / 180.0), coslibc = cos (i * M_PI / 180.0), sints = sine (i), costs = cose (i), serr = fabs (sinlibc - sints), cerr = fabs (coslibc - costs); if (serr > EMAX) fprintf (stderr, "sine error exceeds limit of %e\n" "%3d %11.8Lf %11.8Lf %Le\n", EMAX, i, sinlibc, sints, serr); if (cerr > EMAX) fprintf (stderr, "cose error exceeds limit of %e\n" "%3d %11.8Lf %11.8Lf %Le\n", EMAX, i, coslibc, costs, cerr); } return 0; } 

If you check, you will find that for anything less than DEPTH 20 (10 members in each extension) the error will exceed 1.0e-8 for higher angles. Surprisingly, the decompositions are very accurate in the first quadrant with DEPTH values ​​below 12 (6 terms).

---

Addemdum - Improved Taylor Series Accuracy Using 0-90 and Quadrants

In the normal expansion of the Taylor series, the error increases with increasing angle. And ... because some just can't mess around, I would like to further compare the accuracy between libc sin/cos and the Taylor series if the calculations were limited to 0-90 degrees and the rest of the period from 90-360 processed by quadrant ( 2, 3 & 4 ), reflecting results from 0-90 . It works - amazingly.

For example, the results of transferring only angles 0-90 and brace angles between 90 - 180 , 180 - 270 and 270 - 360 with an initial angle % 360 give results comparable to the libc math lib functions. The maximum error between the libc and 8 and 10 Taylor-Series extensions was apparently:

Max. error from libc sin/cos

With TSLIM 16

 sine_ts max err at : 90.00 deg -- 6.023182e-12 cose_ts max err at : 270.00 deg -- 6.513370e-11 

With TSLIM 20

 sine_ts max err at : 357.00 deg -- 5.342948e-16 cose_ts max err at : 270.00 deg -- 3.557149e-15 

(with a lot of angles that don't differ at all)

The enhanced versions of sine and cose with the Taylor-Series were as follows:

 double sine (const double deg) { double fp = deg - (int64_t)deg, /* save fractional part of deg */ qdeg = (int64_t)deg % 360, /* get equivalent 0-359 deg angle */ rad, sine_deg; /* radians, sine_deg */ int pos_quad = 1, /* positive quadrant flag 1,2 */ sign = -1; /* taylor series term sign */ qdeg += fp; /* add fractional part back to angle */ /* get equivalent 0-90 degree angle, set pos_quad flag */ if (90 < qdeg && qdeg <= 180) /* in 2nd quadrant */ qdeg = 180 - qdeg; else if (180 < qdeg && qdeg <= 270) { /* in 3rd quadrant */ qdeg = qdeg - 180; pos_quad = 0; } else if (270 < qdeg && qdeg <= 360) { /* in 4th quadrant */ qdeg = 360 - qdeg; pos_quad = 0; } rad = qdeg * M_PI / 180.0; /* convert to radians */ sine_deg = rad; /* save copy for computation */ /* compute Taylor-Series expansion for sine for TSLIM / 2 terms */ for (int n = 3; n < TSLIM; n += 2, sign *= -1) { double p = rad; uint64_t f = n; for (int i = 1; i < n; i++) /* pow */ p *= rad; for (int i = 1; i < n; i++) /* nfact */ f *= i; sine_deg += sign * p / f; /* Taylor-series term */ } return pos_quad ? sine_deg : -sine_deg; } 

and for cos

 double cose (const double deg) { double fp = deg - (int64_t)deg, /* save fractional part of deg */ qdeg = (int64_t)deg % 360, /* get equivalent 0-359 deg angle */ rad, cose_deg = 1.0; /* radians, cose_deg */ int pos_quad = 1, /* positive quadrant flag 1,4 */ sign = -1; /* taylor series term sign */ qdeg += fp; /* add fractional part back to angle */ /* get equivalent 0-90 degree angle, set pos_quad flag */ if (90 < qdeg && qdeg <= 180) { /* in 2nd quadrant */ qdeg = 180 - qdeg; pos_quad = 0; } else if (180 < qdeg && qdeg <= 270) { /* in 3rd quadrant */ qdeg = qdeg - 180; pos_quad = 0; } else if (270 < qdeg && qdeg <= 360) /* in 4th quadrant */ qdeg = 360 - qdeg; rad = qdeg * M_PI / 180.0; /* convert to radians */ /* compute Taylor-Series expansion for sine for TSLIM / 2 terms */ for (int n = 2; n < TSLIM; n += 2, sign *= -1) { double p = rad; uint64_t f = n; for (int i = 1; i < n; i++) /* pow */ p *= rad; for (int i = 1; i < n; i++) /* nfact */ f *= i; cose_deg += sign * p / f; /* Taylor-series term */ } return pos_quad ? cose_deg : -cose_deg; } 

The tip of the rabbit trail was found ...

Answer 5

Changing the angle range from -90 to 90 will still cover the entire sine range. but as the Taylor series starts from scratch, the DEPTH value can be reduced to 7. As mentioned earlier, using the 64 bit unsigned fact function will fix the 67 degree problem.