]> git.ipfire.org Git - thirdparty/glibc.git/blame - sysdeps/ieee754/dbl-64/s_sin.c
projects / thirdparty / glibc.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Consolidate definition of constant t22
[thirdparty/glibc.git] / sysdeps / ieee754 / dbl-64 / s_sin.c
CommitLineData
f7eac6eb 1/*
e4d82761 2 * IBM Accurate Mathematical Library
aeb25823 3 * written by International Business Machines Corp.
568035b7 4 * Copyright (C) 2001-2013 Free Software Foundation, Inc.
f7eac6eb 5 *
e4d82761
UD		 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
cc7375ce 8 * the Free Software Foundation; either version 2.1 of the License, or
e4d82761 9 * (at your option) any later version.
f7eac6eb 10 *
e4d82761
UD		 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
c6c6dd48 14 * GNU Lesser General Public License for more details.
f7eac6eb 15 *
e4d82761 16 * You should have received a copy of the GNU Lesser General Public License
59ba27a6 17 * along with this program; if not, see <http://www.gnu.org/licenses/>.
f7eac6eb 18 */
e4d82761
UD		 19/****************************************************************************/
20/* */
21/* MODULE_NAME:usncs.c */
22/* */
23/* FUNCTIONS: usin */
24/* ucos */
25/* slow */
26/* slow1 */
27/* slow2 */
28/* sloww */
29/* sloww1 */
30/* sloww2 */
31/* bsloww */
32/* bsloww1 */
33/* bsloww2 */
34/* cslow2 */
35/* csloww */
36/* csloww1 */
37/* csloww2 */
38/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */
39/* branred.c sincos32.c dosincos.c mpa.c */
40/* sincos.tbl */
41/* */
42/* An ultimate sin and routine. Given an IEEE double machine number x */
43/* it computes the correctly rounded (to nearest) value of sin(x) or cos(x) */
44/* Assumption: Machine arithmetic operations are performed in */
45/* round to nearest mode of IEEE 754 standard. */
46/* */
47/****************************************************************************/
f7eac6eb 48
f7eac6eb 49
0c59a196 50#include <errno.h>
e4d82761
UD		 51#include "endian.h"
52#include "mydefs.h"
53#include "usncs.h"
54#include "MathLib.h"
1ed0291c 55#include <math_private.h>
804360ed 56#include <fenv.h>
e4d82761 57
4aafb73c 58/* Helper macros to compute sin of the input values. */
196f7f5d 59#define POLYNOMIAL2(xx) ((((s5 * (xx) + s4) * (xx) + s3) * (xx) + s2) * (xx))
4aafb73c 60
196f7f5d 61#define POLYNOMIAL(xx) (POLYNOMIAL2 (xx) + s1)
4aafb73c
SP		 62
63/* The computed polynomial is a variation of the Taylor series expansion for
64 sin(a):
65
66 a - a^3/3! + a^5/5! - a^7/7! + a^9/9! + (1 - a^2) * da / 2
67
68 The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so
69 on. The result is returned to LHS and correction in COR. */
70#define TAYLOR_SINCOS(xx, a, da, cor) \
71({ \
72 double t = ((POLYNOMIAL (xx) * (a) - 0.5 * (da)) * (xx) + (da)); \
73 double res = (a) + t; \
74 (cor) = ((a) - res) + t; \
75 res; \
76})
77
78/* This is again a variation of the Taylor series expansion with the term
79 x^3/3! expanded into the following for better accuracy:
80
81 bb * x ^ 3 + 3 * aa * x * x1 * x2 + aa * x1 ^ 3 + aa * x2 ^ 3
82
83 The correction term is dx and bb + aa = -1/3!
84 */
85#define TAYLOR_SLOW(x0, dx, cor) \
86({ \
87 static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ \
88 double xx = (x0) * (x0); \
89 double x1 = ((x0) + th2_36) - th2_36; \
196f7f5d 90 double y = aa * x1 * x1 * x1; \
4aafb73c
SP		 91 double r = (x0) + y; \
92 double x2 = ((x0) - x1) + (dx); \
196f7f5d
SP		 93 double t = (((POLYNOMIAL2 (xx) + bb) * xx + 3.0 * aa * x1 * x2) \
94 * (x0) + aa * x2 * x2 * x2 + (dx)); \
4aafb73c
SP		 95 t = (((x0) - r) + y) + t; \
96 double res = r + t; \
97 (cor) = (r - res) + t; \
98 res; \
99})
100
b348e1e3
SP		 101#define SINCOS_TABLE_LOOKUP(u, sn, ssn, cs, ccs) \
102({ \
103 int4 k = u.i[LOW_HALF] << 2; \
104 sn = __sincostab.x[k]; \
105 ssn = __sincostab.x[k + 1]; \
106 cs = __sincostab.x[k + 2]; \
107 ccs = __sincostab.x[k + 3]; \
108})
109
31d3cc00
UD		 110#ifndef SECTION
111# define SECTION
112#endif
113
af968f62
UD		 114extern const union
115{
116 int4 i[880];
117 double x[440];
118} __sincostab attribute_hidden;
119
e4d82761 120static const double
6dbe713d
SP		 121 sn3 = -1.66666666666664880952546298448555E-01,
122 sn5 = 8.33333214285722277379541354343671E-03,
123 cs2 = 4.99999999999999999999950396842453E-01,
124 cs4 = -4.16666666666664434524222570944589E-02,
125 cs6 = 1.38888874007937613028114285595617E-03;
126
7a74607f
SP		 127static const double t22 = 0x1.8p22;
128
6dbe713d
SP		 129void __dubsin (double x, double dx, double w[]);
130void __docos (double x, double dx, double w[]);
09544cbc
SP		 131double __mpsin (double x, double dx, bool reduce_range);
132double __mpcos (double x, double dx, bool reduce_range);
6dbe713d
SP		 133static double slow (double x);
134static double slow1 (double x);
135static double slow2 (double x);
136static double sloww (double x, double dx, double orig);
137static double sloww1 (double x, double dx, double orig);
138static double sloww2 (double x, double dx, double orig, int n);
139static double bsloww (double x, double dx, double orig, int n);
140static double bsloww1 (double x, double dx, double orig, int n);
141static double bsloww2 (double x, double dx, double orig, int n);
142int __branred (double x, double *a, double *aa);
143static double cslow2 (double x);
144static double csloww (double x, double dx, double orig);
145static double csloww1 (double x, double dx, double orig);
146static double csloww2 (double x, double dx, double orig, int n);
147
6cce25f8
SP		 148/* Reduce range of X and compute sin of a + da. K is the amount by which to
149 rotate the quadrants. This allows us to use the same routine to compute cos
150 by simply rotating the quadrants by 1. */
151static inline double
152__always_inline
153reduce_and_compute (double x, double a, double da, unsigned int k)
154{
155 double retval = 0;
156 unsigned int n = __branred (x, &a, &da);
157 k = (n + k) % 4;
158 switch (k)
159 {
160 case 0:
161 if (a * a < 0.01588)
162 retval = bsloww (a, da, x, n);
163 else
164 retval = bsloww1 (a, da, x, n);
165 break;
166 case 2:
167 if (a * a < 0.01588)
168 retval = bsloww (-a, -da, x, n);
169 else
170 retval = bsloww1 (-a, -da, x, n);
171 break;
172
173 case 1:
174 case 3:
175 retval = bsloww2 (a, da, x, n);
176 break;
177 }
178 return retval;
179}
180
e4d82761
UD		 181/*******************************************************************/
182/* An ultimate sin routine. Given an IEEE double machine number x */
183/* it computes the correctly rounded (to nearest) value of sin(x) */
184/*******************************************************************/
31d3cc00
UD		 185double
186SECTION
6dbe713d
SP		 187__sin (double x)
188{
189 double xx, res, t, cor, y, s, c, sn, ssn, cs, ccs, xn, a, da, db, eps, xn1,
190 xn2;
191 mynumber u, v;
192 int4 k, m, n;
193 double retval = 0;
194
195 SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
196
197 u.x = x;
198 m = u.i[HIGH_HALF];
199 k = 0x7fffffff & m; /* no sign */
200 if (k < 0x3e500000) /* if x->0 =>sin(x)=x */
5eea0404 201 retval = x;
e4d82761 202 /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/
6dbe713d
SP		 203 else if (k < 0x3fd00000)
204 {
205 xx = x * x;
4aafb73c
SP		 206 /* Taylor series. */
207 t = POLYNOMIAL (xx) * (xx * x);
6dbe713d
SP		 208 res = x + t;
209 cor = (x - res) + t;
210 retval = (res == res + 1.07 * cor) ? res : slow (x);
6dbe713d 211 } /* else if (k < 0x3fd00000) */
e4d82761 212/*---------------------------- 0.25<|x|< 0.855469---------------------- */
6dbe713d
SP		 213 else if (k < 0x3feb6000)
214 {
196f7f5d
SP		 215 u.x = (m > 0) ? big + x : big - x;
216 y = (m > 0) ? x - (u.x - big) : x + (u.x - big);
6dbe713d
SP		 217 xx = y * y;
218 s = y + y * xx * (sn3 + xx * sn5);
219 c = xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3
SP		 220 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
221 if (m <= 0)
222 {
223 sn = -sn;
224 ssn = -ssn;
225 }
6dbe713d
SP		 226 cor = (ssn + s * ccs - sn * c) + cs * s;
227 res = sn + cor;
228 cor = (sn - res) + cor;
229 retval = (res == res + 1.096 * cor) ? res : slow1 (x);
6dbe713d 230 } /* else if (k < 0x3feb6000) */
e4d82761
UD		 231
232/*----------------------- 0.855469 <|x|<2.426265 ----------------------*/
6dbe713d
SP		 233 else if (k < 0x400368fd)
234 {
235
196f7f5d 236 y = (m > 0) ? hp0 - x : hp0 + x;
6dbe713d
SP		 237 if (y >= 0)
238 {
196f7f5d
SP		 239 u.x = big + y;
240 y = (y - (u.x - big)) + hp1;
6dbe713d
SP		 241 }
242 else
243 {
196f7f5d
SP		 244 u.x = big - y;
245 y = (-hp1) - (y + (u.x - big));
6dbe713d
SP		 246 }
247 xx = y * y;
248 s = y + y * xx * (sn3 + xx * sn5);
249 c = xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 250 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 251 cor = (ccs - s * ssn - cs * c) - sn * s;
252 res = cs + cor;
253 cor = (cs - res) + cor;
254 retval = (res == res + 1.020 * cor) ? ((m > 0) ? res : -res) : slow2 (x);
6dbe713d 255 } /* else if (k < 0x400368fd) */
e4d82761
UD		 256
257/*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
6dbe713d
SP		 258 else if (k < 0x419921FB)
259 {
196f7f5d
SP		 260 t = (x * hpinv + toint);
261 xn = t - toint;
6dbe713d 262 v.x = t;
196f7f5d 263 y = (x - xn * mp1) - xn * mp2;
6dbe713d 264 n = v.i[LOW_HALF] & 3;
196f7f5d 265 da = xn * mp3;
6dbe713d
SP		 266 a = y - da;
267 da = (y - a) - da;
268 eps = ABS (x) * 1.2e-30;
269
270 switch (n)
271 { /* quarter of unit circle */
272 case 0:
273 case 2:
274 xx = a * a;
275 if (n)
276 {
277 a = -a;
278 da = -da;
279 }
280 if (xx < 0.01588)
281 {
4aafb73c
SP		 282 /* Taylor series. */
283 res = TAYLOR_SINCOS (xx, a, da, cor);
6dbe713d
SP		 284 cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps;
285 retval = (res == res + cor) ? res : sloww (a, da, x);
f7eac6eb 286 }
6dbe713d
SP		 287 else
288 {
289 if (a > 0)
290 {
291 m = 1;
292 t = a;
293 db = da;
294 }
e4d82761 295 else
6dbe713d
SP		 296 {
297 m = 0;
298 t = -a;
299 db = -da;
300 }
196f7f5d
SP		 301 u.x = big + t;
302 y = t - (u.x - big);
6dbe713d
SP		 303 xx = y * y;
304 s = y + (db + y * xx * (sn3 + xx * sn5));
305 c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 306 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 307 cor = (ssn + s * ccs - sn * c) + cs * s;
308 res = sn + cor;
309 cor = (sn - res) + cor;
310 cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps;
311 retval = ((res == res + cor) ? ((m) ? res : -res)
312 : sloww1 (a, da, x));
e4d82761 313 }
6dbe713d
SP		 314 break;
315
316 case 1:
317 case 3:
318 if (a < 0)
319 {
320 a = -a;
321 da = -da;
322 }
196f7f5d
SP		 323 u.x = big + a;
324 y = a - (u.x - big) + da;
6dbe713d 325 xx = y * y;
b348e1e3 326 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 327 s = y + y * xx * (sn3 + xx * sn5);
328 c = xx * (cs2 + xx * (cs4 + xx * cs6));
329 cor = (ccs - s * ssn - cs * c) - sn * s;
330 res = cs + cor;
331 cor = (cs - res) + cor;
332 cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
333 retval = ((res == res + cor) ? ((n & 2) ? -res : res)
334 : sloww2 (a, da, x, n));
6dbe713d
SP		 335 break;
336 }
6dbe713d 337 } /* else if (k < 0x419921FB ) */
e4d82761
UD		 338
339/*---------------------105414350 <|x|< 281474976710656 --------------------*/
6dbe713d
SP		 340 else if (k < 0x42F00000)
341 {
196f7f5d
SP		 342 t = (x * hpinv + toint);
343 xn = t - toint;
6dbe713d
SP		 344 v.x = t;
345 xn1 = (xn + 8.0e22) - 8.0e22;
346 xn2 = xn - xn1;
196f7f5d 347 y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2);
6dbe713d 348 n = v.i[LOW_HALF] & 3;
196f7f5d 349 da = xn1 * pp3;
6dbe713d
SP		 350 t = y - da;
351 da = (y - t) - da;
196f7f5d 352 da = (da - xn2 * pp3) - xn * pp4;
6dbe713d
SP		 353 a = t + da;
354 da = (t - a) + da;
355 eps = 1.0e-24;
356
357 switch (n)
358 {
359 case 0:
360 case 2:
361 xx = a * a;
362 if (n)
363 {
364 a = -a;
365 da = -da;
366 }
367 if (xx < 0.01588)
368 {
4aafb73c
SP		 369 /* Taylor series. */
370 res = TAYLOR_SINCOS (xx, a, da, cor);
6dbe713d
SP		 371 cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps;
372 retval = (res == res + cor) ? res : bsloww (a, da, x, n);
e4d82761 373 }
6dbe713d
SP		 374 else
375 {
376 if (a > 0)
377 {
378 m = 1;
379 t = a;
380 db = da;
381 }
382 else
383 {
384 m = 0;
385 t = -a;
386 db = -da;
387 }
196f7f5d
SP		 388 u.x = big + t;
389 y = t - (u.x - big);
6dbe713d
SP		 390 xx = y * y;
391 s = y + (db + y * xx * (sn3 + xx * sn5));
392 c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 393 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 394 cor = (ssn + s * ccs - sn * c) + cs * s;
395 res = sn + cor;
396 cor = (sn - res) + cor;
397 cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps;
398 retval = ((res == res + cor) ? ((m) ? res : -res)
399 : bsloww1 (a, da, x, n));
6dbe713d
SP		 400 }
401 break;
402
403 case 1:
404 case 3:
405 if (a < 0)
406 {
407 a = -a;
408 da = -da;
409 }
196f7f5d
SP		 410 u.x = big + a;
411 y = a - (u.x - big) + da;
6dbe713d 412 xx = y * y;
b348e1e3 413 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 414 s = y + y * xx * (sn3 + xx * sn5);
415 c = xx * (cs2 + xx * (cs4 + xx * cs6));
416 cor = (ccs - s * ssn - cs * c) - sn * s;
417 res = cs + cor;
418 cor = (cs - res) + cor;
419 cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
420 retval = ((res == res + cor) ? ((n & 2) ? -res : res)
421 : bsloww2 (a, da, x, n));
6dbe713d
SP		 422 break;
423 }
424 } /* else if (k < 0x42F00000 ) */
e4d82761
UD		 425
426/* -----------------281474976710656 <|x| <2^1024----------------------------*/
6dbe713d 427 else if (k < 0x7ff00000)
6cce25f8 428 retval = reduce_and_compute (x, a, da, 0);
6dbe713d
SP		 429
430/*--------------------- |x| > 2^1024 ----------------------------------*/
431 else
432 {
433 if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
434 __set_errno (EDOM);
435 retval = x / x;
6dbe713d 436 }
804360ed 437
6dbe713d 438 return retval;
e4d82761
UD		 439}
440
441
442/*******************************************************************/
443/* An ultimate cos routine. Given an IEEE double machine number x */
444/* it computes the correctly rounded (to nearest) value of cos(x) */
445/*******************************************************************/
446
31d3cc00
UD		 447double
448SECTION
6dbe713d 449__cos (double x)
e4d82761 450{
6dbe713d
SP		 451 double y, xx, res, t, cor, s, c, sn, ssn, cs, ccs, xn, a, da, db, eps, xn1,
452 xn2;
453 mynumber u, v;
454 int4 k, m, n;
e4d82761 455
804360ed
JM		 456 double retval = 0;
457
eb92c487 458 SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
804360ed 459
e4d82761
UD		 460 u.x = x;
461 m = u.i[HIGH_HALF];
6dbe713d 462 k = 0x7fffffff & m;
e4d82761 463
5eea0404 464 /* |x|<2^-27 => cos(x)=1 */
6dbe713d 465 if (k < 0x3e400000)
5eea0404 466 retval = 1.0;
6dbe713d
SP		 467
468 else if (k < 0x3feb6000)
469 { /* 2^-27 < |x| < 0.855469 */
470 y = ABS (x);
196f7f5d
SP		 471 u.x = big + y;
472 y = y - (u.x - big);
6dbe713d
SP		 473 xx = y * y;
474 s = y + y * xx * (sn3 + xx * sn5);
475 c = xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 476 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 477 cor = (ccs - s * ssn - cs * c) - sn * s;
478 res = cs + cor;
479 cor = (cs - res) + cor;
480 retval = (res == res + 1.020 * cor) ? res : cslow2 (x);
6dbe713d
SP		 481 } /* else if (k < 0x3feb6000) */
482
483 else if (k < 0x400368fd)
484 { /* 0.855469 <|x|<2.426265 */ ;
196f7f5d
SP		 485 y = hp0 - ABS (x);
486 a = y + hp1;
487 da = (y - a) + hp1;
6dbe713d
SP		 488 xx = a * a;
489 if (xx < 0.01588)
490 {
4aafb73c 491 res = TAYLOR_SINCOS (xx, a, da, cor);
6dbe713d
SP		 492 cor = (cor > 0) ? 1.02 * cor + 1.0e-31 : 1.02 * cor - 1.0e-31;
493 retval = (res == res + cor) ? res : csloww (a, da, x);
6dbe713d
SP		 494 }
495 else
496 {
497 if (a > 0)
498 {
499 m = 1;
500 t = a;
501 db = da;
502 }
503 else
504 {
505 m = 0;
506 t = -a;
507 db = -da;
508 }
196f7f5d
SP		 509 u.x = big + t;
510 y = t - (u.x - big);
6dbe713d
SP		 511 xx = y * y;
512 s = y + (db + y * xx * (sn3 + xx * sn5));
513 c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 514 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 515 cor = (ssn + s * ccs - sn * c) + cs * s;
516 res = sn + cor;
517 cor = (sn - res) + cor;
518 cor = (cor > 0) ? 1.035 * cor + 1.0e-31 : 1.035 * cor - 1.0e-31;
519 retval = ((res == res + cor) ? ((m) ? res : -res)
520 : csloww1 (a, da, x));
6dbe713d 521 }
e4d82761 522
6dbe713d 523 } /* else if (k < 0x400368fd) */
e4d82761 524
e4d82761 525
6dbe713d
SP		 526 else if (k < 0x419921FB)
527 { /* 2.426265<|x|< 105414350 */
196f7f5d
SP		 528 t = (x * hpinv + toint);
529 xn = t - toint;
6dbe713d 530 v.x = t;
196f7f5d 531 y = (x - xn * mp1) - xn * mp2;
6dbe713d 532 n = v.i[LOW_HALF] & 3;
196f7f5d 533 da = xn * mp3;
6dbe713d
SP		 534 a = y - da;
535 da = (y - a) - da;
536 eps = ABS (x) * 1.2e-30;
537
538 switch (n)
539 {
540 case 1:
541 case 3:
542 xx = a * a;
543 if (n == 1)
544 {
545 a = -a;
546 da = -da;
547 }
548 if (xx < 0.01588)
549 {
4aafb73c 550 res = TAYLOR_SINCOS (xx, a, da, cor);
6dbe713d
SP		 551 cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps;
552 retval = (res == res + cor) ? res : csloww (a, da, x);
6dbe713d
SP		 553 }
554 else
555 {
556 if (a > 0)
557 {
558 m = 1;
559 t = a;
560 db = da;
561 }
562 else
563 {
564 m = 0;
565 t = -a;
566 db = -da;
567 }
196f7f5d
SP		 568 u.x = big + t;
569 y = t - (u.x - big);
6dbe713d
SP		 570 xx = y * y;
571 s = y + (db + y * xx * (sn3 + xx * sn5));
572 c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 573 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 574 cor = (ssn + s * ccs - sn * c) + cs * s;
575 res = sn + cor;
576 cor = (sn - res) + cor;
577 cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps;
578 retval = ((res == res + cor) ? ((m) ? res : -res)
579 : csloww1 (a, da, x));
6dbe713d
SP		 580 }
581 break;
582
583 case 0:
584 case 2:
585 if (a < 0)
586 {
587 a = -a;
588 da = -da;
589 }
196f7f5d
SP		 590 u.x = big + a;
591 y = a - (u.x - big) + da;
6dbe713d 592 xx = y * y;
b348e1e3 593 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 594 s = y + y * xx * (sn3 + xx * sn5);
595 c = xx * (cs2 + xx * (cs4 + xx * cs6));
596 cor = (ccs - s * ssn - cs * c) - sn * s;
597 res = cs + cor;
598 cor = (cs - res) + cor;
599 cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
600 retval = ((res == res + cor) ? ((n) ? -res : res)
601 : csloww2 (a, da, x, n));
6dbe713d
SP		 602 break;
603 }
604 } /* else if (k < 0x419921FB ) */
e4d82761 605
6dbe713d
SP		 606 else if (k < 0x42F00000)
607 {
196f7f5d
SP		 608 t = (x * hpinv + toint);
609 xn = t - toint;
6dbe713d
SP		 610 v.x = t;
611 xn1 = (xn + 8.0e22) - 8.0e22;
612 xn2 = xn - xn1;
196f7f5d 613 y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2);
6dbe713d 614 n = v.i[LOW_HALF] & 3;
196f7f5d 615 da = xn1 * pp3;
6dbe713d
SP		 616 t = y - da;
617 da = (y - t) - da;
196f7f5d 618 da = (da - xn2 * pp3) - xn * pp4;
6dbe713d
SP		 619 a = t + da;
620 da = (t - a) + da;
621 eps = 1.0e-24;
622
623 switch (n)
624 {
625 case 1:
626 case 3:
627 xx = a * a;
628 if (n == 1)
629 {
630 a = -a;
631 da = -da;
632 }
633 if (xx < 0.01588)
634 {
4aafb73c 635 res = TAYLOR_SINCOS (xx, a, da, cor);
6dbe713d
SP		 636 cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps;
637 retval = (res == res + cor) ? res : bsloww (a, da, x, n);
6dbe713d
SP		 638 }
639 else
640 {
641 if (a > 0)
642 {
643 m = 1;
644 t = a;
645 db = da;
646 }
647 else
648 {
649 m = 0;
650 t = -a;
651 db = -da;
652 }
196f7f5d
SP		 653 u.x = big + t;
654 y = t - (u.x - big);
6dbe713d
SP		 655 xx = y * y;
656 s = y + (db + y * xx * (sn3 + xx * sn5));
657 c = y * db + xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 658 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 659 cor = (ssn + s * ccs - sn * c) + cs * s;
660 res = sn + cor;
661 cor = (sn - res) + cor;
662 cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps;
663 retval = ((res == res + cor) ? ((m) ? res : -res)
664 : bsloww1 (a, da, x, n));
6dbe713d
SP		 665 }
666 break;
667
668 case 0:
669 case 2:
670 if (a < 0)
671 {
672 a = -a;
673 da = -da;
674 }
196f7f5d
SP		 675 u.x = big + a;
676 y = a - (u.x - big) + da;
6dbe713d 677 xx = y * y;
b348e1e3 678 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 679 s = y + y * xx * (sn3 + xx * sn5);
680 c = xx * (cs2 + xx * (cs4 + xx * cs6));
681 cor = (ccs - s * ssn - cs * c) - sn * s;
682 res = cs + cor;
683 cor = (cs - res) + cor;
684 cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
685 retval = ((res == res + cor) ? ((n) ? -res : res)
686 : bsloww2 (a, da, x, n));
6dbe713d
SP		 687 break;
688 }
689 } /* else if (k < 0x42F00000 ) */
690
6cce25f8 691 /* 281474976710656 <|x| <2^1024 */
6dbe713d 692 else if (k < 0x7ff00000)
6cce25f8 693 retval = reduce_and_compute (x, a, da, 1);
e4d82761 694
6dbe713d
SP		 695 else
696 {
697 if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
698 __set_errno (EDOM);
699 retval = x / x; /* |x| > 2^1024 */
e4d82761
UD		 700 }
701
804360ed 702 return retval;
e4d82761
UD		 703}
704
705/************************************************************************/
706/* Routine compute sin(x) for 2^-26 < |x|< 0.25 by Taylor with more */
707/* precision and if still doesn't accurate enough by mpsin or dubsin */
708/************************************************************************/
709
31d3cc00
UD		 710static double
711SECTION
6dbe713d
SP		 712slow (double x)
713{
4aafb73c
SP		 714 double res, cor, w[2];
715 res = TAYLOR_SLOW (x, 0, cor);
6dbe713d
SP		 716 if (res == res + 1.0007 * cor)
717 return res;
718 else
719 {
720 __dubsin (ABS (x), 0, w);
721 if (w[0] == w[0] + 1.000000001 * w[1])
722 return (x > 0) ? w[0] : -w[0];
723 else
09544cbc 724 return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false);
6dbe713d 725 }
e4d82761 726}
6dbe713d 727
e4d82761 728/*******************************************************************************/
c5d5d574 729/* Routine compute sin(x) for 0.25<|x|< 0.855469 by __sincostab.tbl and Taylor */
e4d82761
UD		 730/* and if result still doesn't accurate enough by mpsin or dubsin */
731/*******************************************************************************/
732
31d3cc00
UD		 733static double
734SECTION
6dbe713d
SP		 735slow1 (double x)
736{
e4d82761 737 mynumber u;
6dbe713d 738 double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, c1, c2, xx, cor, res;
6dbe713d 739 y = ABS (x);
196f7f5d
SP		 740 u.x = big + y;
741 y = y - (u.x - big);
6dbe713d
SP		 742 xx = y * y;
743 s = y * xx * (sn3 + xx * sn5);
744 c = xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 745 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d 746 y1 = (y + t22) - t22;
e4d82761 747 y2 = y - y1;
6dbe713d
SP		 748 c1 = (cs + t22) - t22;
749 c2 = (cs - c1) + ccs;
750 cor = (ssn + s * ccs + cs * s + c2 * y + c1 * y2) - sn * c;
751 y = sn + c1 * y1;
752 cor = cor + ((sn - y) + c1 * y1);
753 res = y + cor;
754 cor = (y - res) + cor;
755 if (res == res + 1.0005 * cor)
756 return (x > 0) ? res : -res;
757 else
758 {
759 __dubsin (ABS (x), 0, w);
760 if (w[0] == w[0] + 1.000000005 * w[1])
761 return (x > 0) ? w[0] : -w[0];
762 else
09544cbc 763 return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false);
6dbe713d 764 }
e4d82761 765}
6dbe713d 766
e4d82761 767/**************************************************************************/
af968f62 768/* Routine compute sin(x) for 0.855469 <|x|<2.426265 by __sincostab.tbl */
e4d82761
UD		 769/* and if result still doesn't accurate enough by mpsin or dubsin */
770/**************************************************************************/
31d3cc00
UD		 771static double
772SECTION
6dbe713d
SP		 773slow2 (double x)
774{
e4d82761 775 mynumber u;
6dbe713d 776 double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, e1, e2, xx, cor, res, del;
7a74607f 777
6dbe713d 778 y = ABS (x);
196f7f5d 779 y = hp0 - y;
6dbe713d
SP		 780 if (y >= 0)
781 {
196f7f5d
SP		 782 u.x = big + y;
783 y = y - (u.x - big);
784 del = hp1;
6dbe713d
SP		 785 }
786 else
787 {
196f7f5d
SP		 788 u.x = big - y;
789 y = -(y + (u.x - big));
790 del = -hp1;
6dbe713d
SP		 791 }
792 xx = y * y;
793 s = y * xx * (sn3 + xx * sn5);
794 c = y * del + xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 795 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 796 y1 = (y + t22) - t22;
797 y2 = (y - y1) + del;
798 e1 = (sn + t22) - t22;
799 e2 = (sn - e1) + ssn;
800 cor = (ccs - cs * c - e1 * y2 - e2 * y) - sn * s;
801 y = cs - e1 * y1;
802 cor = cor + ((cs - y) - e1 * y1);
803 res = y + cor;
804 cor = (y - res) + cor;
805 if (res == res + 1.0005 * cor)
806 return (x > 0) ? res : -res;
807 else
808 {
196f7f5d
SP		 809 y = ABS (x) - hp0;
810 y1 = y - hp1;
811 y2 = (y - y1) - hp1;
6dbe713d
SP		 812 __docos (y1, y2, w);
813 if (w[0] == w[0] + 1.000000005 * w[1])
814 return (x > 0) ? w[0] : -w[0];
815 else
09544cbc 816 return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false);
6dbe713d 817 }
e4d82761 818}
6dbe713d 819
e4d82761
UD		 820/***************************************************************************/
821/* Routine compute sin(x+dx) (Double-Length number) where x is small enough*/
822/* to use Taylor series around zero and (x+dx) */
823/* in first or third quarter of unit circle.Routine receive also */
824/* (right argument) the original value of x for computing error of */
825/* result.And if result not accurate enough routine calls mpsin1 or dubsin */
826/***************************************************************************/
827
31d3cc00
UD		 828static double
829SECTION
6dbe713d
SP		 830sloww (double x, double dx, double orig)
831{
4aafb73c 832 double y, t, res, cor, w[2], a, da, xn;
6dbe713d
SP		 833 union
834 {
835 int4 i[2];
836 double x;
837 } v;
e4d82761 838 int4 n;
4aafb73c 839 res = TAYLOR_SLOW (x, dx, cor);
6dbe713d
SP		 840 cor =
841 (cor >
842 0) ? 1.0005 * cor + ABS (orig) * 3.1e-30 : 1.0005 * cor -
843 ABS (orig) * 3.1e-30;
844 if (res == res + cor)
845 return res;
846 else
847 {
848 (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w);
849 if (w[1] > 0)
850 cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-30;
851 else
852 cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-30;
853
854 if (w[0] == w[0] + cor)
855 return (x > 0) ? w[0] : -w[0];
856 else
857 {
196f7f5d
SP		 858 t = (orig * hpinv + toint);
859 xn = t - toint;
6dbe713d 860 v.x = t;
196f7f5d 861 y = (orig - xn * mp1) - xn * mp2;
6dbe713d 862 n = v.i[LOW_HALF] & 3;
196f7f5d 863 da = xn * pp3;
6dbe713d
SP		 864 t = y - da;
865 da = (y - t) - da;
196f7f5d 866 y = xn * pp4;
6dbe713d
SP		 867 a = t - y;
868 da = ((t - a) - y) + da;
869 if (n & 2)
870 {
871 a = -a;
872 da = -da;
873 }
874 (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w);
875 if (w[1] > 0)
876 cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-40;
877 else
878 cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-40;
879
880 if (w[0] == w[0] + cor)
881 return (a > 0) ? w[0] : -w[0];
882 else
09544cbc 883 return __mpsin (orig, 0, true);
6dbe713d 884 }
e4d82761 885 }
e4d82761 886}
6dbe713d 887
e4d82761
UD		 888/***************************************************************************/
889/* Routine compute sin(x+dx) (Double-Length number) where x in first or */
890/* third quarter of unit circle.Routine receive also (right argument) the */
891/* original value of x for computing error of result.And if result not */
892/* accurate enough routine calls mpsin1 or dubsin */
893/***************************************************************************/
894
31d3cc00
UD		 895static double
896SECTION
6dbe713d
SP		 897sloww1 (double x, double dx, double orig)
898{
e4d82761 899 mynumber u;
6dbe713d 900 double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, c1, c2, xx, cor, res;
6dbe713d
SP		 901
902 y = ABS (x);
196f7f5d
SP		 903 u.x = big + y;
904 y = y - (u.x - big);
6dbe713d
SP		 905 dx = (x > 0) ? dx : -dx;
906 xx = y * y;
907 s = y * xx * (sn3 + xx * sn5);
908 c = xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 909 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 910 y1 = (y + t22) - t22;
911 y2 = (y - y1) + dx;
912 c1 = (cs + t22) - t22;
913 c2 = (cs - c1) + ccs;
914 cor = (ssn + s * ccs + cs * s + c2 * y + c1 * y2 - sn * y * dx) - sn * c;
915 y = sn + c1 * y1;
916 cor = cor + ((sn - y) + c1 * y1);
917 res = y + cor;
918 cor = (y - res) + cor;
919
920 if (cor > 0)
921 cor = 1.0005 * cor + 3.1e-30 * ABS (orig);
922 else
923 cor = 1.0005 * cor - 3.1e-30 * ABS (orig);
924
925 if (res == res + cor)
926 return (x > 0) ? res : -res;
927 else
928 {
929 __dubsin (ABS (x), dx, w);
930
931 if (w[1] > 0)
932 cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig);
933 else
934 cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig);
935
936 if (w[0] == w[0] + cor)
937 return (x > 0) ? w[0] : -w[0];
938 else
09544cbc 939 return __mpsin (orig, 0, true);
6dbe713d 940 }
e4d82761 941}
6dbe713d 942
e4d82761
UD		 943/***************************************************************************/
944/* Routine compute sin(x+dx) (Double-Length number) where x in second or */
945/* fourth quarter of unit circle.Routine receive also the original value */
946/* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
947/* accurate enough routine calls mpsin1 or dubsin */
948/***************************************************************************/
949
31d3cc00
UD		 950static double
951SECTION
6dbe713d
SP		 952sloww2 (double x, double dx, double orig, int n)
953{
e4d82761 954 mynumber u;
6dbe713d 955 double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, e1, e2, xx, cor, res;
6dbe713d
SP		 956
957 y = ABS (x);
196f7f5d
SP		 958 u.x = big + y;
959 y = y - (u.x - big);
6dbe713d
SP		 960 dx = (x > 0) ? dx : -dx;
961 xx = y * y;
962 s = y * xx * (sn3 + xx * sn5);
963 c = y * dx + xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 964 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 965
966 y1 = (y + t22) - t22;
967 y2 = (y - y1) + dx;
968 e1 = (sn + t22) - t22;
969 e2 = (sn - e1) + ssn;
970 cor = (ccs - cs * c - e1 * y2 - e2 * y) - sn * s;
971 y = cs - e1 * y1;
972 cor = cor + ((cs - y) - e1 * y1);
973 res = y + cor;
974 cor = (y - res) + cor;
975
976 if (cor > 0)
977 cor = 1.0005 * cor + 3.1e-30 * ABS (orig);
978 else
979 cor = 1.0005 * cor - 3.1e-30 * ABS (orig);
980
981 if (res == res + cor)
982 return (n & 2) ? -res : res;
983 else
984 {
985 __docos (ABS (x), dx, w);
986
987 if (w[1] > 0)
988 cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig);
989 else
990 cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig);
991
992 if (w[0] == w[0] + cor)
993 return (n & 2) ? -w[0] : w[0];
994 else
09544cbc 995 return __mpsin (orig, 0, true);
6dbe713d 996 }
e4d82761 997}
6dbe713d 998
e4d82761
UD		 999/***************************************************************************/
1000/* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */
1001/* is small enough to use Taylor series around zero and (x+dx) */
1002/* in first or third quarter of unit circle.Routine receive also */
1003/* (right argument) the original value of x for computing error of */
1004/* result.And if result not accurate enough routine calls other routines */
1005/***************************************************************************/
1006
31d3cc00
UD		 1007static double
1008SECTION
6dbe713d
SP		 1009bsloww (double x, double dx, double orig, int n)
1010{
4aafb73c
SP		 1011 double res, cor, w[2];
1012
1013 res = TAYLOR_SLOW (x, dx, cor);
6dbe713d
SP		 1014 cor = (cor > 0) ? 1.0005 * cor + 1.1e-24 : 1.0005 * cor - 1.1e-24;
1015 if (res == res + cor)
1016 return res;
1017 else
1018 {
1019 (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w);
1020 if (w[1] > 0)
1021 cor = 1.000000001 * w[1] + 1.1e-24;
1022 else
1023 cor = 1.000000001 * w[1] - 1.1e-24;
1024 if (w[0] == w[0] + cor)
1025 return (x > 0) ? w[0] : -w[0];
1026 else
09544cbc 1027 return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
6dbe713d 1028 }
e4d82761
UD		 1029}
1030
1031/***************************************************************************/
1032/* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */
1033/* in first or third quarter of unit circle.Routine receive also */
1034/* (right argument) the original value of x for computing error of result.*/
1035/* And if result not accurate enough routine calls other routines */
1036/***************************************************************************/
1037
31d3cc00
UD		 1038static double
1039SECTION
6dbe713d
SP		 1040bsloww1 (double x, double dx, double orig, int n)
1041{
1042 mynumber u;
1043 double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, c1, c2, xx, cor, res;
6dbe713d
SP		 1044
1045 y = ABS (x);
196f7f5d
SP		 1046 u.x = big + y;
1047 y = y - (u.x - big);
6dbe713d
SP		 1048 dx = (x > 0) ? dx : -dx;
1049 xx = y * y;
1050 s = y * xx * (sn3 + xx * sn5);
1051 c = xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 1052 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 1053 y1 = (y + t22) - t22;
1054 y2 = (y - y1) + dx;
1055 c1 = (cs + t22) - t22;
1056 c2 = (cs - c1) + ccs;
1057 cor = (ssn + s * ccs + cs * s + c2 * y + c1 * y2 - sn * y * dx) - sn * c;
1058 y = sn + c1 * y1;
1059 cor = cor + ((sn - y) + c1 * y1);
1060 res = y + cor;
1061 cor = (y - res) + cor;
1062 cor = (cor > 0) ? 1.0005 * cor + 1.1e-24 : 1.0005 * cor - 1.1e-24;
1063 if (res == res + cor)
1064 return (x > 0) ? res : -res;
1065 else
1066 {
1067 __dubsin (ABS (x), dx, w);
1068
1069 if (w[1] > 0)
1070 cor = 1.000000005 * w[1] + 1.1e-24;
1071 else
1072 cor = 1.000000005 * w[1] - 1.1e-24;
1073
1074 if (w[0] == w[0] + cor)
1075 return (x > 0) ? w[0] : -w[0];
1076 else
09544cbc 1077 return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
6dbe713d 1078 }
e4d82761
UD		 1079}
1080
1081/***************************************************************************/
1082/* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */
1083/* in second or fourth quarter of unit circle.Routine receive also the */
1084/* original value and quarter(n= 1or 3)of x for computing error of result. */
1085/* And if result not accurate enough routine calls other routines */
1086/***************************************************************************/
1087
31d3cc00
UD		 1088static double
1089SECTION
6dbe713d
SP		 1090bsloww2 (double x, double dx, double orig, int n)
1091{
1092 mynumber u;
1093 double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, e1, e2, xx, cor, res;
6dbe713d
SP		 1094
1095 y = ABS (x);
196f7f5d
SP		 1096 u.x = big + y;
1097 y = y - (u.x - big);
6dbe713d
SP		 1098 dx = (x > 0) ? dx : -dx;
1099 xx = y * y;
1100 s = y * xx * (sn3 + xx * sn5);
1101 c = y * dx + xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 1102 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 1103
1104 y1 = (y + t22) - t22;
1105 y2 = (y - y1) + dx;
1106 e1 = (sn + t22) - t22;
1107 e2 = (sn - e1) + ssn;
1108 cor = (ccs - cs * c - e1 * y2 - e2 * y) - sn * s;
1109 y = cs - e1 * y1;
1110 cor = cor + ((cs - y) - e1 * y1);
1111 res = y + cor;
1112 cor = (y - res) + cor;
1113 cor = (cor > 0) ? 1.0005 * cor + 1.1e-24 : 1.0005 * cor - 1.1e-24;
1114 if (res == res + cor)
1115 return (n & 2) ? -res : res;
1116 else
1117 {
1118 __docos (ABS (x), dx, w);
1119
1120 if (w[1] > 0)
1121 cor = 1.000000005 * w[1] + 1.1e-24;
1122 else
1123 cor = 1.000000005 * w[1] - 1.1e-24;
1124
1125 if (w[0] == w[0] + cor)
1126 return (n & 2) ? -w[0] : w[0];
1127 else
09544cbc 1128 return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true);
6dbe713d 1129 }
e4d82761
UD		 1130}
1131
1132/************************************************************************/
1133/* Routine compute cos(x) for 2^-27 < |x|< 0.25 by Taylor with more */
1134/* precision and if still doesn't accurate enough by mpcos or docos */
1135/************************************************************************/
1136
31d3cc00
UD		 1137static double
1138SECTION
6dbe713d
SP		 1139cslow2 (double x)
1140{
e4d82761 1141 mynumber u;
6dbe713d 1142 double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, e1, e2, xx, cor, res;
6dbe713d
SP		 1143
1144 y = ABS (x);
196f7f5d
SP		 1145 u.x = big + y;
1146 y = y - (u.x - big);
6dbe713d
SP		 1147 xx = y * y;
1148 s = y * xx * (sn3 + xx * sn5);
1149 c = xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 1150 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d 1151 y1 = (y + t22) - t22;
e4d82761 1152 y2 = y - y1;
6dbe713d
SP		 1153 e1 = (sn + t22) - t22;
1154 e2 = (sn - e1) + ssn;
1155 cor = (ccs - cs * c - e1 * y2 - e2 * y) - sn * s;
1156 y = cs - e1 * y1;
1157 cor = cor + ((cs - y) - e1 * y1);
1158 res = y + cor;
1159 cor = (y - res) + cor;
1160 if (res == res + 1.0005 * cor)
e4d82761 1161 return res;
6dbe713d
SP		 1162 else
1163 {
1164 y = ABS (x);
1165 __docos (y, 0, w);
1166 if (w[0] == w[0] + 1.000000005 * w[1])
1167 return w[0];
1168 else
09544cbc 1169 return __mpcos (x, 0, false);
6dbe713d 1170 }
e4d82761
UD		 1171}
1172
1173/***************************************************************************/
1174/* Routine compute cos(x+dx) (Double-Length number) where x is small enough*/
1175/* to use Taylor series around zero and (x+dx) .Routine receive also */
1176/* (right argument) the original value of x for computing error of */
1177/* result.And if result not accurate enough routine calls other routines */
1178/***************************************************************************/
1179
1180
31d3cc00
UD		 1181static double
1182SECTION
6dbe713d
SP		 1183csloww (double x, double dx, double orig)
1184{
4aafb73c 1185 double y, t, res, cor, w[2], a, da, xn;
6dbe713d
SP		 1186 union
1187 {
1188 int4 i[2];
1189 double x;
1190 } v;
e4d82761 1191 int4 n;
6dbe713d 1192
6dbe713d 1193 /* Taylor series */
4aafb73c 1194 t = TAYLOR_SLOW (x, dx, cor);
6dbe713d
SP		 1195
1196 if (cor > 0)
1197 cor = 1.0005 * cor + ABS (orig) * 3.1e-30;
1198 else
1199 cor = 1.0005 * cor - ABS (orig) * 3.1e-30;
1200
1201 if (res == res + cor)
1202 return res;
1203 else
1204 {
1205 (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w);
1206
1207 if (w[1] > 0)
1208 cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-30;
1209 else
1210 cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-30;
1211
1212 if (w[0] == w[0] + cor)
1213 return (x > 0) ? w[0] : -w[0];
1214 else
1215 {
196f7f5d
SP		 1216 t = (orig * hpinv + toint);
1217 xn = t - toint;
6dbe713d 1218 v.x = t;
196f7f5d 1219 y = (orig - xn * mp1) - xn * mp2;
6dbe713d 1220 n = v.i[LOW_HALF] & 3;
196f7f5d 1221 da = xn * pp3;
6dbe713d
SP		 1222 t = y - da;
1223 da = (y - t) - da;
196f7f5d 1224 y = xn * pp4;
6dbe713d
SP		 1225 a = t - y;
1226 da = ((t - a) - y) + da;
1227 if (n == 1)
1228 {
1229 a = -a;
1230 da = -da;
1231 }
1232 (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w);
1233
1234 if (w[1] > 0)
1235 cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-40;
1236 else
1237 cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-40;
1238
1239 if (w[0] == w[0] + cor)
1240 return (a > 0) ? w[0] : -w[0];
1241 else
09544cbc 1242 return __mpcos (orig, 0, true);
6dbe713d 1243 }
e4d82761 1244 }
e4d82761
UD		 1245}
1246
1247/***************************************************************************/
1248/* Routine compute sin(x+dx) (Double-Length number) where x in first or */
1249/* third quarter of unit circle.Routine receive also (right argument) the */
1250/* original value of x for computing error of result.And if result not */
1251/* accurate enough routine calls other routines */
1252/***************************************************************************/
1253
31d3cc00
UD		 1254static double
1255SECTION
6dbe713d
SP		 1256csloww1 (double x, double dx, double orig)
1257{
e4d82761 1258 mynumber u;
6dbe713d 1259 double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, c1, c2, xx, cor, res;
6dbe713d
SP		 1260
1261 y = ABS (x);
196f7f5d
SP		 1262 u.x = big + y;
1263 y = y - (u.x - big);
6dbe713d
SP		 1264 dx = (x > 0) ? dx : -dx;
1265 xx = y * y;
1266 s = y * xx * (sn3 + xx * sn5);
1267 c = xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 1268 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 1269 y1 = (y + t22) - t22;
1270 y2 = (y - y1) + dx;
1271 c1 = (cs + t22) - t22;
1272 c2 = (cs - c1) + ccs;
1273 cor = (ssn + s * ccs + cs * s + c2 * y + c1 * y2 - sn * y * dx) - sn * c;
1274 y = sn + c1 * y1;
1275 cor = cor + ((sn - y) + c1 * y1);
1276 res = y + cor;
1277 cor = (y - res) + cor;
1278
1279 if (cor > 0)
1280 cor = 1.0005 * cor + 3.1e-30 * ABS (orig);
1281 else
1282 cor = 1.0005 * cor - 3.1e-30 * ABS (orig);
1283
1284 if (res == res + cor)
1285 return (x > 0) ? res : -res;
1286 else
1287 {
1288 __dubsin (ABS (x), dx, w);
1289 if (w[1] > 0)
1290 cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig);
1291 else
1292 cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig);
1293 if (w[0] == w[0] + cor)
1294 return (x > 0) ? w[0] : -w[0];
1295 else
09544cbc 1296 return __mpcos (orig, 0, true);
6dbe713d 1297 }
e4d82761
UD		 1298}
1299
1300
1301/***************************************************************************/
1302/* Routine compute sin(x+dx) (Double-Length number) where x in second or */
1303/* fourth quarter of unit circle.Routine receive also the original value */
1304/* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
1305/* accurate enough routine calls other routines */
1306/***************************************************************************/
1307
31d3cc00
UD		 1308static double
1309SECTION
6dbe713d
SP		 1310csloww2 (double x, double dx, double orig, int n)
1311{
e4d82761 1312 mynumber u;
6dbe713d 1313 double sn, ssn, cs, ccs, s, c, w[2], y, y1, y2, e1, e2, xx, cor, res;
6dbe713d
SP		 1314
1315 y = ABS (x);
196f7f5d
SP		 1316 u.x = big + y;
1317 y = y - (u.x - big);
6dbe713d
SP		 1318 dx = (x > 0) ? dx : -dx;
1319 xx = y * y;
1320 s = y * xx * (sn3 + xx * sn5);
1321 c = y * dx + xx * (cs2 + xx * (cs4 + xx * cs6));
b348e1e3 1322 SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
6dbe713d
SP		 1323
1324 y1 = (y + t22) - t22;
1325 y2 = (y - y1) + dx;
1326 e1 = (sn + t22) - t22;
1327 e2 = (sn - e1) + ssn;
1328 cor = (ccs - cs * c - e1 * y2 - e2 * y) - sn * s;
1329 y = cs - e1 * y1;
1330 cor = cor + ((cs - y) - e1 * y1);
1331 res = y + cor;
1332 cor = (y - res) + cor;
1333
1334 if (cor > 0)
1335 cor = 1.0005 * cor + 3.1e-30 * ABS (orig);
1336 else
1337 cor = 1.0005 * cor - 3.1e-30 * ABS (orig);
1338
1339 if (res == res + cor)
1340 return (n) ? -res : res;
1341 else
1342 {
1343 __docos (ABS (x), dx, w);
1344 if (w[1] > 0)
1345 cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig);
1346 else
1347 cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig);
1348 if (w[0] == w[0] + cor)
1349 return (n) ? -w[0] : w[0];
1350 else
09544cbc 1351 return __mpcos (orig, 0, true);
6dbe713d 1352 }
e4d82761
UD		 1353}
1354
af968f62 1355#ifndef __cos
ca58f1db 1356weak_alias (__cos, cos)
af968f62
UD		 1357# ifdef NO_LONG_DOUBLE
1358strong_alias (__cos, __cosl)
1359weak_alias (__cos, cosl)
1360# endif
1361#endif
1362#ifndef __sin
ca58f1db 1363weak_alias (__sin, sin)
af968f62 1364# ifdef NO_LONG_DOUBLE
ca58f1db
UD		 1365strong_alias (__sin, __sinl)
1366weak_alias (__sin, sinl)
af968f62 1367# endif
cccda09f 1368#endif
A mirror of the official glibc repository