1 /*
2 * Copyright (c) 2006-2018, RT-Thread Development Team
3 *
4 * SPDX-License-Identifier: Apache-2.0
5 *
6 * Change Logs:
7 * Date Author Notes
8 */
9 #include <math.h>
10
11 /*
12 * COPYRIGHT: See COPYING in the top level directory
13 * PROJECT: ReactOS CRT
14 * FILE: lib/crt/math/cos.c
15 * PURPOSE: Generic C Implementation of cos
16 * PROGRAMMER: Timo Kreuzer ([email protected])
17 */
18
19 #define PRECISION 9
20
21 static double cos_off_tbl[] = {0.0, -M_PI/2., 0, -M_PI/2.};
22 static double cos_sign_tbl[] = {1,-1,-1,1};
23
24 static double sin_off_tbl[] = {0.0, -M_PI/2., 0, -M_PI/2.};
25 static double sin_sign_tbl[] = {1,-1,-1,1};
26
sin(double x)27 double sin(double x)
28 {
29 int quadrant;
30 double x2, result;
31
32 /* Calculate the quadrant */
33 quadrant = x * (2./M_PI);
34
35 /* Get offset inside quadrant */
36 x = x - quadrant * (M_PI/2.);
37
38 /* Normalize quadrant to [0..3] */
39 quadrant = (quadrant - 1) & 0x3;
40
41 /* Fixup value for the generic function */
42 x += sin_off_tbl[quadrant];
43
44 /* Calculate the negative of the square of x */
45 x2 = - (x * x);
46
47 /* This is an unrolled taylor series using <PRECISION> iterations
48 * Example with 4 iterations:
49 * result = 1 - x^2/2! + x^4/4! - x^6/6! + x^8/8!
50 * To save multiplications and to keep the precision high, it's performed
51 * like this:
52 * result = 1 - x^2 * (1/2! - x^2 * (1/4! - x^2 * (1/6! - x^2 * (1/8!))))
53 */
54
55 /* Start with 0, compiler will optimize this away */
56 result = 0;
57
58 #if (PRECISION >= 10)
59 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20);
60 result *= x2;
61 #endif
62 #if (PRECISION >= 9)
63 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18);
64 result *= x2;
65 #endif
66 #if (PRECISION >= 8)
67 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16);
68 result *= x2;
69 #endif
70 #if (PRECISION >= 7)
71 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14);
72 result *= x2;
73 #endif
74 #if (PRECISION >= 6)
75 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12);
76 result *= x2;
77 #endif
78 #if (PRECISION >= 5)
79 result += 1./(1.*2*3*4*5*6*7*8*9*10);
80 result *= x2;
81 #endif
82 result += 1./(1.*2*3*4*5*6*7*8);
83 result *= x2;
84
85 result += 1./(1.*2*3*4*5*6);
86 result *= x2;
87
88 result += 1./(1.*2*3*4);
89 result *= x2;
90
91 result += 1./(1.*2);
92 result *= x2;
93
94 result += 1;
95
96 /* Apply correct sign */
97 result *= sin_sign_tbl[quadrant];
98
99 return result;
100 }
101
cos(double x)102 double cos(double x)
103 {
104 int quadrant;
105 double x2, result;
106
107 /* Calculate the quadrant */
108 quadrant = x * (2./M_PI);
109
110 /* Get offset inside quadrant */
111 x = x - quadrant * (M_PI/2.);
112
113 /* Normalize quadrant to [0..3] */
114 quadrant = quadrant & 0x3;
115
116 /* Fixup value for the generic function */
117 x += cos_off_tbl[quadrant];
118
119 /* Calculate the negative of the square of x */
120 x2 = - (x * x);
121
122 /* This is an unrolled taylor series using <PRECISION> iterations
123 * Example with 4 iterations:
124 * result = 1 - x^2/2! + x^4/4! - x^6/6! + x^8/8!
125 * To save multiplications and to keep the precision high, it's performed
126 * like this:
127 * result = 1 - x^2 * (1/2! - x^2 * (1/4! - x^2 * (1/6! - x^2 * (1/8!))))
128 */
129
130 /* Start with 0, compiler will optimize this away */
131 result = 0;
132
133 #if (PRECISION >= 10)
134 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20);
135 result *= x2;
136 #endif
137 #if (PRECISION >= 9)
138 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18);
139 result *= x2;
140 #endif
141 #if (PRECISION >= 8)
142 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16);
143 result *= x2;
144 #endif
145 #if (PRECISION >= 7)
146 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12*13*14);
147 result *= x2;
148 #endif
149 #if (PRECISION >= 6)
150 result += 1./(1.*2*3*4*5*6*7*8*9*10*11*12);
151 result *= x2;
152 #endif
153 #if (PRECISION >= 5)
154 result += 1./(1.*2*3*4*5*6*7*8*9*10);
155 result *= x2;
156 #endif
157 result += 1./(1.*2*3*4*5*6*7*8);
158 result *= x2;
159
160 result += 1./(1.*2*3*4*5*6);
161 result *= x2;
162
163 result += 1./(1.*2*3*4);
164 result *= x2;
165
166 result += 1./(1.*2);
167 result *= x2;
168
169 result += 1;
170
171 /* Apply correct sign */
172 result *= cos_sign_tbl[quadrant];
173
174 return result;
175 }
176
177