Back to index

glibc  2.9
s_sin.c
Go to the documentation of this file.
00001 /*
00002  * IBM Accurate Mathematical Library
00003  * written by International Business Machines Corp.
00004  * Copyright (C) 2001 Free Software Foundation
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU Lesser General Public License as published by
00008  * the Free Software Foundation; either version 2.1 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU  Lesser General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00019  */
00020 /****************************************************************************/
00021 /*                                                                          */
00022 /* MODULE_NAME:usncs.c                                                      */
00023 /*                                                                          */
00024 /* FUNCTIONS: usin                                                          */
00025 /*            ucos                                                          */
00026 /*            slow                                                          */
00027 /*            slow1                                                         */
00028 /*            slow2                                                         */
00029 /*            sloww                                                         */
00030 /*            sloww1                                                        */
00031 /*            sloww2                                                        */
00032 /*            bsloww                                                        */
00033 /*            bsloww1                                                       */
00034 /*            bsloww2                                                       */
00035 /*            cslow2                                                        */
00036 /*            csloww                                                        */
00037 /*            csloww1                                                       */
00038 /*            csloww2                                                       */
00039 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h  usncs.h                     */
00040 /*               branred.c sincos32.c dosincos.c mpa.c                      */
00041 /*               sincos.tbl                                                 */
00042 /*                                                                          */
00043 /* An ultimate sin and  routine. Given an IEEE double machine number x       */
00044 /* it computes the correctly rounded (to nearest) value of sin(x) or cos(x) */
00045 /* Assumption: Machine arithmetic operations are performed in               */
00046 /* round to nearest mode of IEEE 754 standard.                              */
00047 /*                                                                          */
00048 /****************************************************************************/
00049 
00050 
00051 #include "endian.h"
00052 #include "mydefs.h"
00053 #include "usncs.h"
00054 #include "MathLib.h"
00055 #include "sincos.tbl"
00056 #include "math_private.h"
00057 
00058 static const double
00059           sn3 = -1.66666666666664880952546298448555E-01,
00060           sn5 =  8.33333214285722277379541354343671E-03,
00061           cs2 =  4.99999999999999999999950396842453E-01,
00062           cs4 = -4.16666666666664434524222570944589E-02,
00063           cs6 =  1.38888874007937613028114285595617E-03;
00064 
00065 void __dubsin(double x, double dx, double w[]);
00066 void __docos(double x, double dx, double w[]);
00067 double __mpsin(double x, double dx);
00068 double __mpcos(double x, double dx);
00069 double __mpsin1(double x);
00070 double __mpcos1(double x);
00071 static double slow(double x);
00072 static double slow1(double x);
00073 static double slow2(double x);
00074 static double sloww(double x, double dx, double orig);
00075 static double sloww1(double x, double dx, double orig);
00076 static double sloww2(double x, double dx, double orig, int n);
00077 static double bsloww(double x, double dx, double orig, int n);
00078 static double bsloww1(double x, double dx, double orig, int n);
00079 static double bsloww2(double x, double dx, double orig, int n);
00080 int __branred(double x, double *a, double *aa);
00081 static double cslow2(double x);
00082 static double csloww(double x, double dx, double orig);
00083 static double csloww1(double x, double dx, double orig);
00084 static double csloww2(double x, double dx, double orig, int n);
00085 /*******************************************************************/
00086 /* An ultimate sin routine. Given an IEEE double machine number x   */
00087 /* it computes the correctly rounded (to nearest) value of sin(x)  */
00088 /*******************************************************************/
00089 double __sin(double x){
00090        double xx,res,t,cor,y,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2;
00091 #if 0
00092        double w[2];
00093 #endif
00094        mynumber u,v;
00095        int4 k,m,n;
00096 #if 0
00097        int4 nn;
00098 #endif
00099 
00100        u.x = x;
00101        m = u.i[HIGH_HALF];
00102        k = 0x7fffffff&m;              /* no sign           */
00103        if (k < 0x3e500000)            /* if x->0 =>sin(x)=x */
00104         return x;
00105  /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/
00106        else  if (k < 0x3fd00000){
00107          xx = x*x;
00108          /*Taylor series */
00109          t = ((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*(xx*x);
00110          res = x+t;
00111          cor = (x-res)+t;
00112          return (res == res + 1.07*cor)? res : slow(x);
00113        }    /*  else  if (k < 0x3fd00000)    */
00114 /*---------------------------- 0.25<|x|< 0.855469---------------------- */
00115        else if (k < 0x3feb6000)  {
00116          u.x=(m>0)?big.x+x:big.x-x;
00117          y=(m>0)?x-(u.x-big.x):x+(u.x-big.x);
00118          xx=y*y;
00119          s = y + y*xx*(sn3 +xx*sn5);
00120          c = xx*(cs2 +xx*(cs4 + xx*cs6));
00121          k=u.i[LOW_HALF]<<2;
00122          sn=(m>0)?sincos.x[k]:-sincos.x[k];
00123          ssn=(m>0)?sincos.x[k+1]:-sincos.x[k+1];
00124          cs=sincos.x[k+2];
00125          ccs=sincos.x[k+3];
00126          cor=(ssn+s*ccs-sn*c)+cs*s;
00127          res=sn+cor;
00128          cor=(sn-res)+cor;
00129          return (res==res+1.025*cor)? res : slow1(x);
00130        }    /*   else  if (k < 0x3feb6000)    */
00131 
00132 /*----------------------- 0.855469  <|x|<2.426265  ----------------------*/
00133        else if (k <  0x400368fd ) {
00134 
00135          y = (m>0)? hp0.x-x:hp0.x+x;
00136          if (y>=0) {
00137            u.x = big.x+y;
00138            y = (y-(u.x-big.x))+hp1.x;
00139          }
00140          else {
00141            u.x = big.x-y;
00142            y = (-hp1.x) - (y+(u.x-big.x));
00143          }
00144          xx=y*y;
00145          s = y + y*xx*(sn3 +xx*sn5);
00146          c = xx*(cs2 +xx*(cs4 + xx*cs6));
00147          k=u.i[LOW_HALF]<<2;
00148          sn=sincos.x[k];
00149          ssn=sincos.x[k+1];
00150          cs=sincos.x[k+2];
00151          ccs=sincos.x[k+3];
00152          cor=(ccs-s*ssn-cs*c)-sn*s;
00153          res=cs+cor;
00154          cor=(cs-res)+cor;
00155          return (res==res+1.020*cor)? ((m>0)?res:-res) : slow2(x);
00156        } /*   else  if (k < 0x400368fd)    */
00157 
00158 /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
00159        else if (k < 0x419921FB ) {
00160          t = (x*hpinv.x + toint.x);
00161          xn = t - toint.x;
00162          v.x = t;
00163          y = (x - xn*mp1.x) - xn*mp2.x;
00164          n =v.i[LOW_HALF]&3;
00165          da = xn*mp3.x;
00166          a=y-da;
00167          da = (y-a)-da;
00168          eps = ABS(x)*1.2e-30;
00169 
00170          switch (n) { /* quarter of unit circle */
00171          case 0:
00172          case 2:
00173            xx = a*a;
00174            if (n) {a=-a;da=-da;}
00175            if (xx < 0.01588) {
00176                       /*Taylor series */
00177              t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da;
00178              res = a+t;
00179              cor = (a-res)+t;
00180              cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps;
00181              return (res == res + cor)? res : sloww(a,da,x);
00182            }
00183            else  {
00184              if (a>0)
00185               {m=1;t=a;db=da;}
00186              else
00187               {m=0;t=-a;db=-da;}
00188              u.x=big.x+t;
00189              y=t-(u.x-big.x);
00190              xx=y*y;
00191              s = y + (db+y*xx*(sn3 +xx*sn5));
00192              c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6));
00193              k=u.i[LOW_HALF]<<2;
00194              sn=sincos.x[k];
00195              ssn=sincos.x[k+1];
00196              cs=sincos.x[k+2];
00197              ccs=sincos.x[k+3];
00198              cor=(ssn+s*ccs-sn*c)+cs*s;
00199              res=sn+cor;
00200              cor=(sn-res)+cor;
00201              cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps;
00202              return (res==res+cor)? ((m)?res:-res) : sloww1(a,da,x);
00203            }
00204            break;
00205 
00206          case 1:
00207          case 3:
00208            if (a<0)
00209              {a=-a;da=-da;}
00210            u.x=big.x+a;
00211            y=a-(u.x-big.x)+da;
00212            xx=y*y;
00213            k=u.i[LOW_HALF]<<2;
00214            sn=sincos.x[k];
00215            ssn=sincos.x[k+1];
00216            cs=sincos.x[k+2];
00217            ccs=sincos.x[k+3];
00218            s = y + y*xx*(sn3 +xx*sn5);
00219            c = xx*(cs2 +xx*(cs4 + xx*cs6));
00220            cor=(ccs-s*ssn-cs*c)-sn*s;
00221            res=cs+cor;
00222            cor=(cs-res)+cor;
00223            cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps;
00224            return (res==res+cor)? ((n&2)?-res:res) : sloww2(a,da,x,n);
00225 
00226            break;
00227 
00228          }
00229 
00230        }    /*   else  if (k <  0x419921FB )    */
00231 
00232 /*---------------------105414350 <|x|< 281474976710656 --------------------*/
00233        else if (k < 0x42F00000 ) {
00234          t = (x*hpinv.x + toint.x);
00235          xn = t - toint.x;
00236          v.x = t;
00237          xn1 = (xn+8.0e22)-8.0e22;
00238          xn2 = xn - xn1;
00239          y = ((((x - xn1*mp1.x) - xn1*mp2.x)-xn2*mp1.x)-xn2*mp2.x);
00240          n =v.i[LOW_HALF]&3;
00241          da = xn1*pp3.x;
00242          t=y-da;
00243          da = (y-t)-da;
00244          da = (da - xn2*pp3.x) -xn*pp4.x;
00245          a = t+da;
00246          da = (t-a)+da;
00247          eps = 1.0e-24;
00248 
00249          switch (n) {
00250          case 0:
00251          case 2:
00252            xx = a*a;
00253            if (n) {a=-a;da=-da;}
00254            if (xx < 0.01588) {
00255               /* Taylor series */
00256              t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da;
00257              res = a+t;
00258              cor = (a-res)+t;
00259              cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps;
00260              return (res == res + cor)? res : bsloww(a,da,x,n);
00261            }
00262            else  {
00263              if (a>0) {m=1;t=a;db=da;}
00264              else {m=0;t=-a;db=-da;}
00265              u.x=big.x+t;
00266              y=t-(u.x-big.x);
00267              xx=y*y;
00268              s = y + (db+y*xx*(sn3 +xx*sn5));
00269              c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6));
00270              k=u.i[LOW_HALF]<<2;
00271              sn=sincos.x[k];
00272              ssn=sincos.x[k+1];
00273              cs=sincos.x[k+2];
00274              ccs=sincos.x[k+3];
00275              cor=(ssn+s*ccs-sn*c)+cs*s;
00276              res=sn+cor;
00277              cor=(sn-res)+cor;
00278              cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps;
00279              return (res==res+cor)? ((m)?res:-res) : bsloww1(a,da,x,n);
00280                  }
00281            break;
00282 
00283          case 1:
00284          case 3:
00285            if (a<0)
00286              {a=-a;da=-da;}
00287            u.x=big.x+a;
00288            y=a-(u.x-big.x)+da;
00289            xx=y*y;
00290            k=u.i[LOW_HALF]<<2;
00291            sn=sincos.x[k];
00292            ssn=sincos.x[k+1];
00293            cs=sincos.x[k+2];
00294            ccs=sincos.x[k+3];
00295            s = y + y*xx*(sn3 +xx*sn5);
00296            c = xx*(cs2 +xx*(cs4 + xx*cs6));
00297            cor=(ccs-s*ssn-cs*c)-sn*s;
00298            res=cs+cor;
00299            cor=(cs-res)+cor;
00300            cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps;
00301            return (res==res+cor)? ((n&2)?-res:res) : bsloww2(a,da,x,n);
00302 
00303            break;
00304 
00305          }
00306 
00307        }    /*   else  if (k <  0x42F00000 )   */
00308 
00309 /* -----------------281474976710656 <|x| <2^1024----------------------------*/
00310        else if (k < 0x7ff00000) {
00311 
00312          n = __branred(x,&a,&da);
00313          switch (n) {
00314          case 0:
00315            if (a*a < 0.01588) return bsloww(a,da,x,n);
00316            else return bsloww1(a,da,x,n);
00317            break;
00318          case 2:
00319            if (a*a < 0.01588) return bsloww(-a,-da,x,n);
00320            else return bsloww1(-a,-da,x,n);
00321            break;
00322 
00323          case 1:
00324          case 3:
00325            return  bsloww2(a,da,x,n);
00326            break;
00327          }
00328 
00329        }    /*   else  if (k <  0x7ff00000 )    */
00330 
00331 /*--------------------- |x| > 2^1024 ----------------------------------*/
00332        else return x / x;
00333        return 0;         /* unreachable */
00334 }
00335 
00336 
00337 /*******************************************************************/
00338 /* An ultimate cos routine. Given an IEEE double machine number x   */
00339 /* it computes the correctly rounded (to nearest) value of cos(x)  */
00340 /*******************************************************************/
00341 
00342 double __cos(double x)
00343 {
00344   double y,xx,res,t,cor,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2;
00345   mynumber u,v;
00346   int4 k,m,n;
00347 
00348   u.x = x;
00349   m = u.i[HIGH_HALF];
00350   k = 0x7fffffff&m;
00351 
00352   if (k < 0x3e400000 ) return 1.0; /* |x|<2^-27 => cos(x)=1 */
00353 
00354   else if (k < 0x3feb6000 ) {/* 2^-27 < |x| < 0.855469 */
00355     y=ABS(x);
00356     u.x = big.x+y;
00357     y = y-(u.x-big.x);
00358     xx=y*y;
00359     s = y + y*xx*(sn3 +xx*sn5);
00360     c = xx*(cs2 +xx*(cs4 + xx*cs6));
00361     k=u.i[LOW_HALF]<<2;
00362     sn=sincos.x[k];
00363     ssn=sincos.x[k+1];
00364     cs=sincos.x[k+2];
00365     ccs=sincos.x[k+3];
00366     cor=(ccs-s*ssn-cs*c)-sn*s;
00367     res=cs+cor;
00368     cor=(cs-res)+cor;
00369     return (res==res+1.020*cor)? res : cslow2(x);
00370 
00371 }    /*   else  if (k < 0x3feb6000)    */
00372 
00373   else if (k <  0x400368fd ) {/* 0.855469  <|x|<2.426265  */;
00374     y=hp0.x-ABS(x);
00375     a=y+hp1.x;
00376     da=(y-a)+hp1.x;
00377     xx=a*a;
00378     if (xx < 0.01588) {
00379       t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da;
00380       res = a+t;
00381       cor = (a-res)+t;
00382       cor = (cor>0)? 1.02*cor+1.0e-31 : 1.02*cor -1.0e-31;
00383       return (res == res + cor)? res : csloww(a,da,x);
00384     }
00385     else  {
00386       if (a>0) {m=1;t=a;db=da;}
00387       else {m=0;t=-a;db=-da;}
00388       u.x=big.x+t;
00389       y=t-(u.x-big.x);
00390       xx=y*y;
00391       s = y + (db+y*xx*(sn3 +xx*sn5));
00392       c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6));
00393       k=u.i[LOW_HALF]<<2;
00394       sn=sincos.x[k];
00395       ssn=sincos.x[k+1];
00396       cs=sincos.x[k+2];
00397       ccs=sincos.x[k+3];
00398       cor=(ssn+s*ccs-sn*c)+cs*s;
00399       res=sn+cor;
00400       cor=(sn-res)+cor;
00401       cor = (cor>0)? 1.035*cor+1.0e-31 : 1.035*cor-1.0e-31;
00402       return (res==res+cor)? ((m)?res:-res) : csloww1(a,da,x);
00403 }
00404 
00405 }    /*   else  if (k < 0x400368fd)    */
00406 
00407 
00408   else if (k < 0x419921FB ) {/* 2.426265<|x|< 105414350 */
00409     t = (x*hpinv.x + toint.x);
00410     xn = t - toint.x;
00411     v.x = t;
00412     y = (x - xn*mp1.x) - xn*mp2.x;
00413     n =v.i[LOW_HALF]&3;
00414     da = xn*mp3.x;
00415     a=y-da;
00416     da = (y-a)-da;
00417     eps = ABS(x)*1.2e-30;
00418 
00419     switch (n) {
00420     case 1:
00421     case 3:
00422       xx = a*a;
00423       if (n == 1) {a=-a;da=-da;}
00424       if (xx < 0.01588) {
00425        t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da;
00426        res = a+t;
00427        cor = (a-res)+t;
00428        cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps;
00429        return (res == res + cor)? res : csloww(a,da,x);
00430       }
00431       else  {
00432        if (a>0) {m=1;t=a;db=da;}
00433        else {m=0;t=-a;db=-da;}
00434        u.x=big.x+t;
00435        y=t-(u.x-big.x);
00436        xx=y*y;
00437        s = y + (db+y*xx*(sn3 +xx*sn5));
00438        c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6));
00439        k=u.i[LOW_HALF]<<2;
00440        sn=sincos.x[k];
00441        ssn=sincos.x[k+1];
00442        cs=sincos.x[k+2];
00443        ccs=sincos.x[k+3];
00444        cor=(ssn+s*ccs-sn*c)+cs*s;
00445        res=sn+cor;
00446        cor=(sn-res)+cor;
00447        cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps;
00448        return (res==res+cor)? ((m)?res:-res) : csloww1(a,da,x);
00449       }
00450       break;
00451 
00452   case 0:
00453     case 2:
00454       if (a<0) {a=-a;da=-da;}
00455       u.x=big.x+a;
00456       y=a-(u.x-big.x)+da;
00457       xx=y*y;
00458       k=u.i[LOW_HALF]<<2;
00459       sn=sincos.x[k];
00460       ssn=sincos.x[k+1];
00461       cs=sincos.x[k+2];
00462       ccs=sincos.x[k+3];
00463       s = y + y*xx*(sn3 +xx*sn5);
00464       c = xx*(cs2 +xx*(cs4 + xx*cs6));
00465       cor=(ccs-s*ssn-cs*c)-sn*s;
00466       res=cs+cor;
00467       cor=(cs-res)+cor;
00468       cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps;
00469       return (res==res+cor)? ((n)?-res:res) : csloww2(a,da,x,n);
00470 
00471            break;
00472 
00473     }
00474 
00475   }    /*   else  if (k <  0x419921FB )    */
00476 
00477 
00478   else if (k < 0x42F00000 ) {
00479     t = (x*hpinv.x + toint.x);
00480     xn = t - toint.x;
00481     v.x = t;
00482     xn1 = (xn+8.0e22)-8.0e22;
00483     xn2 = xn - xn1;
00484     y = ((((x - xn1*mp1.x) - xn1*mp2.x)-xn2*mp1.x)-xn2*mp2.x);
00485     n =v.i[LOW_HALF]&3;
00486     da = xn1*pp3.x;
00487     t=y-da;
00488     da = (y-t)-da;
00489     da = (da - xn2*pp3.x) -xn*pp4.x;
00490     a = t+da;
00491     da = (t-a)+da;
00492     eps = 1.0e-24;
00493 
00494     switch (n) {
00495     case 1:
00496     case 3:
00497       xx = a*a;
00498       if (n==1) {a=-a;da=-da;}
00499       if (xx < 0.01588) {
00500        t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da;
00501        res = a+t;
00502        cor = (a-res)+t;
00503        cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps;
00504        return (res == res + cor)? res : bsloww(a,da,x,n);
00505       }
00506       else  {
00507        if (a>0) {m=1;t=a;db=da;}
00508        else {m=0;t=-a;db=-da;}
00509        u.x=big.x+t;
00510        y=t-(u.x-big.x);
00511        xx=y*y;
00512        s = y + (db+y*xx*(sn3 +xx*sn5));
00513        c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6));
00514        k=u.i[LOW_HALF]<<2;
00515        sn=sincos.x[k];
00516        ssn=sincos.x[k+1];
00517        cs=sincos.x[k+2];
00518        ccs=sincos.x[k+3];
00519        cor=(ssn+s*ccs-sn*c)+cs*s;
00520        res=sn+cor;
00521        cor=(sn-res)+cor;
00522        cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps;
00523        return (res==res+cor)? ((m)?res:-res) : bsloww1(a,da,x,n);
00524       }
00525       break;
00526 
00527     case 0:
00528     case 2:
00529       if (a<0) {a=-a;da=-da;}
00530       u.x=big.x+a;
00531       y=a-(u.x-big.x)+da;
00532       xx=y*y;
00533       k=u.i[LOW_HALF]<<2;
00534       sn=sincos.x[k];
00535       ssn=sincos.x[k+1];
00536       cs=sincos.x[k+2];
00537       ccs=sincos.x[k+3];
00538       s = y + y*xx*(sn3 +xx*sn5);
00539       c = xx*(cs2 +xx*(cs4 + xx*cs6));
00540       cor=(ccs-s*ssn-cs*c)-sn*s;
00541       res=cs+cor;
00542       cor=(cs-res)+cor;
00543       cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps;
00544       return (res==res+cor)? ((n)?-res:res) : bsloww2(a,da,x,n);
00545       break;
00546 
00547     }
00548 
00549   }    /*   else  if (k <  0x42F00000 )    */
00550 
00551   else if (k < 0x7ff00000) {/* 281474976710656 <|x| <2^1024 */
00552 
00553     n = __branred(x,&a,&da);
00554     switch (n) {
00555     case 1:
00556       if (a*a < 0.01588) return bsloww(-a,-da,x,n);
00557       else return bsloww1(-a,-da,x,n);
00558       break;
00559               case 3:
00560                 if (a*a < 0.01588) return bsloww(a,da,x,n);
00561                 else return bsloww1(a,da,x,n);
00562                 break;
00563 
00564     case 0:
00565     case 2:
00566       return  bsloww2(a,da,x,n);
00567       break;
00568     }
00569 
00570   }    /*   else  if (k <  0x7ff00000 )    */
00571 
00572 
00573 
00574 
00575   else return x / x; /* |x| > 2^1024 */
00576   return 0;
00577 
00578 }
00579 
00580 /************************************************************************/
00581 /*  Routine compute sin(x) for  2^-26 < |x|< 0.25 by  Taylor with more   */
00582 /* precision  and if still doesn't accurate enough by mpsin   or dubsin */
00583 /************************************************************************/
00584 
00585 static double slow(double x) {
00586 static const double th2_36 = 206158430208.0;   /*    1.5*2**37   */
00587  double y,x1,x2,xx,r,t,res,cor,w[2];
00588  x1=(x+th2_36)-th2_36;
00589  y = aa.x*x1*x1*x1;
00590  r=x+y;
00591  x2=x-x1;
00592  xx=x*x;
00593  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;
00594  t=((x-r)+y)+t;
00595  res=r+t;
00596  cor = (r-res)+t;
00597  if (res == res + 1.0007*cor) return res;
00598  else {
00599    __dubsin(ABS(x),0,w);
00600    if (w[0] == w[0]+1.000000001*w[1]) return (x>0)?w[0]:-w[0];
00601    else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
00602  }
00603 }
00604 /*******************************************************************************/
00605 /* Routine compute sin(x) for   0.25<|x|< 0.855469 by  sincos.tbl   and Taylor */
00606 /* and if result still doesn't accurate enough by mpsin   or dubsin            */
00607 /*******************************************************************************/
00608 
00609 static double slow1(double x) {
00610   mynumber u;
00611   double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res;
00612   static const double t22 = 6291456.0;
00613   int4 k;
00614   y=ABS(x);
00615   u.x=big.x+y;
00616   y=y-(u.x-big.x);
00617   xx=y*y;
00618   s = y*xx*(sn3 +xx*sn5);
00619   c = xx*(cs2 +xx*(cs4 + xx*cs6));
00620   k=u.i[LOW_HALF]<<2;
00621   sn=sincos.x[k];          /* Data          */
00622   ssn=sincos.x[k+1];       /*  from         */
00623   cs=sincos.x[k+2];        /*   tables      */
00624   ccs=sincos.x[k+3];       /*    sincos.tbl */
00625   y1 = (y+t22)-t22;
00626   y2 = y - y1;
00627   c1 = (cs+t22)-t22;
00628   c2=(cs-c1)+ccs;
00629   cor=(ssn+s*ccs+cs*s+c2*y+c1*y2)-sn*c;
00630   y=sn+c1*y1;
00631   cor = cor+((sn-y)+c1*y1);
00632   res=y+cor;
00633   cor=(y-res)+cor;
00634   if (res == res+1.0005*cor) return (x>0)?res:-res;
00635   else {
00636     __dubsin(ABS(x),0,w);
00637     if (w[0] == w[0]+1.000000005*w[1]) return (x>0)?w[0]:-w[0];
00638     else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
00639   }
00640 }
00641 /**************************************************************************/
00642 /*  Routine compute sin(x) for   0.855469  <|x|<2.426265  by  sincos.tbl  */
00643 /* and if result still doesn't accurate enough by mpsin   or dubsin       */
00644 /**************************************************************************/
00645 static double slow2(double x) {
00646   mynumber u;
00647   double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res,del;
00648   static const double t22 = 6291456.0;
00649   int4 k;
00650   y=ABS(x);
00651   y = hp0.x-y;
00652   if (y>=0) {
00653     u.x = big.x+y;
00654     y = y-(u.x-big.x);
00655     del = hp1.x;
00656   }
00657   else {
00658     u.x = big.x-y;
00659     y = -(y+(u.x-big.x));
00660     del = -hp1.x;
00661   }
00662   xx=y*y;
00663   s = y*xx*(sn3 +xx*sn5);
00664   c = y*del+xx*(cs2 +xx*(cs4 + xx*cs6));
00665   k=u.i[LOW_HALF]<<2;
00666   sn=sincos.x[k];
00667   ssn=sincos.x[k+1];
00668   cs=sincos.x[k+2];
00669   ccs=sincos.x[k+3];
00670   y1 = (y+t22)-t22;
00671   y2 = (y - y1)+del;
00672   e1 = (sn+t22)-t22;
00673   e2=(sn-e1)+ssn;
00674   cor=(ccs-cs*c-e1*y2-e2*y)-sn*s;
00675   y=cs-e1*y1;
00676   cor = cor+((cs-y)-e1*y1);
00677   res=y+cor;
00678   cor=(y-res)+cor;
00679   if (res == res+1.0005*cor) return (x>0)?res:-res;
00680   else {
00681     y=ABS(x)-hp0.x;
00682     y1=y-hp1.x;
00683     y2=(y-y1)-hp1.x;
00684     __docos(y1,y2,w);
00685     if (w[0] == w[0]+1.000000005*w[1]) return (x>0)?w[0]:-w[0];
00686     else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
00687   }
00688 }
00689 /***************************************************************************/
00690 /*  Routine compute sin(x+dx) (Double-Length number) where x is small enough*/
00691 /* to use Taylor series around zero and   (x+dx)                            */
00692 /* in first or third quarter of unit circle.Routine receive also            */
00693 /* (right argument) the  original   value of x for computing error of      */
00694 /* result.And if result not accurate enough routine calls mpsin1 or dubsin */
00695 /***************************************************************************/
00696 
00697 static double sloww(double x,double dx, double orig) {
00698   static const double th2_36 = 206158430208.0;   /*    1.5*2**37   */
00699   double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn;
00700   union {int4 i[2]; double x;} v;
00701   int4 n;
00702   x1=(x+th2_36)-th2_36;
00703   y = aa.x*x1*x1*x1;
00704   r=x+y;
00705   x2=(x-x1)+dx;
00706   xx=x*x;
00707   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;
00708   t=((x-r)+y)+t;
00709   res=r+t;
00710   cor = (r-res)+t;
00711   cor = (cor>0)? 1.0005*cor+ABS(orig)*3.1e-30 : 1.0005*cor-ABS(orig)*3.1e-30;
00712   if (res == res + cor) return res;
00713   else {
00714     (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w);
00715     cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-30 : 1.000000001*w[1] - ABS(orig)*1.1e-30;
00716     if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
00717     else {
00718       t = (orig*hpinv.x + toint.x);
00719       xn = t - toint.x;
00720       v.x = t;
00721       y = (orig - xn*mp1.x) - xn*mp2.x;
00722       n =v.i[LOW_HALF]&3;
00723       da = xn*pp3.x;
00724       t=y-da;
00725       da = (y-t)-da;
00726       y = xn*pp4.x;
00727       a = t - y;
00728       da = ((t-a)-y)+da;
00729       if (n&2) {a=-a; da=-da;}
00730       (a>0)? __dubsin(a,da,w) : __dubsin(-a,-da,w);
00731       cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-40 : 1.000000001*w[1] - ABS(orig)*1.1e-40;
00732       if (w[0] == w[0]+cor) return (a>0)?w[0]:-w[0];
00733       else return __mpsin1(orig);
00734     }
00735   }
00736 }
00737 /***************************************************************************/
00738 /*  Routine compute sin(x+dx)   (Double-Length number) where x in first or  */
00739 /*  third quarter of unit circle.Routine receive also (right argument) the  */
00740 /*  original   value of x for computing error of result.And if result not  */
00741 /* accurate enough routine calls  mpsin1   or dubsin                       */
00742 /***************************************************************************/
00743 
00744 static double sloww1(double x, double dx, double orig) {
00745   mynumber u;
00746   double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res;
00747   static const double t22 = 6291456.0;
00748   int4 k;
00749   y=ABS(x);
00750   u.x=big.x+y;
00751   y=y-(u.x-big.x);
00752   dx=(x>0)?dx:-dx;
00753   xx=y*y;
00754   s = y*xx*(sn3 +xx*sn5);
00755   c = xx*(cs2 +xx*(cs4 + xx*cs6));
00756   k=u.i[LOW_HALF]<<2;
00757   sn=sincos.x[k];
00758   ssn=sincos.x[k+1];
00759   cs=sincos.x[k+2];
00760   ccs=sincos.x[k+3];
00761   y1 = (y+t22)-t22;
00762   y2 = (y - y1)+dx;
00763   c1 = (cs+t22)-t22;
00764   c2=(cs-c1)+ccs;
00765   cor=(ssn+s*ccs+cs*s+c2*y+c1*y2-sn*y*dx)-sn*c;
00766   y=sn+c1*y1;
00767   cor = cor+((sn-y)+c1*y1);
00768   res=y+cor;
00769   cor=(y-res)+cor;
00770   cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
00771   if (res == res + cor) return (x>0)?res:-res;
00772   else {
00773     __dubsin(ABS(x),dx,w);
00774     cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
00775     if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
00776   else  return __mpsin1(orig);
00777   }
00778 }
00779 /***************************************************************************/
00780 /*  Routine compute sin(x+dx)   (Double-Length number) where x in second or */
00781 /*  fourth quarter of unit circle.Routine receive also  the  original value */
00782 /* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
00783 /* accurate enough routine calls  mpsin1   or dubsin                       */
00784 /***************************************************************************/
00785 
00786 static double sloww2(double x, double dx, double orig, int n) {
00787   mynumber u;
00788   double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res;
00789   static const double t22 = 6291456.0;
00790   int4 k;
00791   y=ABS(x);
00792   u.x=big.x+y;
00793   y=y-(u.x-big.x);
00794   dx=(x>0)?dx:-dx;
00795   xx=y*y;
00796   s = y*xx*(sn3 +xx*sn5);
00797   c = y*dx+xx*(cs2 +xx*(cs4 + xx*cs6));
00798   k=u.i[LOW_HALF]<<2;
00799   sn=sincos.x[k];
00800   ssn=sincos.x[k+1];
00801   cs=sincos.x[k+2];
00802   ccs=sincos.x[k+3];
00803 
00804   y1 = (y+t22)-t22;
00805   y2 = (y - y1)+dx;
00806   e1 = (sn+t22)-t22;
00807   e2=(sn-e1)+ssn;
00808   cor=(ccs-cs*c-e1*y2-e2*y)-sn*s;
00809   y=cs-e1*y1;
00810   cor = cor+((cs-y)-e1*y1);
00811   res=y+cor;
00812   cor=(y-res)+cor;
00813   cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
00814   if (res == res + cor) return (n&2)?-res:res;
00815   else {
00816    __docos(ABS(x),dx,w);
00817    cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
00818    if (w[0] == w[0]+cor) return (n&2)?-w[0]:w[0];
00819    else  return __mpsin1(orig);
00820   }
00821 }
00822 /***************************************************************************/
00823 /*  Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x   */
00824 /* is small enough to use Taylor series around zero and   (x+dx)            */
00825 /* in first or third quarter of unit circle.Routine receive also            */
00826 /* (right argument) the  original   value of x for computing error of      */
00827 /* result.And if result not accurate enough routine calls other routines    */
00828 /***************************************************************************/
00829 
00830 static double bsloww(double x,double dx, double orig,int n) {
00831   static const double th2_36 = 206158430208.0;   /*    1.5*2**37   */
00832   double y,x1,x2,xx,r,t,res,cor,w[2];
00833 #if 0
00834   double a,da,xn;
00835   union {int4 i[2]; double x;} v;
00836 #endif
00837   x1=(x+th2_36)-th2_36;
00838   y = aa.x*x1*x1*x1;
00839   r=x+y;
00840   x2=(x-x1)+dx;
00841   xx=x*x;
00842   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;
00843   t=((x-r)+y)+t;
00844   res=r+t;
00845   cor = (r-res)+t;
00846   cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24;
00847   if (res == res + cor) return res;
00848   else {
00849     (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w);
00850     cor = (w[1]>0)? 1.000000001*w[1] + 1.1e-24 : 1.000000001*w[1] - 1.1e-24;
00851     if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
00852     else return (n&1)?__mpcos1(orig):__mpsin1(orig);
00853   }
00854 }
00855 
00856 /***************************************************************************/
00857 /*  Routine compute sin(x+dx)  or cos(x+dx) (Double-Length number) where x  */
00858 /* in first or third quarter of unit circle.Routine receive also            */
00859 /* (right argument) the original  value of x for computing error of result.*/
00860 /* And if result not  accurate enough routine calls  other routines         */
00861 /***************************************************************************/
00862 
00863 static double bsloww1(double x, double dx, double orig,int n) {
00864 mynumber u;
00865  double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res;
00866  static const double t22 = 6291456.0;
00867  int4 k;
00868  y=ABS(x);
00869  u.x=big.x+y;
00870  y=y-(u.x-big.x);
00871  dx=(x>0)?dx:-dx;
00872  xx=y*y;
00873  s = y*xx*(sn3 +xx*sn5);
00874  c = xx*(cs2 +xx*(cs4 + xx*cs6));
00875  k=u.i[LOW_HALF]<<2;
00876  sn=sincos.x[k];
00877  ssn=sincos.x[k+1];
00878  cs=sincos.x[k+2];
00879  ccs=sincos.x[k+3];
00880  y1 = (y+t22)-t22;
00881  y2 = (y - y1)+dx;
00882  c1 = (cs+t22)-t22;
00883  c2=(cs-c1)+ccs;
00884  cor=(ssn+s*ccs+cs*s+c2*y+c1*y2-sn*y*dx)-sn*c;
00885  y=sn+c1*y1;
00886  cor = cor+((sn-y)+c1*y1);
00887  res=y+cor;
00888  cor=(y-res)+cor;
00889  cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24;
00890  if (res == res + cor) return (x>0)?res:-res;
00891  else {
00892    __dubsin(ABS(x),dx,w);
00893    cor = (w[1]>0)? 1.000000005*w[1]+1.1e-24: 1.000000005*w[1]-1.1e-24;
00894    if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
00895    else  return (n&1)?__mpcos1(orig):__mpsin1(orig);
00896  }
00897 }
00898 
00899 /***************************************************************************/
00900 /*  Routine compute sin(x+dx)  or cos(x+dx) (Double-Length number) where x  */
00901 /* in second or fourth quarter of unit circle.Routine receive also  the     */
00902 /* original value and quarter(n= 1or 3)of x for computing error of result.  */
00903 /* And if result not accurate enough routine calls  other routines          */
00904 /***************************************************************************/
00905 
00906 static double bsloww2(double x, double dx, double orig, int n) {
00907 mynumber u;
00908  double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res;
00909  static const double t22 = 6291456.0;
00910  int4 k;
00911  y=ABS(x);
00912  u.x=big.x+y;
00913  y=y-(u.x-big.x);
00914  dx=(x>0)?dx:-dx;
00915  xx=y*y;
00916  s = y*xx*(sn3 +xx*sn5);
00917  c = y*dx+xx*(cs2 +xx*(cs4 + xx*cs6));
00918  k=u.i[LOW_HALF]<<2;
00919  sn=sincos.x[k];
00920  ssn=sincos.x[k+1];
00921  cs=sincos.x[k+2];
00922  ccs=sincos.x[k+3];
00923 
00924  y1 = (y+t22)-t22;
00925  y2 = (y - y1)+dx;
00926  e1 = (sn+t22)-t22;
00927  e2=(sn-e1)+ssn;
00928  cor=(ccs-cs*c-e1*y2-e2*y)-sn*s;
00929  y=cs-e1*y1;
00930  cor = cor+((cs-y)-e1*y1);
00931  res=y+cor;
00932  cor=(y-res)+cor;
00933  cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24;
00934  if (res == res + cor) return (n&2)?-res:res;
00935  else {
00936    __docos(ABS(x),dx,w);
00937    cor = (w[1]>0)? 1.000000005*w[1]+1.1e-24 : 1.000000005*w[1]-1.1e-24;
00938    if (w[0] == w[0]+cor) return (n&2)?-w[0]:w[0];
00939    else  return (n&1)?__mpsin1(orig):__mpcos1(orig);
00940  }
00941 }
00942 
00943 /************************************************************************/
00944 /*  Routine compute cos(x) for  2^-27 < |x|< 0.25 by  Taylor with more   */
00945 /* precision  and if still doesn't accurate enough by mpcos   or docos  */
00946 /************************************************************************/
00947 
00948 static double cslow2(double x) {
00949   mynumber u;
00950   double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res;
00951   static const double t22 = 6291456.0;
00952   int4 k;
00953   y=ABS(x);
00954   u.x = big.x+y;
00955   y = y-(u.x-big.x);
00956   xx=y*y;
00957   s = y*xx*(sn3 +xx*sn5);
00958   c = xx*(cs2 +xx*(cs4 + xx*cs6));
00959   k=u.i[LOW_HALF]<<2;
00960   sn=sincos.x[k];
00961   ssn=sincos.x[k+1];
00962   cs=sincos.x[k+2];
00963   ccs=sincos.x[k+3];
00964   y1 = (y+t22)-t22;
00965   y2 = y - y1;
00966   e1 = (sn+t22)-t22;
00967   e2=(sn-e1)+ssn;
00968   cor=(ccs-cs*c-e1*y2-e2*y)-sn*s;
00969   y=cs-e1*y1;
00970   cor = cor+((cs-y)-e1*y1);
00971   res=y+cor;
00972   cor=(y-res)+cor;
00973   if (res == res+1.0005*cor)
00974     return res;
00975   else {
00976     y=ABS(x);
00977     __docos(y,0,w);
00978     if (w[0] == w[0]+1.000000005*w[1]) return w[0];
00979     else return __mpcos(x,0);
00980   }
00981 }
00982 
00983 /***************************************************************************/
00984 /*  Routine compute cos(x+dx) (Double-Length number) where x is small enough*/
00985 /* to use Taylor series around zero and   (x+dx) .Routine receive also      */
00986 /* (right argument) the  original   value of x for computing error of      */
00987 /* result.And if result not accurate enough routine calls other routines    */
00988 /***************************************************************************/
00989 
00990 
00991 static double csloww(double x,double dx, double orig) {
00992   static const double th2_36 = 206158430208.0;   /*    1.5*2**37   */
00993   double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn;
00994   union {int4 i[2]; double x;} v;
00995   int4 n;
00996   x1=(x+th2_36)-th2_36;
00997   y = aa.x*x1*x1*x1;
00998   r=x+y;
00999   x2=(x-x1)+dx;
01000   xx=x*x;
01001     /* Taylor series */
01002   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;
01003   t=((x-r)+y)+t;
01004   res=r+t;
01005   cor = (r-res)+t;
01006   cor = (cor>0)? 1.0005*cor+ABS(orig)*3.1e-30 : 1.0005*cor-ABS(orig)*3.1e-30;
01007   if (res == res + cor) return res;
01008   else {
01009     (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w);
01010     cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-30 : 1.000000001*w[1] - ABS(orig)*1.1e-30;
01011     if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
01012     else {
01013       t = (orig*hpinv.x + toint.x);
01014       xn = t - toint.x;
01015       v.x = t;
01016       y = (orig - xn*mp1.x) - xn*mp2.x;
01017       n =v.i[LOW_HALF]&3;
01018       da = xn*pp3.x;
01019       t=y-da;
01020       da = (y-t)-da;
01021       y = xn*pp4.x;
01022       a = t - y;
01023       da = ((t-a)-y)+da;
01024       if (n==1) {a=-a; da=-da;}
01025       (a>0)? __dubsin(a,da,w) : __dubsin(-a,-da,w);
01026       cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-40 : 1.000000001*w[1] - ABS(orig)*1.1e-40;
01027       if (w[0] == w[0]+cor) return (a>0)?w[0]:-w[0];
01028       else return __mpcos1(orig);
01029     }
01030   }
01031 }
01032 
01033 /***************************************************************************/
01034 /*  Routine compute sin(x+dx)   (Double-Length number) where x in first or  */
01035 /*  third quarter of unit circle.Routine receive also (right argument) the  */
01036 /*  original   value of x for computing error of result.And if result not  */
01037 /* accurate enough routine calls  other routines                            */
01038 /***************************************************************************/
01039 
01040 static double csloww1(double x, double dx, double orig) {
01041   mynumber u;
01042   double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res;
01043   static const double t22 = 6291456.0;
01044   int4 k;
01045   y=ABS(x);
01046   u.x=big.x+y;
01047   y=y-(u.x-big.x);
01048   dx=(x>0)?dx:-dx;
01049   xx=y*y;
01050   s = y*xx*(sn3 +xx*sn5);
01051   c = xx*(cs2 +xx*(cs4 + xx*cs6));
01052   k=u.i[LOW_HALF]<<2;
01053   sn=sincos.x[k];
01054   ssn=sincos.x[k+1];
01055   cs=sincos.x[k+2];
01056   ccs=sincos.x[k+3];
01057   y1 = (y+t22)-t22;
01058   y2 = (y - y1)+dx;
01059   c1 = (cs+t22)-t22;
01060   c2=(cs-c1)+ccs;
01061   cor=(ssn+s*ccs+cs*s+c2*y+c1*y2-sn*y*dx)-sn*c;
01062   y=sn+c1*y1;
01063   cor = cor+((sn-y)+c1*y1);
01064   res=y+cor;
01065   cor=(y-res)+cor;
01066   cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
01067   if (res == res + cor) return (x>0)?res:-res;
01068   else {
01069     __dubsin(ABS(x),dx,w);
01070     cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
01071     if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
01072     else  return __mpcos1(orig);
01073   }
01074 }
01075 
01076 
01077 /***************************************************************************/
01078 /*  Routine compute sin(x+dx)   (Double-Length number) where x in second or */
01079 /*  fourth quarter of unit circle.Routine receive also  the  original value */
01080 /* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
01081 /* accurate enough routine calls  other routines                            */
01082 /***************************************************************************/
01083 
01084 static double csloww2(double x, double dx, double orig, int n) {
01085   mynumber u;
01086   double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res;
01087   static const double t22 = 6291456.0;
01088   int4 k;
01089   y=ABS(x);
01090   u.x=big.x+y;
01091   y=y-(u.x-big.x);
01092   dx=(x>0)?dx:-dx;
01093   xx=y*y;
01094   s = y*xx*(sn3 +xx*sn5);
01095   c = y*dx+xx*(cs2 +xx*(cs4 + xx*cs6));
01096   k=u.i[LOW_HALF]<<2;
01097   sn=sincos.x[k];
01098   ssn=sincos.x[k+1];
01099   cs=sincos.x[k+2];
01100   ccs=sincos.x[k+3];
01101 
01102   y1 = (y+t22)-t22;
01103   y2 = (y - y1)+dx;
01104   e1 = (sn+t22)-t22;
01105   e2=(sn-e1)+ssn;
01106   cor=(ccs-cs*c-e1*y2-e2*y)-sn*s;
01107   y=cs-e1*y1;
01108   cor = cor+((cs-y)-e1*y1);
01109   res=y+cor;
01110   cor=(y-res)+cor;
01111   cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
01112   if (res == res + cor) return (n)?-res:res;
01113   else {
01114     __docos(ABS(x),dx,w);
01115     cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
01116     if (w[0] == w[0]+cor) return (n)?-w[0]:w[0];
01117     else  return __mpcos1(orig);
01118   }
01119 }
01120 
01121 weak_alias (__cos, cos)
01122 weak_alias (__sin, sin)
01123 
01124 #ifdef NO_LONG_DOUBLE
01125 strong_alias (__sin, __sinl)
01126 weak_alias (__sin, sinl)
01127 strong_alias (__cos, __cosl)
01128 weak_alias (__cos, cosl)
01129 #endif