1 | /* Used by sinf, cosf and sincosf functions. |
2 | Copyright (C) 2017-2018 Free Software Foundation, Inc. |
3 | This file is part of the GNU C Library. |
4 | |
5 | The GNU C Library is free software; you can redistribute it and/or |
6 | modify it under the terms of the GNU Lesser General Public |
7 | License as published by the Free Software Foundation; either |
8 | version 2.1 of the License, or (at your option) any later version. |
9 | |
10 | The GNU C Library is distributed in the hope that it will be useful, |
11 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
13 | Lesser General Public License for more details. |
14 | |
15 | You should have received a copy of the GNU Lesser General Public |
16 | License along with the GNU C Library; if not, see |
17 | <http://www.gnu.org/licenses/>. */ |
18 | |
19 | /* Chebyshev constants for cos, range -PI/4 - PI/4. */ |
20 | static const double C0 = -0x1.ffffffffe98aep-2; |
21 | static const double C1 = 0x1.55555545c50c7p-5; |
22 | static const double C2 = -0x1.6c16b348b6874p-10; |
23 | static const double C3 = 0x1.a00eb9ac43ccp-16; |
24 | static const double C4 = -0x1.23c97dd8844d7p-22; |
25 | |
26 | /* Chebyshev constants for sin, range -PI/4 - PI/4. */ |
27 | static const double S0 = -0x1.5555555551cd9p-3; |
28 | static const double S1 = 0x1.1111110c2688bp-7; |
29 | static const double S2 = -0x1.a019f8b4bd1f9p-13; |
30 | static const double S3 = 0x1.71d7264e6b5b4p-19; |
31 | static const double S4 = -0x1.a947e1674b58ap-26; |
32 | |
33 | /* Chebyshev constants for sin, range 2^-27 - 2^-5. */ |
34 | static const double SS0 = -0x1.555555543d49dp-3; |
35 | static const double SS1 = 0x1.110f475cec8c5p-7; |
36 | |
37 | /* Chebyshev constants for cos, range 2^-27 - 2^-5. */ |
38 | static const double CC0 = -0x1.fffffff5cc6fdp-2; |
39 | static const double CC1 = 0x1.55514b178dac5p-5; |
40 | |
41 | /* PI/2 with 98 bits of accuracy. */ |
42 | static const double PI_2_hi = 0x1.921fb544p+0; |
43 | static const double PI_2_lo = 0x1.0b4611a626332p-34; |
44 | |
45 | static const double SMALL = 0x1p-50; /* 2^-50. */ |
46 | static const double inv_PI_4 = 0x1.45f306dc9c883p+0; /* 4/PI. */ |
47 | |
48 | #define FLOAT_EXPONENT_SHIFT 23 |
49 | #define FLOAT_EXPONENT_BIAS 127 |
50 | |
51 | static const double pio2_table[] = { |
52 | 0 * M_PI_2, |
53 | 1 * M_PI_2, |
54 | 2 * M_PI_2, |
55 | 3 * M_PI_2, |
56 | 4 * M_PI_2, |
57 | 5 * M_PI_2 |
58 | }; |
59 | |
60 | static const double invpio4_table[] = { |
61 | 0x0p+0, |
62 | 0x1.45f306cp+0, |
63 | 0x1.c9c882ap-28, |
64 | 0x1.4fe13a8p-58, |
65 | 0x1.f47d4dp-85, |
66 | 0x1.bb81b6cp-112, |
67 | 0x1.4acc9ep-142, |
68 | 0x1.0e4107cp-169 |
69 | }; |
70 | |
71 | static const double ones[] = { 1.0, -1.0 }; |
72 | |
73 | /* Compute the sine value using Chebyshev polynomials where |
74 | THETA is the range reduced absolute value of the input |
75 | and it is less than Pi/4, |
76 | N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide |
77 | whether a sine or cosine approximation is more accurate and |
78 | SIGNBIT is used to add the correct sign after the Chebyshev |
79 | polynomial is computed. */ |
80 | static inline float |
81 | reduced_sin (const double theta, const unsigned int n, |
82 | const unsigned int signbit) |
83 | { |
84 | double sx; |
85 | const double theta2 = theta * theta; |
86 | /* We are operating on |x|, so we need to add back the original |
87 | signbit for sinf. */ |
88 | double sign; |
89 | /* Determine positive or negative primary interval. */ |
90 | sign = ones[((n >> 2) & 1) ^ signbit]; |
91 | /* Are we in the primary interval of sin or cos? */ |
92 | if ((n & 2) == 0) |
93 | { |
94 | /* Here sinf() is calculated using sin Chebyshev polynomial: |
95 | x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */ |
96 | sx = S3 + theta2 * S4; /* S3+x^2*S4. */ |
97 | sx = S2 + theta2 * sx; /* S2+x^2*(S3+x^2*S4). */ |
98 | sx = S1 + theta2 * sx; /* S1+x^2*(S2+x^2*(S3+x^2*S4)). */ |
99 | sx = S0 + theta2 * sx; /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))). */ |
100 | sx = theta + theta * theta2 * sx; |
101 | } |
102 | else |
103 | { |
104 | /* Here sinf() is calculated using cos Chebyshev polynomial: |
105 | 1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))). */ |
106 | sx = C3 + theta2 * C4; /* C3+x^2*C4. */ |
107 | sx = C2 + theta2 * sx; /* C2+x^2*(C3+x^2*C4). */ |
108 | sx = C1 + theta2 * sx; /* C1+x^2*(C2+x^2*(C3+x^2*C4)). */ |
109 | sx = C0 + theta2 * sx; /* C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4))). */ |
110 | sx = 1.0 + theta2 * sx; |
111 | } |
112 | |
113 | /* Add in the signbit and assign the result. */ |
114 | return sign * sx; |
115 | } |
116 | |
117 | /* Compute the cosine value using Chebyshev polynomials where |
118 | THETA is the range reduced absolute value of the input |
119 | and it is less than Pi/4, |
120 | N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide |
121 | whether a sine or cosine approximation is more accurate and |
122 | the sign of the result. */ |
123 | static inline float |
124 | reduced_cos (double theta, unsigned int n) |
125 | { |
126 | double sign, cx; |
127 | const double theta2 = theta * theta; |
128 | |
129 | /* Determine positive or negative primary interval. */ |
130 | n += 2; |
131 | sign = ones[(n >> 2) & 1]; |
132 | |
133 | /* Are we in the primary interval of sin or cos? */ |
134 | if ((n & 2) == 0) |
135 | { |
136 | /* Here cosf() is calculated using sin Chebyshev polynomial: |
137 | x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */ |
138 | cx = S3 + theta2 * S4; |
139 | cx = S2 + theta2 * cx; |
140 | cx = S1 + theta2 * cx; |
141 | cx = S0 + theta2 * cx; |
142 | cx = theta + theta * theta2 * cx; |
143 | } |
144 | else |
145 | { |
146 | /* Here cosf() is calculated using cos Chebyshev polynomial: |
147 | 1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))). */ |
148 | cx = C3 + theta2 * C4; |
149 | cx = C2 + theta2 * cx; |
150 | cx = C1 + theta2 * cx; |
151 | cx = C0 + theta2 * cx; |
152 | cx = 1. + theta2 * cx; |
153 | } |
154 | return sign * cx; |
155 | } |
156 | |