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
31float
32SINF_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
155libm_alias_float (__sin, sin)
156#endif
157