sysdeps/ieee754/dbl-64/s_sin.c

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