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