sysdeps/ieee754/dbl-64/s_sin.c

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