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