sysdeps/ieee754/dbl-64/e_asin.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001, 2011 Free Software Foundation
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation; either version 2.1 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public License
  17  * along with this program; if not, see <http://www.gnu.org/licenses/>.
  18  */
  19 /******************************************************************/
  20 /*     MODULE_NAME:uasncs.c                                       */
  21 /*                                                                */
  22 /*     FUNCTIONS: uasin                                           */
  23 /*                uacos                                           */
  24 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h  usncs.h           */
  25 /*               doasin.c sincos32.c dosincos.c mpa.c             */
  26 /*               sincos.tbl  asincos.tbl  powtwo.tbl root.tbl     */
  27 /*                                                                */
  28 /* Ultimate asin/acos routines. Given an IEEE double machine      */
  29 /* number x, compute the correctly rounded value of               */
  30 /* arcsin(x)or arccos(x)  according to the function called.       */
  31 /* Assumption: Machine arithmetic operations are performed in     */
  32 /* round to nearest mode of IEEE 754 standard.                    */
  33 /*                                                                */
  34 /******************************************************************/
  35 #include "endian.h"
  36 #include "mydefs.h"
  37 #include "asincos.tbl"
  38 #include "root.tbl"
  39 #include "powtwo.tbl"
  40 #include "MathLib.h"
  41 #include "uasncs.h"
  42 #include "math_private.h"
  43
  44 #ifndef SECTION
  45 # define SECTION
  46 #endif
  47
  48 void __doasin(double x, double dx, double w[]);
  49 void __dubsin(double x, double dx, double v[]);
  50 void __dubcos(double x, double dx, double v[]);
  51 void __docos(double x, double dx, double v[]);
  52 double __sin32(double x, double res, double res1);
  53 double __cos32(double x, double res, double res1);
  54
  55 /***************************************************************************/
  56 /* An ultimate asin routine. Given an IEEE double machine number x         */
  57 /* it computes the correctly rounded (to nearest) value of arcsin(x)       */
  58 /***************************************************************************/
  59 double
  60 SECTION
  61 __ieee754_asin(double x){
  62   double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2];
  63   mynumber u,v;
  64   int4 k,m,n;
  65 #if 0
  66   int4 nn;
  67 #endif
  68
  69   u.x = x;
  70   m = u.i[HIGH_HALF];
  71   k = 0x7fffffff&m;              /* no sign */
  72
  73   if (k < 0x3e500000) return x;  /* for x->0 => sin(x)=x */
  74   /*----------------------2^-26 <= |x| < 2^ -3    -----------------*/
  75   else
  76   if (k < 0x3fc00000) {
  77     x2 = x*x;
  78     t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
  79     res = x+t;         /*  res=arcsin(x) according to Taylor series  */
  80     cor = (x-res)+t;
  81     if (res == res+1.025*cor) return res;
  82     else {
  83       x1 = x+big;
  84       xx = x*x;
  85       x1 -= big;
  86       x2 = x - x1;
  87       p = x1*x1*x1;
  88       s1 = a1.x*p;
  89       s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
  90              ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
  91       res1 = x+s1;
  92       s2 = ((x-res1)+s1)+s2;
  93       res = res1+s2;
  94       cor = (res1-res)+s2;
  95       if (res == res+1.00014*cor) return res;
  96       else {
  97         __doasin(x,0,w);
  98         if (w[0]==(w[0]+1.00000001*w[1])) return w[0];
  99         else {
 100           y=ABS(x);
 101           res=ABS(w[0]);
 102           res1=ABS(w[0]+1.1*w[1]);
 103           return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
 104         }
 105       }
 106     }
 107   }
 108   /*---------------------0.125 <= |x| < 0.5 -----------------------------*/
 109   else if (k < 0x3fe00000) {
 110     if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
 111     else n = 11*((k&0x000fffff)>>14)+352;
 112     if (m>0) xx = x - asncs.x[n];
 113     else xx = -x - asncs.x[n];
 114     t = asncs.x[n+1]*xx;
 115     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
 116      +xx*asncs.x[n+6]))))+asncs.x[n+7];
 117     t+=p;
 118     res =asncs.x[n+8] +t;
 119     cor = (asncs.x[n+8]-res)+t;
 120     if (res == res+1.05*cor) return (m>0)?res:-res;
 121     else {
 122       r=asncs.x[n+8]+xx*asncs.x[n+9];
 123       t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
 124       res = r+t;
 125       cor = (r-res)+t;
 126       if (res == res+1.0005*cor) return (m>0)?res:-res;
 127       else {
 128         res1=res+1.1*cor;
 129         z=0.5*(res1-res);
 130         __dubsin(res,z,w);
 131         z=(w[0]-ABS(x))+w[1];
 132         if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
 133         else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
 134         else {
 135           y=ABS(x);
 136           return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
 137         }
 138       }
 139     }
 140   }    /*   else  if (k < 0x3fe00000)    */
 141   /*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/
 142   else
 143   if (k < 0x3fe80000) {
 144     n = 1056+((k&0x000fe000)>>11)*3;
 145     if (m>0) xx = x - asncs.x[n];
 146     else xx = -x - asncs.x[n];
 147     t = asncs.x[n+1]*xx;
 148     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
 149            +xx*(asncs.x[n+6]+xx*asncs.x[n+7])))))+asncs.x[n+8];
 150     t+=p;
 151     res =asncs.x[n+9] +t;
 152     cor = (asncs.x[n+9]-res)+t;
 153     if (res == res+1.01*cor) return (m>0)?res:-res;
 154     else {
 155       r=asncs.x[n+9]+xx*asncs.x[n+10];
 156       t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
 157       res = r+t;
 158       cor = (r-res)+t;
 159       if (res == res+1.0005*cor) return (m>0)?res:-res;
 160       else {
 161         res1=res+1.1*cor;
 162         z=0.5*(res1-res);
 163         __dubsin(res,z,w);
 164         z=(w[0]-ABS(x))+w[1];
 165         if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
 166         else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
 167         else {
 168           y=ABS(x);
 169           return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
 170         }
 171       }
 172     }
 173   }    /*   else  if (k < 0x3fe80000)    */
 174   /*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/
 175   else
 176   if (k < 0x3fed8000) {
 177     n = 992+((k&0x000fe000)>>13)*13;
 178     if (m>0) xx = x - asncs.x[n];
 179     else xx = -x - asncs.x[n];
 180     t = asncs.x[n+1]*xx;
 181     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
 182      +xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+xx*asncs.x[n+8]))))))+asncs.x[n+9];
 183     t+=p;
 184     res =asncs.x[n+10] +t;
 185     cor = (asncs.x[n+10]-res)+t;
 186     if (res == res+1.01*cor) return (m>0)?res:-res;
 187     else {
 188       r=asncs.x[n+10]+xx*asncs.x[n+11];
 189       t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
 190       res = r+t;
 191       cor = (r-res)+t;
 192       if (res == res+1.0008*cor) return (m>0)?res:-res;
 193       else {
 194         res1=res+1.1*cor;
 195         z=0.5*(res1-res);
 196         y=hp0.x-res;
 197         z=((hp0.x-y)-res)+(hp1.x-z);
 198         __dubcos(y,z,w);
 199         z=(w[0]-ABS(x))+w[1];
 200         if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
 201         else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
 202         else {
 203           y=ABS(x);
 204           return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
 205         }
 206       }
 207     }
 208   }    /*   else  if (k < 0x3fed8000)    */
 209   /*-------------------0.921875 <= |x| < 0.953125 ------------------------*/
 210   else
 211   if (k < 0x3fee8000) {
 212     n = 884+((k&0x000fe000)>>13)*14;
 213     if (m>0) xx = x - asncs.x[n];
 214     else xx = -x - asncs.x[n];
 215     t = asncs.x[n+1]*xx;
 216     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 217                       xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
 218                       +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
 219                       xx*asncs.x[n+9])))))))+asncs.x[n+10];
 220     t+=p;
 221     res =asncs.x[n+11] +t;
 222     cor = (asncs.x[n+11]-res)+t;
 223     if (res == res+1.01*cor) return (m>0)?res:-res;
 224     else {
 225       r=asncs.x[n+11]+xx*asncs.x[n+12];
 226       t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
 227       res = r+t;
 228       cor = (r-res)+t;
 229       if (res == res+1.0007*cor) return (m>0)?res:-res;
 230       else {
 231         res1=res+1.1*cor;
 232         z=0.5*(res1-res);
 233         y=(hp0.x-res)-z;
 234         z=y+hp1.x;
 235         y=(y-z)+hp1.x;
 236         __dubcos(z,y,w);
 237         z=(w[0]-ABS(x))+w[1];
 238         if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
 239         else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
 240         else {
 241           y=ABS(x);
 242           return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
 243         }
 244       }
 245     }
 246   }    /*   else  if (k < 0x3fee8000)    */
 247
 248   /*--------------------0.953125 <= |x| < 0.96875 ------------------------*/
 249   else
 250   if (k < 0x3fef0000) {
 251     n = 768+((k&0x000fe000)>>13)*15;
 252     if (m>0) xx = x - asncs.x[n];
 253     else xx = -x - asncs.x[n];
 254     t = asncs.x[n+1]*xx;
 255     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 256                          xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
 257                          +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
 258                     xx*(asncs.x[n+9]+xx*asncs.x[n+10]))))))))+asncs.x[n+11];
 259     t+=p;
 260     res =asncs.x[n+12] +t;
 261     cor = (asncs.x[n+12]-res)+t;
 262     if (res == res+1.01*cor) return (m>0)?res:-res;
 263     else {
 264       r=asncs.x[n+12]+xx*asncs.x[n+13];
 265       t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
 266       res = r+t;
 267       cor = (r-res)+t;
 268       if (res == res+1.0007*cor) return (m>0)?res:-res;
 269       else {
 270         res1=res+1.1*cor;
 271         z=0.5*(res1-res);
 272         y=(hp0.x-res)-z;
 273         z=y+hp1.x;
 274         y=(y-z)+hp1.x;
 275         __dubcos(z,y,w);
 276         z=(w[0]-ABS(x))+w[1];
 277         if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
 278         else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
 279         else {
 280           y=ABS(x);
 281           return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
 282         }
 283       }
 284     }
 285   }    /*   else  if (k < 0x3fef0000)    */
 286   /*--------------------0.96875 <= |x| < 1 --------------------------------*/
 287   else
 288   if (k<0x3ff00000)  {
 289     z = 0.5*((m>0)?(1.0-x):(1.0+x));
 290     v.x=z;
 291     k=v.i[HIGH_HALF];
 292     t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
 293     r=1.0-t*t*z;
 294     t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
 295     c=t*z;
 296     t=c*(1.5-0.5*t*c);
 297     y=(c+t24)-t24;
 298     cc = (z-y*y)/(t+y);
 299     p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
 300     cor = (hp1.x - 2.0*cc)-2.0*(y+cc)*p;
 301     res1 = hp0.x - 2.0*y;
 302     res =res1 + cor;
 303     if (res == res+1.003*((res1-res)+cor)) return (m>0)?res:-res;
 304     else {
 305       c=y+cc;
 306       cc=(y-c)+cc;
 307       __doasin(c,cc,w);
 308       res1=hp0.x-2.0*w[0];
 309       cor=((hp0.x-res1)-2.0*w[0])+(hp1.x-2.0*w[1]);
 310       res = res1+cor;
 311       cor = (res1-res)+cor;
 312       if (res==(res+1.0000001*cor)) return (m>0)?res:-res;
 313       else {
 314         y=ABS(x);
 315         res1=res+1.1*cor;
 316         return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
 317       }
 318     }
 319   }    /*   else  if (k < 0x3ff00000)    */
 320   /*---------------------------- |x|>=1 -------------------------------*/
 321   else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?hp0.x:-hp0.x;
 322   else
 323   if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x;
 324   else {
 325     u.i[HIGH_HALF]=0x7ff00000;
 326     v.i[HIGH_HALF]=0x7ff00000;
 327     u.i[LOW_HALF]=0;
 328     v.i[LOW_HALF]=0;
 329     return u.x/v.x;  /* NaN */
 330  }
 331 }
 332 #ifndef __ieee754_asin
 333 strong_alias (__ieee754_asin, __asin_finite)
 334 #endif
 335
 336 /*******************************************************************/
 337 /*                                                                 */
 338 /*         End of arcsine,  below is arccosine                     */
 339 /*                                                                 */
 340 /*******************************************************************/
 341
 342 double
 343 SECTION
 344 __ieee754_acos(double x)
 345 {
 346   double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2],eps;
 347 #if 0
 348   double fc;
 349 #endif
 350   mynumber u,v;
 351   int4 k,m,n;
 352 #if 0
 353   int4 nn;
 354 #endif
 355   u.x = x;
 356   m = u.i[HIGH_HALF];
 357   k = 0x7fffffff&m;
 358   /*-------------------  |x|<2.77556*10^-17 ----------------------*/
 359   if (k < 0x3c880000) return hp0.x;
 360
 361   /*-----------------  2.77556*10^-17 <= |x| < 2^-3 --------------*/
 362   else
 363   if (k < 0x3fc00000) {
 364     x2 = x*x;
 365     t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
 366     r=hp0.x-x;
 367     cor=(((hp0.x-r)-x)+hp1.x)-t;
 368     res = r+cor;
 369     cor = (r-res)+cor;
 370     if (res == res+1.004*cor) return res;
 371     else {
 372       x1 = x+big;
 373       xx = x*x;
 374       x1 -= big;
 375       x2 = x - x1;
 376       p = x1*x1*x1;
 377       s1 = a1.x*p;
 378       s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
 379             ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
 380       res1 = x+s1;
 381       s2 = ((x-res1)+s1)+s2;
 382       r=hp0.x-res1;
 383       cor=(((hp0.x-r)-res1)+hp1.x)-s2;
 384       res = r+cor;
 385       cor = (r-res)+cor;
 386       if (res == res+1.00004*cor) return res;
 387       else {
 388         __doasin(x,0,w);
 389         r=hp0.x-w[0];
 390         cor=((hp0.x-r)-w[0])+(hp1.x-w[1]);
 391         res=r+cor;
 392         cor=(r-res)+cor;
 393         if (res ==(res +1.00000001*cor)) return res;
 394         else {
 395           res1=res+1.1*cor;
 396           return __cos32(x,res,res1);
 397         }
 398       }
 399     }
 400   }    /*   else  if (k < 0x3fc00000)    */
 401   /*----------------------  0.125 <= |x| < 0.5 --------------------*/
 402   else
 403   if (k < 0x3fe00000) {
 404     if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
 405     else n = 11*((k&0x000fffff)>>14)+352;
 406     if (m>0) xx = x - asncs.x[n];
 407     else xx = -x - asncs.x[n];
 408     t = asncs.x[n+1]*xx;
 409     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 410                    xx*(asncs.x[n+5]+xx*asncs.x[n+6]))))+asncs.x[n+7];
 411     t+=p;
 412     y = (m>0)?(hp0.x-asncs.x[n+8]):(hp0.x+asncs.x[n+8]);
 413     t = (m>0)?(hp1.x-t):(hp1.x+t);
 414     res = y+t;
 415     if (res == res+1.02*((y-res)+t)) return res;
 416     else {
 417       r=asncs.x[n+8]+xx*asncs.x[n+9];
 418       t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
 419       if (m>0)
 420         {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; }
 421       else
 422         {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); }
 423       res = p+t;
 424       cor = (p-res)+t;
 425       if (res == (res+1.0002*cor)) return res;
 426       else {
 427         res1=res+1.1*cor;
 428         z=0.5*(res1-res);
 429         __docos(res,z,w);
 430         z=(w[0]-x)+w[1];
 431         if (z>1.0e-27) return max(res,res1);
 432         else if (z<-1.0e-27) return min(res,res1);
 433         else return __cos32(x,res,res1);
 434       }
 435     }
 436   }    /*   else  if (k < 0x3fe00000)    */
 437
 438   /*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/
 439   else
 440   if (k < 0x3fe80000) {
 441     n = 1056+((k&0x000fe000)>>11)*3;
 442     if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
 443     else {xx = -x - asncs.x[n]; eps=1.02; }
 444     t = asncs.x[n+1]*xx;
 445     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 446                    xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+
 447                    xx*asncs.x[n+7])))))+asncs.x[n+8];
 448     t+=p;
 449    y = (m>0)?(hp0.x-asncs.x[n+9]):(hp0.x+asncs.x[n+9]);
 450    t = (m>0)?(hp1.x-t):(hp1.x+t);
 451    res = y+t;
 452    if (res == res+eps*((y-res)+t)) return res;
 453    else {
 454      r=asncs.x[n+9]+xx*asncs.x[n+10];
 455      t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
 456      if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0004; }
 457      else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0002; }
 458      res = p+t;
 459      cor = (p-res)+t;
 460      if (res == (res+eps*cor)) return res;
 461      else {
 462        res1=res+1.1*cor;
 463        z=0.5*(res1-res);
 464        __docos(res,z,w);
 465        z=(w[0]-x)+w[1];
 466        if (z>1.0e-27) return max(res,res1);
 467        else if (z<-1.0e-27) return min(res,res1);
 468        else return __cos32(x,res,res1);
 469      }
 470    }
 471   }    /*   else  if (k < 0x3fe80000)    */
 472
 473 /*------------------------- 0.75 <= |x| < 0.921875 -------------*/
 474   else
 475   if (k < 0x3fed8000) {
 476     n = 992+((k&0x000fe000)>>13)*13;
 477     if (m>0) {xx = x - asncs.x[n]; eps = 1.04; }
 478     else {xx = -x - asncs.x[n]; eps = 1.01; }
 479     t = asncs.x[n+1]*xx;
 480     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 481                       xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+
 482                       xx*asncs.x[n+8]))))))+asncs.x[n+9];
 483     t+=p;
 484     y = (m>0)?(hp0.x-asncs.x[n+10]):(hp0.x+asncs.x[n+10]);
 485     t = (m>0)?(hp1.x-t):(hp1.x+t);
 486     res = y+t;
 487     if (res == res+eps*((y-res)+t)) return res;
 488     else {
 489       r=asncs.x[n+10]+xx*asncs.x[n+11];
 490       t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
 491       if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0032; }
 492       else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0008; }
 493       res = p+t;
 494       cor = (p-res)+t;
 495       if (res == (res+eps*cor)) return res;
 496       else {
 497         res1=res+1.1*cor;
 498         z=0.5*(res1-res);
 499         __docos(res,z,w);
 500         z=(w[0]-x)+w[1];
 501         if (z>1.0e-27) return max(res,res1);
 502         else if (z<-1.0e-27) return min(res,res1);
 503         else return __cos32(x,res,res1);
 504       }
 505     }
 506   }    /*   else  if (k < 0x3fed8000)    */
 507
 508 /*-------------------0.921875 <= |x| < 0.953125 ------------------*/
 509   else
 510   if (k < 0x3fee8000) {
 511     n = 884+((k&0x000fe000)>>13)*14;
 512     if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
 513     else {xx = -x - asncs.x[n]; eps =1.005; }
 514     t = asncs.x[n+1]*xx;
 515     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 516                    xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
 517                    +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
 518                    xx*asncs.x[n+9])))))))+asncs.x[n+10];
 519     t+=p;
 520     y = (m>0)?(hp0.x-asncs.x[n+11]):(hp0.x+asncs.x[n+11]);
 521     t = (m>0)?(hp1.x-t):(hp1.x+t);
 522     res = y+t;
 523     if (res == res+eps*((y-res)+t)) return res;
 524     else {
 525       r=asncs.x[n+11]+xx*asncs.x[n+12];
 526       t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
 527       if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
 528       else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
 529       res = p+t;
 530       cor = (p-res)+t;
 531       if (res == (res+eps*cor)) return res;
 532       else {
 533         res1=res+1.1*cor;
 534         z=0.5*(res1-res);
 535         __docos(res,z,w);
 536         z=(w[0]-x)+w[1];
 537         if (z>1.0e-27) return max(res,res1);
 538         else if (z<-1.0e-27) return min(res,res1);
 539         else return __cos32(x,res,res1);
 540       }
 541     }
 542   }    /*   else  if (k < 0x3fee8000)    */
 543
 544   /*--------------------0.953125 <= |x| < 0.96875 ----------------*/
 545   else
 546   if (k < 0x3fef0000) {
 547     n = 768+((k&0x000fe000)>>13)*15;
 548     if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
 549     else {xx = -x - asncs.x[n]; eps=1.005;}
 550     t = asncs.x[n+1]*xx;
 551     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 552             xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
 553             +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+xx*(asncs.x[n+9]+
 554             xx*asncs.x[n+10]))))))))+asncs.x[n+11];
 555     t+=p;
 556     y = (m>0)?(hp0.x-asncs.x[n+12]):(hp0.x+asncs.x[n+12]);
 557    t = (m>0)?(hp1.x-t):(hp1.x+t);
 558    res = y+t;
 559    if (res == res+eps*((y-res)+t)) return res;
 560    else {
 561      r=asncs.x[n+12]+xx*asncs.x[n+13];
 562      t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
 563      if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
 564      else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
 565      res = p+t;
 566      cor = (p-res)+t;
 567      if (res == (res+eps*cor)) return res;
 568      else {
 569        res1=res+1.1*cor;
 570        z=0.5*(res1-res);
 571        __docos(res,z,w);
 572        z=(w[0]-x)+w[1];
 573        if (z>1.0e-27) return max(res,res1);
 574        else if (z<-1.0e-27) return min(res,res1);
 575        else return __cos32(x,res,res1);
 576      }
 577    }
 578   }    /*   else  if (k < 0x3fef0000)    */
 579   /*-----------------0.96875 <= |x| < 1 ---------------------------*/
 580
 581   else
 582   if (k<0x3ff00000)  {
 583     z = 0.5*((m>0)?(1.0-x):(1.0+x));
 584     v.x=z;
 585     k=v.i[HIGH_HALF];
 586     t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
 587     r=1.0-t*t*z;
 588     t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
 589     c=t*z;
 590     t=c*(1.5-0.5*t*c);
 591     y = (t27*c+c)-t27*c;
 592     cc = (z-y*y)/(t+y);
 593     p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
 594     if (m<0) {
 595       cor = (hp1.x - cc)-(y+cc)*p;
 596       res1 = hp0.x - y;
 597       res =res1 + cor;
 598       if (res == res+1.002*((res1-res)+cor)) return (res+res);
 599       else {
 600         c=y+cc;
 601         cc=(y-c)+cc;
 602         __doasin(c,cc,w);
 603         res1=hp0.x-w[0];
 604         cor=((hp0.x-res1)-w[0])+(hp1.x-w[1]);
 605         res = res1+cor;
 606         cor = (res1-res)+cor;
 607         if (res==(res+1.000001*cor)) return (res+res);
 608         else {
 609           res=res+res;
 610           res1=res+1.2*cor;
 611           return __cos32(x,res,res1);
 612         }
 613       }
 614     }
 615     else {
 616       cor = cc+p*(y+cc);
 617       res = y + cor;
 618       if (res == res+1.03*((y-res)+cor)) return (res+res);
 619       else {
 620         c=y+cc;
 621         cc=(y-c)+cc;
 622         __doasin(c,cc,w);
 623         res = w[0];
 624         cor=w[1];
 625         if (res==(res+1.000001*cor)) return (res+res);
 626         else {
 627           res=res+res;
 628           res1=res+1.2*cor;
 629           return __cos32(x,res,res1);
 630         }
 631       }
 632     }
 633   }    /*   else  if (k < 0x3ff00000)    */
 634
 635   /*---------------------------- |x|>=1 -----------------------*/
 636   else
 637   if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?0:2.0*hp0.x;
 638   else
 639   if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x;
 640   else {
 641     u.i[HIGH_HALF]=0x7ff00000;
 642     v.i[HIGH_HALF]=0x7ff00000;
 643     u.i[LOW_HALF]=0;
 644     v.i[LOW_HALF]=0;
 645     return u.x/v.x;
 646   }
 647 }
 648 #ifndef __ieee754_acos
 649 strong_alias (__ieee754_acos, __acos_finite)
 650 #endif