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