1 | /* Compute sine of argument. |
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 | #include <errno.h> |
20 | #include <math.h> |
21 | #include <math_private.h> |
22 | #include <libm-alias-float.h> |
23 | #include "s_sincosf.h" |
24 | |
25 | #ifndef SINF |
26 | # define SINF_FUNC __sinf |
27 | #else |
28 | # define SINF_FUNC SINF |
29 | #endif |
30 | |
31 | float |
32 | SINF_FUNC (float x) |
33 | { |
34 | double cx; |
35 | double theta = x; |
36 | double abstheta = fabs (theta); |
37 | /* If |x|< Pi/4. */ |
38 | if (isless (abstheta, M_PI_4)) |
39 | { |
40 | if (abstheta >= 0x1p-5) /* |x| >= 2^-5. */ |
41 | { |
42 | const double theta2 = theta * theta; |
43 | /* Chebyshev polynomial of the form for sin |
44 | x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */ |
45 | cx = S3 + theta2 * S4; |
46 | cx = S2 + theta2 * cx; |
47 | cx = S1 + theta2 * cx; |
48 | cx = S0 + theta2 * cx; |
49 | cx = theta + theta * theta2 * cx; |
50 | return cx; |
51 | } |
52 | else if (abstheta >= 0x1p-27) /* |x| >= 2^-27. */ |
53 | { |
54 | /* A simpler Chebyshev approximation is close enough for this range: |
55 | for sin: x+x^3*(SS0+x^2*SS1). */ |
56 | const double theta2 = theta * theta; |
57 | cx = SS0 + theta2 * SS1; |
58 | cx = theta + theta * theta2 * cx; |
59 | return cx; |
60 | } |
61 | else |
62 | { |
63 | /* Handle some special cases. */ |
64 | if (theta) |
65 | return theta - (theta * SMALL); |
66 | else |
67 | return theta; |
68 | } |
69 | } |
70 | else /* |x| >= Pi/4. */ |
71 | { |
72 | unsigned int signbit = isless (x, 0); |
73 | if (isless (abstheta, 9 * M_PI_4)) /* |x| < 9*Pi/4. */ |
74 | { |
75 | /* There are cases where FE_UPWARD rounding mode can |
76 | produce a result of abstheta * inv_PI_4 == 9, |
77 | where abstheta < 9pi/4, so the domain for |
78 | pio2_table must go to 5 (9 / 2 + 1). */ |
79 | unsigned int n = (abstheta * inv_PI_4) + 1; |
80 | theta = abstheta - pio2_table[n / 2]; |
81 | return reduced_sin (theta, n, signbit); |
82 | } |
83 | else if (isless (abstheta, INFINITY)) |
84 | { |
85 | if (abstheta < 0x1p+23) /* |x| < 2^23. */ |
86 | { |
87 | unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1; |
88 | double x = n / 2; |
89 | theta = (abstheta - x * PI_2_hi) - x * PI_2_lo; |
90 | /* Argument reduction needed. */ |
91 | return reduced_sin (theta, n, signbit); |
92 | } |
93 | else /* |x| >= 2^23. */ |
94 | { |
95 | x = fabsf (x); |
96 | int exponent; |
97 | GET_FLOAT_WORD (exponent, x); |
98 | exponent |
99 | = (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS; |
100 | exponent += 3; |
101 | exponent /= 28; |
102 | double a = invpio4_table[exponent] * x; |
103 | double b = invpio4_table[exponent + 1] * x; |
104 | double c = invpio4_table[exponent + 2] * x; |
105 | double d = invpio4_table[exponent + 3] * x; |
106 | uint64_t l = a; |
107 | l &= ~0x7; |
108 | a -= l; |
109 | double e = a + b; |
110 | l = e; |
111 | e = a - l; |
112 | if (l & 1) |
113 | { |
114 | e -= 1.0; |
115 | e += b; |
116 | e += c; |
117 | e += d; |
118 | e *= M_PI_4; |
119 | return reduced_sin (e, l + 1, signbit); |
120 | } |
121 | else |
122 | { |
123 | e += b; |
124 | e += c; |
125 | e += d; |
126 | if (e <= 1.0) |
127 | { |
128 | e *= M_PI_4; |
129 | return reduced_sin (e, l + 1, signbit); |
130 | } |
131 | else |
132 | { |
133 | l++; |
134 | e -= 2.0; |
135 | e *= M_PI_4; |
136 | return reduced_sin (e, l + 1, signbit); |
137 | } |
138 | } |
139 | } |
140 | } |
141 | else |
142 | { |
143 | int32_t ix; |
144 | /* High word of x. */ |
145 | GET_FLOAT_WORD (ix, abstheta); |
146 | /* Sin(Inf or NaN) is NaN. */ |
147 | if (ix == 0x7f800000) |
148 | __set_errno (EDOM); |
149 | return x - x; |
150 | } |
151 | } |
152 | } |
153 | |
154 | #ifndef SINF |
155 | libm_alias_float (__sin, sin) |
156 | #endif |
157 | |