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