Sync with trunk (48237)
[reactos.git] / lib / sdk / crt / math / sin.c
1 /*
2 * COPYRIGHT: See COPYING in the top level directory
3 * PROJECT: ReactOS CRT
4 * FILE: lib/crt/math/sin.c
5 * PURPOSE: Generic C Implementation of sin
6 * PROGRAMMER: Timo Kreuzer (timo.kreuzer@reactos.org)
7 */
8
9 #define PRECISION 9
10 #define M_PI 3.141592653589793238462643
11
12 static double sin_off_tbl[] = {0.0, -M_PI/2., 0, -M_PI/2.};
13 static double sin_sign_tbl[] = {1,-1,-1,1};
14
15 double
16 sin(double x)
17 {
18 int quadrant;
19 double x2, result;
20
21 /* Calculate the quadrant */
22 quadrant = x * (2./M_PI);
23
24 /* Get offset inside quadrant */
25 x = x - quadrant * (M_PI/2.);
26
27 /* Normalize quadrant to [0..3] */
28 quadrant = (quadrant - 1) & 0x3;
29
30 /* Fixup value for the generic function */
31 x += sin_off_tbl[quadrant];
32
33 /* Calculate the negative of the square of x */
34 x2 = - (x * x);
35
36 /* This is an unrolled taylor series using <PRECISION> iterations
37 * Example with 4 iterations:
38 * result = 1 - x^2/2! + x^4/4! - x^6/6! + x^8/8!
39 * To save multiplications and to keep the precision high, it's performed
40 * like this:
41 * result = 1 - x^2 * (1/2! - x^2 * (1/4! - x^2 * (1/6! - x^2 * (1/8!))))
42 */
43
44 /* Start with 0, compiler will optimize this away */
45 result = 0;
46
47 #if (PRECISION >= 10)
48 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20);
49 result *= x2;
50 #endif
51 #if (PRECISION >= 9)
52 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18);
53 result *= x2;
54 #endif
55 #if (PRECISION >= 8)
56 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16);
57 result *= x2;
58 #endif
59 #if (PRECISION >= 7)
60 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14);
61 result *= x2;
62 #endif
63 #if (PRECISION >= 6)
64 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12);
65 result *= x2;
66 #endif
67 #if (PRECISION >= 5)
68 result += 1./(1.*2*3*4*5*6*7*8*9*10);
69 result *= x2;
70 #endif
71 result += 1./(1.*2*3*4*5*6*7*8);
72 result *= x2;
73
74 result += 1./(1.*2*3*4*5*6);
75 result *= x2;
76
77 result += 1./(1.*2*3*4);
78 result *= x2;
79
80 result += 1./(1.*2);
81 result *= x2;
82
83 result += 1;
84
85 /* Apply correct sign */
86 result *= sin_sign_tbl[quadrant];
87
88 return result;
89 }