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