Back to index

glibc  2.9
e_asin.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 /*     MODULE_NAME:uasncs.c                                       */
00022 /*                                                                */
00023 /*     FUNCTIONS: uasin                                           */
00024 /*                uacos                                           */
00025 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h  usncs.h           */
00026 /*               doasin.c sincos32.c dosincos.c mpa.c             */
00027 /*               sincos.tbl  asincos.tbl  powtwo.tbl root.tbl     */
00028 /*                                                                */
00029 /* Ultimate asin/acos routines. Given an IEEE double machine      */
00030 /* number x, compute the correctly rounded value of               */
00031 /* arcsin(x)or arccos(x)  according to the function called.       */
00032 /* Assumption: Machine arithmetic operations are performed in     */
00033 /* round to nearest mode of IEEE 754 standard.                    */
00034 /*                                                                */
00035 /******************************************************************/
00036 #include "endian.h"
00037 #include "mydefs.h"
00038 #include "asincos.tbl"
00039 #include "root.tbl"
00040 #include "powtwo.tbl"
00041 #include "MathLib.h"
00042 #include "uasncs.h"
00043 #include "math_private.h"
00044 
00045 void __doasin(double x, double dx, double w[]);
00046 void __dubsin(double x, double dx, double v[]);
00047 void __dubcos(double x, double dx, double v[]);
00048 void __docos(double x, double dx, double v[]);
00049 double __sin32(double x, double res, double res1);
00050 double __cos32(double x, double res, double res1);
00051 
00052 /***************************************************************************/
00053 /* An ultimate asin routine. Given an IEEE double machine number x         */
00054 /* it computes the correctly rounded (to nearest) value of arcsin(x)       */
00055 /***************************************************************************/
00056 double __ieee754_asin(double x){
00057   double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2];
00058   mynumber u,v;
00059   int4 k,m,n;
00060 #if 0
00061   int4 nn;
00062 #endif
00063 
00064   u.x = x;
00065   m = u.i[HIGH_HALF];
00066   k = 0x7fffffff&m;              /* no sign */
00067 
00068   if (k < 0x3e500000) return x;  /* for x->0 => sin(x)=x */
00069   /*----------------------2^-26 <= |x| < 2^ -3    -----------------*/
00070   else
00071   if (k < 0x3fc00000) {
00072     x2 = x*x;
00073     t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
00074     res = x+t;         /*  res=arcsin(x) according to Taylor series  */
00075     cor = (x-res)+t;
00076     if (res == res+1.025*cor) return res;
00077     else {
00078       x1 = x+big;
00079       xx = x*x;
00080       x1 -= big;
00081       x2 = x - x1;
00082       p = x1*x1*x1;
00083       s1 = a1.x*p;
00084       s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
00085             ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
00086       res1 = x+s1;
00087       s2 = ((x-res1)+s1)+s2;
00088       res = res1+s2;
00089       cor = (res1-res)+s2;
00090       if (res == res+1.00014*cor) return res;
00091       else {
00092        __doasin(x,0,w);
00093        if (w[0]==(w[0]+1.00000001*w[1])) return w[0];
00094        else {
00095          y=ABS(x);
00096          res=ABS(w[0]);
00097          res1=ABS(w[0]+1.1*w[1]);
00098          return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
00099        }
00100       }
00101     }
00102   }
00103   /*---------------------0.125 <= |x| < 0.5 -----------------------------*/
00104   else if (k < 0x3fe00000) {
00105     if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
00106     else n = 11*((k&0x000fffff)>>14)+352;
00107     if (m>0) xx = x - asncs.x[n];
00108     else xx = -x - asncs.x[n];
00109     t = asncs.x[n+1]*xx;
00110     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
00111      +xx*asncs.x[n+6]))))+asncs.x[n+7];
00112     t+=p;
00113     res =asncs.x[n+8] +t;
00114     cor = (asncs.x[n+8]-res)+t;
00115     if (res == res+1.05*cor) return (m>0)?res:-res;
00116     else {
00117       r=asncs.x[n+8]+xx*asncs.x[n+9];
00118       t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
00119       res = r+t;
00120       cor = (r-res)+t;
00121       if (res == res+1.0005*cor) return (m>0)?res:-res;
00122       else {
00123        res1=res+1.1*cor;
00124        z=0.5*(res1-res);
00125        __dubsin(res,z,w);
00126        z=(w[0]-ABS(x))+w[1];
00127        if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
00128        else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
00129        else {
00130          y=ABS(x);
00131          return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
00132        }
00133       }
00134     }
00135   }    /*   else  if (k < 0x3fe00000)    */
00136   /*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/
00137   else
00138   if (k < 0x3fe80000) {
00139     n = 1056+((k&0x000fe000)>>11)*3;
00140     if (m>0) xx = x - asncs.x[n];
00141     else xx = -x - asncs.x[n];
00142     t = asncs.x[n+1]*xx;
00143     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
00144           +xx*(asncs.x[n+6]+xx*asncs.x[n+7])))))+asncs.x[n+8];
00145     t+=p;
00146     res =asncs.x[n+9] +t;
00147     cor = (asncs.x[n+9]-res)+t;
00148     if (res == res+1.01*cor) return (m>0)?res:-res;
00149     else {
00150       r=asncs.x[n+9]+xx*asncs.x[n+10];
00151       t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
00152       res = r+t;
00153       cor = (r-res)+t;
00154       if (res == res+1.0005*cor) return (m>0)?res:-res;
00155       else {
00156        res1=res+1.1*cor;
00157        z=0.5*(res1-res);
00158        __dubsin(res,z,w);
00159        z=(w[0]-ABS(x))+w[1];
00160        if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
00161        else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
00162        else {
00163          y=ABS(x);
00164          return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
00165        }
00166       }
00167     }
00168   }    /*   else  if (k < 0x3fe80000)    */
00169   /*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/
00170   else
00171   if (k < 0x3fed8000) {
00172     n = 992+((k&0x000fe000)>>13)*13;
00173     if (m>0) xx = x - asncs.x[n];
00174     else xx = -x - asncs.x[n];
00175     t = asncs.x[n+1]*xx;
00176     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
00177      +xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+xx*asncs.x[n+8]))))))+asncs.x[n+9];
00178     t+=p;
00179     res =asncs.x[n+10] +t;
00180     cor = (asncs.x[n+10]-res)+t;
00181     if (res == res+1.01*cor) return (m>0)?res:-res;
00182     else {
00183       r=asncs.x[n+10]+xx*asncs.x[n+11];
00184       t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
00185       res = r+t;
00186       cor = (r-res)+t;
00187       if (res == res+1.0008*cor) return (m>0)?res:-res;
00188       else {
00189        res1=res+1.1*cor;
00190        z=0.5*(res1-res);
00191        y=hp0.x-res;
00192        z=((hp0.x-y)-res)+(hp1.x-z);
00193        __dubcos(y,z,w);
00194        z=(w[0]-ABS(x))+w[1];
00195        if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
00196        else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
00197        else {
00198          y=ABS(x);
00199          return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
00200        }
00201       }
00202     }
00203   }    /*   else  if (k < 0x3fed8000)    */
00204   /*-------------------0.921875 <= |x| < 0.953125 ------------------------*/
00205   else
00206   if (k < 0x3fee8000) {
00207     n = 884+((k&0x000fe000)>>13)*14;
00208     if (m>0) xx = x - asncs.x[n];
00209     else xx = -x - asncs.x[n];
00210     t = asncs.x[n+1]*xx;
00211     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
00212                       xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
00213                     +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
00214                       xx*asncs.x[n+9])))))))+asncs.x[n+10];
00215     t+=p;
00216     res =asncs.x[n+11] +t;
00217     cor = (asncs.x[n+11]-res)+t;
00218     if (res == res+1.01*cor) return (m>0)?res:-res;
00219     else {
00220       r=asncs.x[n+11]+xx*asncs.x[n+12];
00221       t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
00222       res = r+t;
00223       cor = (r-res)+t;
00224       if (res == res+1.0007*cor) return (m>0)?res:-res;
00225       else {
00226        res1=res+1.1*cor;
00227        z=0.5*(res1-res);
00228        y=(hp0.x-res)-z;
00229        z=y+hp1.x;
00230        y=(y-z)+hp1.x;
00231        __dubcos(z,y,w);
00232        z=(w[0]-ABS(x))+w[1];
00233        if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
00234        else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
00235        else {
00236          y=ABS(x);
00237          return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
00238        }
00239       }
00240     }
00241   }    /*   else  if (k < 0x3fee8000)    */
00242 
00243   /*--------------------0.953125 <= |x| < 0.96875 ------------------------*/
00244   else
00245   if (k < 0x3fef0000) {
00246     n = 768+((k&0x000fe000)>>13)*15;
00247     if (m>0) xx = x - asncs.x[n];
00248     else xx = -x - asncs.x[n];
00249     t = asncs.x[n+1]*xx;
00250     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
00251                          xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
00252                       +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
00253                     xx*(asncs.x[n+9]+xx*asncs.x[n+10]))))))))+asncs.x[n+11];
00254     t+=p;
00255     res =asncs.x[n+12] +t;
00256     cor = (asncs.x[n+12]-res)+t;
00257     if (res == res+1.01*cor) return (m>0)?res:-res;
00258     else {
00259       r=asncs.x[n+12]+xx*asncs.x[n+13];
00260       t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
00261       res = r+t;
00262       cor = (r-res)+t;
00263       if (res == res+1.0007*cor) return (m>0)?res:-res;
00264       else {
00265        res1=res+1.1*cor;
00266        z=0.5*(res1-res);
00267        y=(hp0.x-res)-z;
00268        z=y+hp1.x;
00269        y=(y-z)+hp1.x;
00270        __dubcos(z,y,w);
00271        z=(w[0]-ABS(x))+w[1];
00272        if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
00273        else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
00274        else {
00275          y=ABS(x);
00276          return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
00277        }
00278       }
00279     }
00280   }    /*   else  if (k < 0x3fef0000)    */
00281   /*--------------------0.96875 <= |x| < 1 --------------------------------*/
00282   else
00283   if (k<0x3ff00000)  {
00284     z = 0.5*((m>0)?(1.0-x):(1.0+x));
00285     v.x=z;
00286     k=v.i[HIGH_HALF];
00287     t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
00288     r=1.0-t*t*z;
00289     t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
00290     c=t*z;
00291     t=c*(1.5-0.5*t*c);
00292     y=(c+t24)-t24;
00293     cc = (z-y*y)/(t+y);
00294     p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
00295     cor = (hp1.x - 2.0*cc)-2.0*(y+cc)*p;
00296     res1 = hp0.x - 2.0*y;
00297     res =res1 + cor;
00298     if (res == res+1.003*((res1-res)+cor)) return (m>0)?res:-res;
00299     else {
00300       c=y+cc;
00301       cc=(y-c)+cc;
00302       __doasin(c,cc,w);
00303       res1=hp0.x-2.0*w[0];
00304       cor=((hp0.x-res1)-2.0*w[0])+(hp1.x-2.0*w[1]);
00305       res = res1+cor;
00306       cor = (res1-res)+cor;
00307       if (res==(res+1.0000001*cor)) return (m>0)?res:-res;
00308       else {
00309        y=ABS(x);
00310        res1=res+1.1*cor;
00311        return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
00312       }
00313     }
00314   }    /*   else  if (k < 0x3ff00000)    */
00315   /*---------------------------- |x|>=1 -------------------------------*/
00316   else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?hp0.x:-hp0.x;
00317   else
00318   if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x;
00319   else {
00320     u.i[HIGH_HALF]=0x7ff00000;
00321     v.i[HIGH_HALF]=0x7ff00000;
00322     u.i[LOW_HALF]=0;
00323     v.i[LOW_HALF]=0;
00324     return u.x/v.x;  /* NaN */
00325  }
00326 }
00327 
00328 /*******************************************************************/
00329 /*                                                                 */
00330 /*         End of arcsine,  below is arccosine                     */
00331 /*                                                                 */
00332 /*******************************************************************/
00333 
00334 double __ieee754_acos(double x)
00335 {
00336   double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2],eps;
00337 #if 0
00338   double fc;
00339 #endif
00340   mynumber u,v;
00341   int4 k,m,n;
00342 #if 0
00343   int4 nn;
00344 #endif
00345   u.x = x;
00346   m = u.i[HIGH_HALF];
00347   k = 0x7fffffff&m;
00348   /*-------------------  |x|<2.77556*10^-17 ----------------------*/
00349   if (k < 0x3c880000) return hp0.x;
00350 
00351   /*-----------------  2.77556*10^-17 <= |x| < 2^-3 --------------*/
00352   else
00353   if (k < 0x3fc00000) {
00354     x2 = x*x;
00355     t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
00356     r=hp0.x-x;
00357     cor=(((hp0.x-r)-x)+hp1.x)-t;
00358     res = r+cor;
00359     cor = (r-res)+cor;
00360     if (res == res+1.004*cor) return res;
00361     else {
00362       x1 = x+big;
00363       xx = x*x;
00364       x1 -= big;
00365       x2 = x - x1;
00366       p = x1*x1*x1;
00367       s1 = a1.x*p;
00368       s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
00369            ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
00370       res1 = x+s1;
00371       s2 = ((x-res1)+s1)+s2;
00372       r=hp0.x-res1;
00373       cor=(((hp0.x-r)-res1)+hp1.x)-s2;
00374       res = r+cor;
00375       cor = (r-res)+cor;
00376       if (res == res+1.00004*cor) return res;
00377       else {
00378        __doasin(x,0,w);
00379        r=hp0.x-w[0];
00380        cor=((hp0.x-r)-w[0])+(hp1.x-w[1]);
00381        res=r+cor;
00382        cor=(r-res)+cor;
00383        if (res ==(res +1.00000001*cor)) return res;
00384        else {
00385          res1=res+1.1*cor;
00386          return __cos32(x,res,res1);
00387        }
00388       }
00389     }
00390   }    /*   else  if (k < 0x3fc00000)    */
00391   /*----------------------  0.125 <= |x| < 0.5 --------------------*/
00392   else
00393   if (k < 0x3fe00000) {
00394     if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
00395     else n = 11*((k&0x000fffff)>>14)+352;
00396     if (m>0) xx = x - asncs.x[n];
00397     else xx = -x - asncs.x[n];
00398     t = asncs.x[n+1]*xx;
00399     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
00400                    xx*(asncs.x[n+5]+xx*asncs.x[n+6]))))+asncs.x[n+7];
00401     t+=p;
00402     y = (m>0)?(hp0.x-asncs.x[n+8]):(hp0.x+asncs.x[n+8]);
00403     t = (m>0)?(hp1.x-t):(hp1.x+t);
00404     res = y+t;
00405     if (res == res+1.02*((y-res)+t)) return res;
00406     else {
00407       r=asncs.x[n+8]+xx*asncs.x[n+9];
00408       t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
00409       if (m>0)
00410        {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; }
00411       else
00412        {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); }
00413       res = p+t;
00414       cor = (p-res)+t;
00415       if (res == (res+1.0002*cor)) return res;
00416       else {
00417        res1=res+1.1*cor;
00418        z=0.5*(res1-res);
00419        __docos(res,z,w);
00420        z=(w[0]-x)+w[1];
00421        if (z>1.0e-27) return max(res,res1);
00422        else if (z<-1.0e-27) return min(res,res1);
00423        else return __cos32(x,res,res1);
00424       }
00425     }
00426   }    /*   else  if (k < 0x3fe00000)    */
00427 
00428   /*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/
00429   else
00430   if (k < 0x3fe80000) {
00431     n = 1056+((k&0x000fe000)>>11)*3;
00432     if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
00433     else {xx = -x - asncs.x[n]; eps=1.02; }
00434     t = asncs.x[n+1]*xx;
00435     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
00436                    xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+
00437                    xx*asncs.x[n+7])))))+asncs.x[n+8];
00438     t+=p;
00439    y = (m>0)?(hp0.x-asncs.x[n+9]):(hp0.x+asncs.x[n+9]);
00440    t = (m>0)?(hp1.x-t):(hp1.x+t);
00441    res = y+t;
00442    if (res == res+eps*((y-res)+t)) return res;
00443    else {
00444      r=asncs.x[n+9]+xx*asncs.x[n+10];
00445      t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
00446      if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0004; }
00447      else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0002; }
00448      res = p+t;
00449      cor = (p-res)+t;
00450      if (res == (res+eps*cor)) return res;
00451      else {
00452        res1=res+1.1*cor;
00453        z=0.5*(res1-res);
00454        __docos(res,z,w);
00455        z=(w[0]-x)+w[1];
00456        if (z>1.0e-27) return max(res,res1);
00457        else if (z<-1.0e-27) return min(res,res1);
00458        else return __cos32(x,res,res1);
00459      }
00460    }
00461   }    /*   else  if (k < 0x3fe80000)    */
00462 
00463 /*------------------------- 0.75 <= |x| < 0.921875 -------------*/
00464   else
00465   if (k < 0x3fed8000) {
00466     n = 992+((k&0x000fe000)>>13)*13;
00467     if (m>0) {xx = x - asncs.x[n]; eps = 1.04; }
00468     else {xx = -x - asncs.x[n]; eps = 1.01; }
00469     t = asncs.x[n+1]*xx;
00470     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
00471                       xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+
00472                       xx*asncs.x[n+8]))))))+asncs.x[n+9];
00473     t+=p;
00474     y = (m>0)?(hp0.x-asncs.x[n+10]):(hp0.x+asncs.x[n+10]);
00475     t = (m>0)?(hp1.x-t):(hp1.x+t);
00476     res = y+t;
00477     if (res == res+eps*((y-res)+t)) return res;
00478     else {
00479       r=asncs.x[n+10]+xx*asncs.x[n+11];
00480       t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
00481       if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0032; }
00482       else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0008; }
00483       res = p+t;
00484       cor = (p-res)+t;
00485       if (res == (res+eps*cor)) return res;
00486       else {
00487        res1=res+1.1*cor;
00488        z=0.5*(res1-res);
00489        __docos(res,z,w);
00490        z=(w[0]-x)+w[1];
00491        if (z>1.0e-27) return max(res,res1);
00492        else if (z<-1.0e-27) return min(res,res1);
00493        else return __cos32(x,res,res1);
00494       }
00495     }
00496   }    /*   else  if (k < 0x3fed8000)    */
00497 
00498 /*-------------------0.921875 <= |x| < 0.953125 ------------------*/
00499   else
00500   if (k < 0x3fee8000) {
00501     n = 884+((k&0x000fe000)>>13)*14;
00502     if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
00503     else {xx = -x - asncs.x[n]; eps =1.005; }
00504     t = asncs.x[n+1]*xx;
00505     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
00506                    xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
00507                  +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
00508                    xx*asncs.x[n+9])))))))+asncs.x[n+10];
00509     t+=p;
00510     y = (m>0)?(hp0.x-asncs.x[n+11]):(hp0.x+asncs.x[n+11]);
00511     t = (m>0)?(hp1.x-t):(hp1.x+t);
00512     res = y+t;
00513     if (res == res+eps*((y-res)+t)) return res;
00514     else {
00515       r=asncs.x[n+11]+xx*asncs.x[n+12];
00516       t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
00517       if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
00518       else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
00519       res = p+t;
00520       cor = (p-res)+t;
00521       if (res == (res+eps*cor)) return res;
00522       else {
00523        res1=res+1.1*cor;
00524        z=0.5*(res1-res);
00525        __docos(res,z,w);
00526        z=(w[0]-x)+w[1];
00527        if (z>1.0e-27) return max(res,res1);
00528        else if (z<-1.0e-27) return min(res,res1);
00529        else return __cos32(x,res,res1);
00530       }
00531     }
00532   }    /*   else  if (k < 0x3fee8000)    */
00533 
00534   /*--------------------0.953125 <= |x| < 0.96875 ----------------*/
00535   else
00536   if (k < 0x3fef0000) {
00537     n = 768+((k&0x000fe000)>>13)*15;
00538     if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
00539     else {xx = -x - asncs.x[n]; eps=1.005;}
00540     t = asncs.x[n+1]*xx;
00541     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
00542             xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
00543            +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+xx*(asncs.x[n+9]+
00544             xx*asncs.x[n+10]))))))))+asncs.x[n+11];
00545     t+=p;
00546     y = (m>0)?(hp0.x-asncs.x[n+12]):(hp0.x+asncs.x[n+12]);
00547    t = (m>0)?(hp1.x-t):(hp1.x+t);
00548    res = y+t;
00549    if (res == res+eps*((y-res)+t)) return res;
00550    else {
00551      r=asncs.x[n+12]+xx*asncs.x[n+13];
00552      t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
00553      if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
00554      else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
00555      res = p+t;
00556      cor = (p-res)+t;
00557      if (res == (res+eps*cor)) return res;
00558      else {
00559        res1=res+1.1*cor;
00560        z=0.5*(res1-res);
00561        __docos(res,z,w);
00562        z=(w[0]-x)+w[1];
00563        if (z>1.0e-27) return max(res,res1);
00564        else if (z<-1.0e-27) return min(res,res1);
00565        else return __cos32(x,res,res1);
00566      }
00567    }
00568   }    /*   else  if (k < 0x3fef0000)    */
00569   /*-----------------0.96875 <= |x| < 1 ---------------------------*/
00570 
00571   else
00572   if (k<0x3ff00000)  {
00573     z = 0.5*((m>0)?(1.0-x):(1.0+x));
00574     v.x=z;
00575     k=v.i[HIGH_HALF];
00576     t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
00577     r=1.0-t*t*z;
00578     t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
00579     c=t*z;
00580     t=c*(1.5-0.5*t*c);
00581     y = (t27*c+c)-t27*c;
00582     cc = (z-y*y)/(t+y);
00583     p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
00584     if (m<0) {
00585       cor = (hp1.x - cc)-(y+cc)*p;
00586       res1 = hp0.x - y;
00587       res =res1 + cor;
00588       if (res == res+1.002*((res1-res)+cor)) return (res+res);
00589       else {
00590        c=y+cc;
00591        cc=(y-c)+cc;
00592        __doasin(c,cc,w);
00593        res1=hp0.x-w[0];
00594        cor=((hp0.x-res1)-w[0])+(hp1.x-w[1]);
00595        res = res1+cor;
00596        cor = (res1-res)+cor;
00597        if (res==(res+1.000001*cor)) return (res+res);
00598        else {
00599          res=res+res;
00600          res1=res+1.2*cor;
00601          return __cos32(x,res,res1);
00602        }
00603       }
00604     }
00605     else {
00606       cor = cc+p*(y+cc);
00607       res = y + cor;
00608       if (res == res+1.03*((y-res)+cor)) return (res+res);
00609       else {
00610        c=y+cc;
00611        cc=(y-c)+cc;
00612        __doasin(c,cc,w);
00613        res = w[0];
00614        cor=w[1];
00615        if (res==(res+1.000001*cor)) return (res+res);
00616        else {
00617          res=res+res;
00618          res1=res+1.2*cor;
00619          return __cos32(x,res,res1);
00620        }
00621       }
00622     }
00623   }    /*   else  if (k < 0x3ff00000)    */
00624 
00625   /*---------------------------- |x|>=1 -----------------------*/
00626   else
00627   if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?0:2.0*hp0.x;
00628   else
00629   if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x;
00630   else {
00631     u.i[HIGH_HALF]=0x7ff00000;
00632     v.i[HIGH_HALF]=0x7ff00000;
00633     u.i[LOW_HALF]=0;
00634     v.i[LOW_HALF]=0;
00635     return u.x/v.x;
00636   }
00637 }