]>
Commit | Line | Data |
---|---|---|
7799b7b3 | 1 | /* Compute sine and cosine of argument. |
688903eb | 2 | Copyright (C) 2017-2018 Free Software Foundation, Inc. |
7799b7b3 | 3 | This file is part of the GNU C Library. |
7799b7b3 UD |
4 | |
5 | The GNU C Library is free software; you can redistribute it and/or | |
41bdb6e2 AJ |
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. | |
7799b7b3 UD |
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 | |
41bdb6e2 | 13 | Lesser General Public License for more details. |
7799b7b3 | 14 | |
41bdb6e2 | 15 | You should have received a copy of the GNU Lesser General Public |
59ba27a6 PE |
16 | License along with the GNU C Library; if not, see |
17 | <http://www.gnu.org/licenses/>. */ | |
7799b7b3 | 18 | |
d435569c | 19 | #include <errno.h> |
7799b7b3 | 20 | #include <math.h> |
1ed0291c | 21 | #include <math_private.h> |
2f49ce7d | 22 | #include <libm-alias-float.h> |
984ae996 | 23 | #include "s_sincosf.h" |
7799b7b3 | 24 | |
22bf5c17 LD |
25 | #ifndef SINCOSF |
26 | # define SINCOSF_FUNC __sincosf | |
27 | #else | |
28 | # define SINCOSF_FUNC SINCOSF | |
29 | #endif | |
7799b7b3 UD |
30 | |
31 | void | |
22bf5c17 | 32 | SINCOSF_FUNC (float x, float *sinx, float *cosx) |
7799b7b3 | 33 | { |
984ae996 RS |
34 | double cx; |
35 | double theta = x; | |
36 | double abstheta = fabs (theta); | |
37 | /* If |x|< Pi/4. */ | |
38 | if (isless (abstheta, M_PI_4)) | |
7799b7b3 | 39 | { |
984ae996 RS |
40 | if (abstheta >= 0x1p-5) /* |x| >= 2^-5. */ |
41 | { | |
42 | const double theta2 = theta * theta; | |
43 | /* Chebyshev polynomial of the form for sin and cos. */ | |
44 | cx = C3 + theta2 * C4; | |
45 | cx = C2 + theta2 * cx; | |
46 | cx = C1 + theta2 * cx; | |
47 | cx = C0 + theta2 * cx; | |
48 | cx = 1.0 + theta2 * cx; | |
49 | *cosx = cx; | |
50 | cx = S3 + theta2 * S4; | |
51 | cx = S2 + theta2 * cx; | |
52 | cx = S1 + theta2 * cx; | |
53 | cx = S0 + theta2 * cx; | |
54 | cx = theta + theta * theta2 * cx; | |
55 | *sinx = cx; | |
56 | } | |
57 | else if (abstheta >= 0x1p-27) /* |x| >= 2^-27. */ | |
58 | { | |
59 | /* A simpler Chebyshev approximation is close enough for this range: | |
60 | for sin: x+x^3*(SS0+x^2*SS1) | |
61 | for cos: 1.0+x^2*(CC0+x^3*CC1). */ | |
62 | const double theta2 = theta * theta; | |
63 | cx = CC0 + theta * theta2 * CC1; | |
64 | cx = 1.0 + theta2 * cx; | |
65 | *cosx = cx; | |
66 | cx = SS0 + theta2 * SS1; | |
67 | cx = theta + theta * theta2 * cx; | |
68 | *sinx = cx; | |
69 | } | |
70 | else | |
71 | { | |
72 | /* Handle some special cases. */ | |
73 | if (theta) | |
74 | *sinx = theta - (theta * SMALL); | |
75 | else | |
76 | *sinx = theta; | |
77 | *cosx = 1.0 - abstheta; | |
78 | } | |
7799b7b3 | 79 | } |
984ae996 | 80 | else /* |x| >= Pi/4. */ |
7799b7b3 | 81 | { |
984ae996 RS |
82 | unsigned int signbit = isless (x, 0); |
83 | if (isless (abstheta, 9 * M_PI_4)) /* |x| < 9*Pi/4. */ | |
84 | { | |
85 | /* There are cases where FE_UPWARD rounding mode can | |
86 | produce a result of abstheta * inv_PI_4 == 9, | |
87 | where abstheta < 9pi/4, so the domain for | |
88 | pio2_table must go to 5 (9 / 2 + 1). */ | |
89 | unsigned int n = (abstheta * inv_PI_4) + 1; | |
90 | theta = abstheta - pio2_table[n / 2]; | |
91 | *sinx = reduced_sin (theta, n, signbit); | |
92 | *cosx = reduced_cos (theta, n); | |
93 | } | |
94 | else if (isless (abstheta, INFINITY)) | |
95 | { | |
96 | if (abstheta < 0x1p+23) /* |x| < 2^23. */ | |
97 | { | |
98 | unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1; | |
99 | double x = n / 2; | |
100 | theta = (abstheta - x * PI_2_hi) - x * PI_2_lo; | |
101 | /* Argument reduction needed. */ | |
102 | *sinx = reduced_sin (theta, n, signbit); | |
103 | *cosx = reduced_cos (theta, n); | |
104 | } | |
105 | else /* |x| >= 2^23. */ | |
106 | { | |
107 | x = fabsf (x); | |
108 | int exponent; | |
109 | GET_FLOAT_WORD (exponent, x); | |
110 | exponent | |
111 | = (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS; | |
112 | exponent += 3; | |
113 | exponent /= 28; | |
114 | double a = invpio4_table[exponent] * x; | |
115 | double b = invpio4_table[exponent + 1] * x; | |
116 | double c = invpio4_table[exponent + 2] * x; | |
117 | double d = invpio4_table[exponent + 3] * x; | |
118 | uint64_t l = a; | |
119 | l &= ~0x7; | |
120 | a -= l; | |
121 | double e = a + b; | |
122 | l = e; | |
123 | e = a - l; | |
124 | if (l & 1) | |
125 | { | |
126 | e -= 1.0; | |
127 | e += b; | |
128 | e += c; | |
129 | e += d; | |
130 | e *= M_PI_4; | |
131 | *sinx = reduced_sin (e, l + 1, signbit); | |
132 | *cosx = reduced_cos (e, l + 1); | |
133 | } | |
134 | else | |
135 | { | |
136 | e += b; | |
137 | e += c; | |
138 | e += d; | |
139 | if (e <= 1.0) | |
140 | { | |
141 | e *= M_PI_4; | |
142 | *sinx = reduced_sin (e, l + 1, signbit); | |
143 | *cosx = reduced_cos (e, l + 1); | |
144 | } | |
145 | else | |
146 | { | |
147 | l++; | |
148 | e -= 2.0; | |
149 | e *= M_PI_4; | |
150 | *sinx = reduced_sin (e, l + 1, signbit); | |
151 | *cosx = reduced_cos (e, l + 1); | |
152 | } | |
153 | } | |
154 | } | |
155 | } | |
156 | else | |
7799b7b3 | 157 | { |
984ae996 RS |
158 | int32_t ix; |
159 | /* High word of x. */ | |
160 | GET_FLOAT_WORD (ix, abstheta); | |
161 | /* sin/cos(Inf or NaN) is NaN. */ | |
162 | *sinx = *cosx = x - x; | |
163 | if (ix == 0x7f800000) | |
164 | __set_errno (EDOM); | |
7799b7b3 UD |
165 | } |
166 | } | |
167 | } | |
22bf5c17 LD |
168 | |
169 | #ifndef SINCOSF | |
2f49ce7d | 170 | libm_alias_float (__sincos, sincos) |
22bf5c17 | 171 | #endif |