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