1 /* Used by sinf, cosf and sincosf functions.
2 Copyright (C) 2017 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
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.
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.
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/>. */
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;
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;
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;
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;
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;
45 static const double SMALL
= 0x1p
-50; /* 2^-50. */
46 static const double inv_PI_4
= 0x1.45f306dc9c883p
+0; /* 4/PI. */
48 #define FLOAT_EXPONENT_SHIFT 23
49 #define FLOAT_EXPONENT_BIAS 127
51 static const double pio2_table
[] = {
60 static const double invpio4_table
[] = {
71 static const double ones
[] = { 1.0, -1.0 };
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. */
81 reduced_sin (const double theta
, const unsigned int n
,
82 const unsigned int signbit
)
85 const double theta2
= theta
* theta
;
86 /* We are operating on |x|, so we need to add back the original
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? */
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
;
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
;
113 /* Add in the signbit and assign the result. */
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. */
124 reduced_cos (double theta
, unsigned int n
)
127 const double theta2
= theta
* theta
;
129 /* Determine positive or negative primary interval. */
131 sign
= ones
[(n
>> 2) & 1];
133 /* Are we in the primary interval of sin or cos? */
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
;
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
;