2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2019 Free Software Foundation, Inc.
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <http://www.gnu.org/licenses/>.
19 /****************************************************************************/
21 /* MODULE_NAME:usncs.c */
25 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */
26 /* branred.c sincos.tbl */
28 /* An ultimate sin and cos routine. Given an IEEE double machine number x */
29 /* it computes sin(x) or cos(x) with ~0.55 ULP. */
30 /* Assumption: Machine arithmetic operations are performed in */
31 /* round to nearest mode of IEEE 754 standard. */
33 /****************************************************************************/
43 #include <math_private.h>
44 #include <fenv_private.h>
45 #include <math-underflow.h>
46 #include <libm-alias-double.h>
49 /* Helper macros to compute sin of the input values. */
50 #define POLYNOMIAL2(xx) ((((s5 * (xx) + s4) * (xx) + s3) * (xx) + s2) * (xx))
52 #define POLYNOMIAL(xx) (POLYNOMIAL2 (xx) + s1)
54 /* The computed polynomial is a variation of the Taylor series expansion for
57 a - a^3/3! + a^5/5! - a^7/7! + a^9/9! + (1 - a^2) * da / 2
59 The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so
60 on. The result is returned to LHS. */
61 #define TAYLOR_SIN(xx, a, da) \
63 double t = ((POLYNOMIAL (xx) * (a) - 0.5 * (da)) * (xx) + (da)); \
64 double res = (a) + t; \
68 #define SINCOS_TABLE_LOOKUP(u, sn, ssn, cs, ccs) \
70 int4 k = u.i[LOW_HALF] << 2; \
71 sn = __sincostab.x[k]; \
72 ssn = __sincostab.x[k + 1]; \
73 cs = __sincostab.x[k + 2]; \
74 ccs = __sincostab.x[k + 3]; \
85 } __sincostab attribute_hidden
;
88 sn3
= -1.66666666666664880952546298448555E-01,
89 sn5
= 8.33333214285722277379541354343671E-03,
90 cs2
= 4.99999999999999999999950396842453E-01,
91 cs4
= -4.16666666666664434524222570944589E-02,
92 cs6
= 1.38888874007937613028114285595617E-03;
94 int __branred (double x
, double *a
, double *aa
);
96 /* Given a number partitioned into X and DX, this function computes the cosine
97 of the number by combining the sin and cos of X (as computed by a variation
98 of the Taylor series) with the values looked up from the sin/cos table to
102 do_cos (double x
, double dx
)
109 u
.x
= big
+ fabs (x
);
110 x
= fabs (x
) - (u
.x
- big
) + dx
;
112 double xx
, s
, sn
, ssn
, c
, cs
, ccs
, cor
;
114 s
= x
+ x
* xx
* (sn3
+ xx
* sn5
);
115 c
= xx
* (cs2
+ xx
* (cs4
+ xx
* cs6
));
116 SINCOS_TABLE_LOOKUP (u
, sn
, ssn
, cs
, ccs
);
117 cor
= (ccs
- s
* ssn
- cs
* c
) - sn
* s
;
121 /* Given a number partitioned into X and DX, this function computes the sine of
122 the number by combining the sin and cos of X (as computed by a variation of
123 the Taylor series) with the values looked up from the sin/cos table to get
127 do_sin (double x
, double dx
)
130 /* Max ULP is 0.501 if |x| < 0.126, otherwise ULP is 0.518. */
131 if (fabs (x
) < 0.126)
132 return TAYLOR_SIN (x
* x
, x
, dx
);
138 u
.x
= big
+ fabs (x
);
139 x
= fabs (x
) - (u
.x
- big
);
141 double xx
, s
, sn
, ssn
, c
, cs
, ccs
, cor
;
143 s
= x
+ (dx
+ x
* xx
* (sn3
+ xx
* sn5
));
144 c
= x
* dx
+ xx
* (cs2
+ xx
* (cs4
+ xx
* cs6
));
145 SINCOS_TABLE_LOOKUP (u
, sn
, ssn
, cs
, ccs
);
146 cor
= (ssn
+ s
* ccs
- sn
* c
) + cs
* s
;
147 return copysign (sn
+ cor
, xold
);
150 /* Reduce range of x to within PI/2 with abs (x) < 105414350. The high part
151 is written to *a, the low part to *da. Range reduction is accurate to 136
152 bits so that when x is large and *a very close to zero, all 53 bits of *a
156 reduce_sincos (double x
, double *a
, double *da
)
160 double t
= (x
* hpinv
+ toint
);
161 double xn
= t
- toint
;
163 double y
= (x
- xn
* mp1
) - xn
* mp2
;
164 int4 n
= v
.i
[LOW_HALF
] & 3;
166 double b
, db
, t1
, t2
;
180 /* Compute sin or cos (A + DA) for the given quadrant N. */
183 do_sincos (double a
, double da
, int4 n
)
188 /* Max ULP is 0.513. */
189 retval
= do_cos (a
, da
);
191 /* Max ULP is 0.501 if xx < 0.01588, otherwise ULP is 0.518. */
192 retval
= do_sin (a
, da
);
194 return (n
& 2) ? -retval
: retval
;
198 /*******************************************************************/
199 /* An ultimate sin routine. Given an IEEE double machine number x */
200 /* it computes the correctly rounded (to nearest) value of sin(x) */
201 /*******************************************************************/
212 SET_RESTORE_ROUND_53BIT (FE_TONEAREST
);
216 k
= 0x7fffffff & m
; /* no sign */
217 if (k
< 0x3e500000) /* if x->0 =>sin(x)=x */
219 math_check_force_underflow (x
);
222 /*--------------------------- 2^-26<|x|< 0.855469---------------------- */
223 else if (k
< 0x3feb6000)
225 /* Max ULP is 0.548. */
226 retval
= do_sin (x
, 0);
227 } /* else if (k < 0x3feb6000) */
229 /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/
230 else if (k
< 0x400368fd)
233 /* Max ULP is 0.51. */
234 retval
= copysign (do_cos (t
, hp1
), x
);
235 } /* else if (k < 0x400368fd) */
237 /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
238 else if (k
< 0x419921FB)
240 n
= reduce_sincos (x
, &a
, &da
);
241 retval
= do_sincos (a
, da
, n
);
242 } /* else if (k < 0x419921FB ) */
244 /* --------------------105414350 <|x| <2^1024------------------------------*/
245 else if (k
< 0x7ff00000)
247 n
= __branred (x
, &a
, &da
);
248 retval
= do_sincos (a
, da
, n
);
250 /*--------------------- |x| > 2^1024 ----------------------------------*/
253 if (k
== 0x7ff00000 && u
.i
[LOW_HALF
] == 0)
262 /*******************************************************************/
263 /* An ultimate cos routine. Given an IEEE double machine number x */
264 /* it computes the correctly rounded (to nearest) value of cos(x) */
265 /*******************************************************************/
277 SET_RESTORE_ROUND_53BIT (FE_TONEAREST
);
283 /* |x|<2^-27 => cos(x)=1 */
287 else if (k
< 0x3feb6000)
288 { /* 2^-27 < |x| < 0.855469 */
289 /* Max ULP is 0.51. */
290 retval
= do_cos (x
, 0);
291 } /* else if (k < 0x3feb6000) */
293 else if (k
< 0x400368fd)
294 { /* 0.855469 <|x|<2.426265 */ ;
298 /* Max ULP is 0.501 if xx < 0.01588 or 0.518 otherwise.
299 Range reduction uses 106 bits here which is sufficient. */
300 retval
= do_sin (a
, da
);
301 } /* else if (k < 0x400368fd) */
303 else if (k
< 0x419921FB)
304 { /* 2.426265<|x|< 105414350 */
305 n
= reduce_sincos (x
, &a
, &da
);
306 retval
= do_sincos (a
, da
, n
+ 1);
307 } /* else if (k < 0x419921FB ) */
309 /* 105414350 <|x| <2^1024 */
310 else if (k
< 0x7ff00000)
312 n
= __branred (x
, &a
, &da
);
313 retval
= do_sincos (a
, da
, n
+ 1);
318 if (k
== 0x7ff00000 && u
.i
[LOW_HALF
] == 0)
320 retval
= x
/ x
; /* |x| > 2^1024 */
327 libm_alias_double (__cos
, cos
)
330 libm_alias_double (__sin
, sin
)