1 /* Compute sine of argument.
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/>. */
21 #include <math_private.h>
22 #include <libm-alias-float.h>
23 #include "s_sincosf.h"
26 # define SINF_FUNC __sinf
28 # define SINF_FUNC SINF
36 double abstheta
= fabs (theta
);
38 if (isless (abstheta
, M_PI_4
))
40 if (abstheta
>= 0x1p
-5) /* |x| >= 2^-5. */
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
;
52 else if (abstheta
>= 0x1p
-27) /* |x| >= 2^-27. */
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
;
63 /* Handle some special cases. */
65 return theta
- (theta
* SMALL
);
70 else /* |x| >= Pi/4. */
72 unsigned int signbit
= isless (x
, 0);
73 if (isless (abstheta
, 9 * M_PI_4
)) /* |x| < 9*Pi/4. */
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
);
83 else if (isless (abstheta
, INFINITY
))
85 if (abstheta
< 0x1p
+23) /* |x| < 2^23. */
87 unsigned int n
= ((unsigned int) (abstheta
* inv_PI_4
)) + 1;
89 theta
= (abstheta
- x
* PI_2_hi
) - x
* PI_2_lo
;
90 /* Argument reduction needed. */
91 return reduced_sin (theta
, n
, signbit
);
93 else /* |x| >= 2^23. */
97 GET_FLOAT_WORD (exponent
, x
);
99 = (exponent
>> FLOAT_EXPONENT_SHIFT
) - FLOAT_EXPONENT_BIAS
;
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
;
119 return reduced_sin (e
, l
+ 1, signbit
);
129 return reduced_sin (e
, l
+ 1, signbit
);
136 return reduced_sin (e
, l
+ 1, signbit
);
144 /* High word of x. */
145 GET_FLOAT_WORD (ix
, abstheta
);
146 /* Sin(Inf or NaN) is NaN. */
147 if (ix
== 0x7f800000)
155 libm_alias_float (__sin
, sin
)