Use protected name.
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / dbl-64 / s_sin.c
1 /*
2  * IBM Accurate Mathematical Library
3  * Copyright (c) International Business Machines Corp., 2001
4  *
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.
9  *
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.
14  *
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.
18  */
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 /****************************************************************************/
48
49
50 #include "endian.h"
51 #include "mydefs.h"
52 #include "usncs.h"
53 #include "MathLib.h"
54 #include "sincos.tbl"
55
56 static 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
63 void __dubsin(double x, double dx, double w[]);
64 void __docos(double x, double dx, double w[]);
65 double __mpsin(double x, double dx);
66 double __mpcos(double x, double dx);
67 double __mpsin1(double x);
68 double __mpcos1(double x);
69 static double slow(double x);
70 static double slow1(double x);
71 static double slow2(double x);
72 static double sloww(double x, double dx, double orig);
73 static double sloww1(double x, double dx, double orig);
74 static double sloww2(double x, double dx, double orig, int n);
75 static double bsloww(double x, double dx, double orig, int n);
76 static double bsloww1(double x, double dx, double orig, int n);
77 static double bsloww2(double x, double dx, double orig, int n);
78 int branred(double x, double *a, double *aa);
79 static double cslow2(double x);
80 static double csloww(double x, double dx, double orig);
81 static double csloww1(double x, double dx, double orig);
82 static 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 /*******************************************************************/
87 double __sin(double x){
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
92         mynumber u,v;
93         int4 k,m,n;
94 #if 0
95         int4 nn;
96 #endif
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);
180             }
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
340 double __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);
401 }
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
583 static double slow(double x) {
584 static 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];
599    else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
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
607 static 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];
636     else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
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 /**************************************************************************/
643 static 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];
684     else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
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
695 static 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];
731       else return __mpsin1(orig);
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
742 static 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];
774   else  return __mpsin1(orig);
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
784 static 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];
817    else  return __mpsin1(orig);
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
828 static double bsloww(double x,double dx, double orig,int n) {
829   static const double th2_36 = 206158430208.0;   /*    1.5*2**37   */
830   double y,x1,x2,xx,r,t,res,cor,w[2];
831 #if 0
832   double a,da,xn;
833   union {int4 i[2]; double x;} v;
834 #endif
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];
850     else return (n&1)?__mpcos1(orig):__mpsin1(orig);
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
861 static double bsloww1(double x, double dx, double orig,int n) {
862 mynumber 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];
893    else  return (n&1)?__mpcos1(orig):__mpsin1(orig);
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
904 static double bsloww2(double x, double dx, double orig, int n) {
905 mynumber 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];
937    else  return (n&1)?__mpsin1(orig):__mpcos1(orig);
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
946 static 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];
977     else return __mpcos(x,0);
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
989 static 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];
1026       else return __mpcos1(orig);
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
1038 static 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];
1070     else  return __mpcos1(orig);
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
1082 static 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];
1115     else  return __mpcos1(orig);
1116   }
1117 }
1118
1119 weak_alias (__cos, cos)
1120 weak_alias (__sin, sin)
1121
1122 #ifdef NO_LONG_DOUBLE
1123 strong_alias (__sin, __sinl)
1124 weak_alias (__sin, sinl)
1125 strong_alias (__cos, __cosl)
1126 weak_alias (__cos, cosl)
1127 #endif