libgcc-math/dbl-64/s_sin.c

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