sysdeps/ieee754/dbl-64/e_asin.c

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