sysdeps/ieee754/dbl-64/e_asin.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2021 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 <https://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 /******************************************************************/
  29 #include "endian.h"
  30 #include "mydefs.h"
  31 #include "asincos.tbl"
  32 #include "root.tbl"
  33 #include "powtwo.tbl"
  34 #include "MathLib.h"
  35 #include "uasncs.h"
  36 #include <float.h>
  37 #include <math.h>
  38 #include <math_private.h>
  39 #include <math-underflow.h>
  40 #include <libm-alias-finite.h>
  41
  42 #ifndef SECTION
  43 # define SECTION
  44 #endif
  45
  46 void __doasin(double x, double dx, double w[]);
  47 void __dubsin(double x, double dx, double v[]);
  48 void __dubcos(double x, double dx, double v[]);
  49 void __docos(double x, double dx, double v[]);
  50
  51 double
  52 SECTION
  53 __ieee754_asin(double x){
  54   double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2];
  55   mynumber u,v;
  56   int4 k,m,n;
  57
  58   u.x = x;
  59   m = u.i[HIGH_HALF];
  60   k = 0x7fffffff&m;              /* no sign */
  61
  62   if (k < 0x3e500000)
  63     {
  64       math_check_force_underflow (x);
  65       return x;  /* for x->0 => sin(x)=x */
  66     }
  67   /*----------------------2^-26 <= |x| < 2^ -3    -----------------*/
  68   else
  69   if (k < 0x3fc00000) {
  70     x2 = x*x;
  71     t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
  72     res = x+t;         /*  res=arcsin(x) according to Taylor series  */
  73     cor = (x-res)+t;
  74     if (res == res+1.025*cor) return res;
  75     else {
  76       x1 = x+big;
  77       xx = x*x;
  78       x1 -= big;
  79       x2 = x - x1;
  80       p = x1*x1*x1;
  81       s1 = a1.x*p;
  82       s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
  83              ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
  84       res1 = x+s1;
  85       s2 = ((x-res1)+s1)+s2;
  86       res = res1+s2;
  87       cor = (res1-res)+s2;
  88       if (res == res+1.00014*cor) return res;
  89       else {
  90         __doasin(x,0,w);
  91         return w[0];
  92       }
  93     }
  94   }
  95   /*---------------------0.125 <= |x| < 0.5 -----------------------------*/
  96   else if (k < 0x3fe00000) {
  97     if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
  98     else n = 11*((k&0x000fffff)>>14)+352;
  99     if (m>0) xx = x - asncs.x[n];
 100     else xx = -x - asncs.x[n];
 101     t = asncs.x[n+1]*xx;
 102     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
 103      +xx*asncs.x[n+6]))))+asncs.x[n+7];
 104     t+=p;
 105     res =asncs.x[n+8] +t;
 106     cor = (asncs.x[n+8]-res)+t;
 107     if (res == res+1.05*cor) return (m>0)?res:-res;
 108     else {
 109       r=asncs.x[n+8]+xx*asncs.x[n+9];
 110       t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
 111       res = r+t;
 112       cor = (r-res)+t;
 113       if (res == res+1.0005*cor) return (m>0)?res:-res;
 114       else {
 115         res1=res+1.1*cor;
 116         z=0.5*(res1-res);
 117         __dubsin(res,z,w);
 118         z=(w[0]-fabs(x))+w[1];
 119         if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
 120         else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
 121         else {
 122           return (m>0)?res:-res;
 123         }
 124       }
 125     }
 126   }    /*   else  if (k < 0x3fe00000)    */
 127   /*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/
 128   else
 129   if (k < 0x3fe80000) {
 130     n = 1056+((k&0x000fe000)>>11)*3;
 131     if (m>0) xx = x - asncs.x[n];
 132     else xx = -x - asncs.x[n];
 133     t = asncs.x[n+1]*xx;
 134     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
 135            +xx*(asncs.x[n+6]+xx*asncs.x[n+7])))))+asncs.x[n+8];
 136     t+=p;
 137     res =asncs.x[n+9] +t;
 138     cor = (asncs.x[n+9]-res)+t;
 139     if (res == res+1.01*cor) return (m>0)?res:-res;
 140     else {
 141       r=asncs.x[n+9]+xx*asncs.x[n+10];
 142       t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
 143       res = r+t;
 144       cor = (r-res)+t;
 145       if (res == res+1.0005*cor) return (m>0)?res:-res;
 146       else {
 147         res1=res+1.1*cor;
 148         z=0.5*(res1-res);
 149         __dubsin(res,z,w);
 150         z=(w[0]-fabs(x))+w[1];
 151         if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
 152         else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
 153         else {
 154           return (m>0)?res:-res;
 155         }
 156       }
 157     }
 158   }    /*   else  if (k < 0x3fe80000)    */
 159   /*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/
 160   else
 161   if (k < 0x3fed8000) {
 162     n = 992+((k&0x000fe000)>>13)*13;
 163     if (m>0) xx = x - asncs.x[n];
 164     else xx = -x - asncs.x[n];
 165     t = asncs.x[n+1]*xx;
 166     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
 167      +xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+xx*asncs.x[n+8]))))))+asncs.x[n+9];
 168     t+=p;
 169     res =asncs.x[n+10] +t;
 170     cor = (asncs.x[n+10]-res)+t;
 171     if (res == res+1.01*cor) return (m>0)?res:-res;
 172     else {
 173       r=asncs.x[n+10]+xx*asncs.x[n+11];
 174       t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
 175       res = r+t;
 176       cor = (r-res)+t;
 177       if (res == res+1.0008*cor) return (m>0)?res:-res;
 178       else {
 179         res1=res+1.1*cor;
 180         z=0.5*(res1-res);
 181         y=hp0.x-res;
 182         z=((hp0.x-y)-res)+(hp1.x-z);
 183         __dubcos(y,z,w);
 184         z=(w[0]-fabs(x))+w[1];
 185         if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
 186         else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
 187         else {
 188           return (m>0)?res:-res;
 189         }
 190       }
 191     }
 192   }    /*   else  if (k < 0x3fed8000)    */
 193   /*-------------------0.921875 <= |x| < 0.953125 ------------------------*/
 194   else
 195   if (k < 0x3fee8000) {
 196     n = 884+((k&0x000fe000)>>13)*14;
 197     if (m>0) xx = x - asncs.x[n];
 198     else xx = -x - asncs.x[n];
 199     t = asncs.x[n+1]*xx;
 200     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 201                       xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
 202                       +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
 203                       xx*asncs.x[n+9])))))))+asncs.x[n+10];
 204     t+=p;
 205     res =asncs.x[n+11] +t;
 206     cor = (asncs.x[n+11]-res)+t;
 207     if (res == res+1.01*cor) return (m>0)?res:-res;
 208     else {
 209       r=asncs.x[n+11]+xx*asncs.x[n+12];
 210       t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
 211       res = r+t;
 212       cor = (r-res)+t;
 213       if (res == res+1.0007*cor) return (m>0)?res:-res;
 214       else {
 215         res1=res+1.1*cor;
 216         z=0.5*(res1-res);
 217         y=(hp0.x-res)-z;
 218         z=y+hp1.x;
 219         y=(y-z)+hp1.x;
 220         __dubcos(z,y,w);
 221         z=(w[0]-fabs(x))+w[1];
 222         if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
 223         else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
 224         else {
 225           return (m>0)?res:-res;
 226         }
 227       }
 228     }
 229   }    /*   else  if (k < 0x3fee8000)    */
 230
 231   /*--------------------0.953125 <= |x| < 0.96875 ------------------------*/
 232   else
 233   if (k < 0x3fef0000) {
 234     n = 768+((k&0x000fe000)>>13)*15;
 235     if (m>0) xx = x - asncs.x[n];
 236     else xx = -x - asncs.x[n];
 237     t = asncs.x[n+1]*xx;
 238     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 239                          xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
 240                          +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
 241                     xx*(asncs.x[n+9]+xx*asncs.x[n+10]))))))))+asncs.x[n+11];
 242     t+=p;
 243     res =asncs.x[n+12] +t;
 244     cor = (asncs.x[n+12]-res)+t;
 245     if (res == res+1.01*cor) return (m>0)?res:-res;
 246     else {
 247       r=asncs.x[n+12]+xx*asncs.x[n+13];
 248       t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
 249       res = r+t;
 250       cor = (r-res)+t;
 251       if (res == res+1.0007*cor) return (m>0)?res:-res;
 252       else {
 253         res1=res+1.1*cor;
 254         z=0.5*(res1-res);
 255         y=(hp0.x-res)-z;
 256         z=y+hp1.x;
 257         y=(y-z)+hp1.x;
 258         __dubcos(z,y,w);
 259         z=(w[0]-fabs(x))+w[1];
 260         if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
 261         else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
 262         else {
 263           return (m>0)?res:-res;
 264         }
 265       }
 266     }
 267   }    /*   else  if (k < 0x3fef0000)    */
 268   /*--------------------0.96875 <= |x| < 1 --------------------------------*/
 269   else
 270   if (k<0x3ff00000)  {
 271     z = 0.5*((m>0)?(1.0-x):(1.0+x));
 272     v.x=z;
 273     k=v.i[HIGH_HALF];
 274     t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
 275     r=1.0-t*t*z;
 276     t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
 277     c=t*z;
 278     t=c*(1.5-0.5*t*c);
 279     y=(c+t24)-t24;
 280     cc = (z-y*y)/(t+y);
 281     p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
 282     cor = (hp1.x - 2.0*cc)-2.0*(y+cc)*p;
 283     res1 = hp0.x - 2.0*y;
 284     res =res1 + cor;
 285     if (res == res+1.003*((res1-res)+cor)) return (m>0)?res:-res;
 286     else {
 287       c=y+cc;
 288       cc=(y-c)+cc;
 289       __doasin(c,cc,w);
 290       res1=hp0.x-2.0*w[0];
 291       cor=((hp0.x-res1)-2.0*w[0])+(hp1.x-2.0*w[1]);
 292       res = res1+cor;
 293       return (m>0)?res:-res;
 294     }
 295   }    /*   else  if (k < 0x3ff00000)    */
 296   /*---------------------------- |x|>=1 -------------------------------*/
 297   else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?hp0.x:-hp0.x;
 298   else
 299   if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x + x;
 300   else {
 301     u.i[HIGH_HALF]=0x7ff00000;
 302     v.i[HIGH_HALF]=0x7ff00000;
 303     u.i[LOW_HALF]=0;
 304     v.i[LOW_HALF]=0;
 305     return u.x/v.x;  /* NaN */
 306  }
 307 }
 308 #ifndef __ieee754_asin
 309 libm_alias_finite (__ieee754_asin, __asin)
 310 #endif
 311
 312 /*******************************************************************/
 313 /*                                                                 */
 314 /*         End of arcsine,  below is arccosine                     */
 315 /*                                                                 */
 316 /*******************************************************************/
 317
 318 double
 319 SECTION
 320 __ieee754_acos(double x)
 321 {
 322   double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2],eps;
 323   mynumber u,v;
 324   int4 k,m,n;
 325   u.x = x;
 326   m = u.i[HIGH_HALF];
 327   k = 0x7fffffff&m;
 328   /*-------------------  |x|<2.77556*10^-17 ----------------------*/
 329   if (k < 0x3c880000) return hp0.x;
 330
 331   /*-----------------  2.77556*10^-17 <= |x| < 2^-3 --------------*/
 332   else
 333   if (k < 0x3fc00000) {
 334     x2 = x*x;
 335     t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
 336     r=hp0.x-x;
 337     cor=(((hp0.x-r)-x)+hp1.x)-t;
 338     res = r+cor;
 339     cor = (r-res)+cor;
 340     if (res == res+1.004*cor) return res;
 341     else {
 342       x1 = x+big;
 343       xx = x*x;
 344       x1 -= big;
 345       x2 = x - x1;
 346       p = x1*x1*x1;
 347       s1 = a1.x*p;
 348       s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
 349             ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
 350       res1 = x+s1;
 351       s2 = ((x-res1)+s1)+s2;
 352       r=hp0.x-res1;
 353       cor=(((hp0.x-r)-res1)+hp1.x)-s2;
 354       res = r+cor;
 355       cor = (r-res)+cor;
 356       if (res == res+1.00004*cor) return res;
 357       else {
 358         __doasin(x,0,w);
 359         r=hp0.x-w[0];
 360         cor=((hp0.x-r)-w[0])+(hp1.x-w[1]);
 361         res=r+cor;
 362         return res;
 363       }
 364     }
 365   }    /*   else  if (k < 0x3fc00000)    */
 366   /*----------------------  0.125 <= |x| < 0.5 --------------------*/
 367   else
 368   if (k < 0x3fe00000) {
 369     if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
 370     else n = 11*((k&0x000fffff)>>14)+352;
 371     if (m>0) xx = x - asncs.x[n];
 372     else xx = -x - asncs.x[n];
 373     t = asncs.x[n+1]*xx;
 374     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 375                    xx*(asncs.x[n+5]+xx*asncs.x[n+6]))))+asncs.x[n+7];
 376     t+=p;
 377     y = (m>0)?(hp0.x-asncs.x[n+8]):(hp0.x+asncs.x[n+8]);
 378     t = (m>0)?(hp1.x-t):(hp1.x+t);
 379     res = y+t;
 380     if (res == res+1.02*((y-res)+t)) return res;
 381     else {
 382       r=asncs.x[n+8]+xx*asncs.x[n+9];
 383       t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
 384       if (m>0)
 385         {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; }
 386       else
 387         {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); }
 388       res = p+t;
 389       cor = (p-res)+t;
 390       if (res == (res+1.0002*cor)) return res;
 391       else {
 392         res1=res+1.1*cor;
 393         z=0.5*(res1-res);
 394         __docos(res,z,w);
 395         z=(w[0]-x)+w[1];
 396         if (z>1.0e-27) return max(res,res1);
 397         else if (z<-1.0e-27) return min(res,res1);
 398         else return res;
 399       }
 400     }
 401   }    /*   else  if (k < 0x3fe00000)    */
 402
 403   /*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/
 404   else
 405   if (k < 0x3fe80000) {
 406     n = 1056+((k&0x000fe000)>>11)*3;
 407     if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
 408     else {xx = -x - asncs.x[n]; eps=1.02; }
 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]+
 412                    xx*asncs.x[n+7])))))+asncs.x[n+8];
 413     t+=p;
 414    y = (m>0)?(hp0.x-asncs.x[n+9]):(hp0.x+asncs.x[n+9]);
 415    t = (m>0)?(hp1.x-t):(hp1.x+t);
 416    res = y+t;
 417    if (res == res+eps*((y-res)+t)) return res;
 418    else {
 419      r=asncs.x[n+9]+xx*asncs.x[n+10];
 420      t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
 421      if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0004; }
 422      else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0002; }
 423      res = p+t;
 424      cor = (p-res)+t;
 425      if (res == (res+eps*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 res;
 434      }
 435    }
 436   }    /*   else  if (k < 0x3fe80000)    */
 437
 438 /*------------------------- 0.75 <= |x| < 0.921875 -------------*/
 439   else
 440   if (k < 0x3fed8000) {
 441     n = 992+((k&0x000fe000)>>13)*13;
 442     if (m>0) {xx = x - asncs.x[n]; eps = 1.04; }
 443     else {xx = -x - asncs.x[n]; eps = 1.01; }
 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]+xx*(asncs.x[n+7]+
 447                       xx*asncs.x[n+8]))))))+asncs.x[n+9];
 448     t+=p;
 449     y = (m>0)?(hp0.x-asncs.x[n+10]):(hp0.x+asncs.x[n+10]);
 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+10]+xx*asncs.x[n+11];
 455       t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
 456       if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0032; }
 457       else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0008; }
 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 res;
 469       }
 470     }
 471   }    /*   else  if (k < 0x3fed8000)    */
 472
 473 /*-------------------0.921875 <= |x| < 0.953125 ------------------*/
 474   else
 475   if (k < 0x3fee8000) {
 476     n = 884+((k&0x000fe000)>>13)*14;
 477     if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
 478     else {xx = -x - asncs.x[n]; eps =1.005; }
 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]
 482                    +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
 483                    xx*asncs.x[n+9])))))))+asncs.x[n+10];
 484     t+=p;
 485     y = (m>0)?(hp0.x-asncs.x[n+11]):(hp0.x+asncs.x[n+11]);
 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+11]+xx*asncs.x[n+12];
 491       t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
 492       if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
 493       else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
 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 res;
 505       }
 506     }
 507   }    /*   else  if (k < 0x3fee8000)    */
 508
 509   /*--------------------0.953125 <= |x| < 0.96875 ----------------*/
 510   else
 511   if (k < 0x3fef0000) {
 512     n = 768+((k&0x000fe000)>>13)*15;
 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]+xx*(asncs.x[n+9]+
 519             xx*asncs.x[n+10]))))))))+asncs.x[n+11];
 520     t+=p;
 521     y = (m>0)?(hp0.x-asncs.x[n+12]):(hp0.x+asncs.x[n+12]);
 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+12]+xx*asncs.x[n+13];
 527      t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
 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 res;
 541      }
 542    }
 543   }    /*   else  if (k < 0x3fef0000)    */
 544   /*-----------------0.96875 <= |x| < 1 ---------------------------*/
 545
 546   else
 547   if (k<0x3ff00000)  {
 548     z = 0.5*((m>0)?(1.0-x):(1.0+x));
 549     v.x=z;
 550     k=v.i[HIGH_HALF];
 551     t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
 552     r=1.0-t*t*z;
 553     t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
 554     c=t*z;
 555     t=c*(1.5-0.5*t*c);
 556     y = (t27*c+c)-t27*c;
 557     cc = (z-y*y)/(t+y);
 558     p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
 559     if (m<0) {
 560       cor = (hp1.x - cc)-(y+cc)*p;
 561       res1 = hp0.x - y;
 562       res =res1 + cor;
 563       if (res == res+1.002*((res1-res)+cor)) return (res+res);
 564       else {
 565         c=y+cc;
 566         cc=(y-c)+cc;
 567         __doasin(c,cc,w);
 568         res1=hp0.x-w[0];
 569         cor=((hp0.x-res1)-w[0])+(hp1.x-w[1]);
 570         res = res1+cor;
 571         return (res+res);
 572       }
 573     }
 574     else {
 575       cor = cc+p*(y+cc);
 576       res = y + cor;
 577       if (res == res+1.03*((y-res)+cor)) return (res+res);
 578       else {
 579         c=y+cc;
 580         cc=(y-c)+cc;
 581         __doasin(c,cc,w);
 582         res = w[0];
 583         return (res+res);
 584       }
 585     }
 586   }    /*   else  if (k < 0x3ff00000)    */
 587
 588   /*---------------------------- |x|>=1 -----------------------*/
 589   else
 590   if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?0:2.0*hp0.x;
 591   else
 592   if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x + x;
 593   else {
 594     u.i[HIGH_HALF]=0x7ff00000;
 595     v.i[HIGH_HALF]=0x7ff00000;
 596     u.i[LOW_HALF]=0;
 597     v.i[LOW_HALF]=0;
 598     return u.x/v.x;
 599   }
 600 }
 601 #ifndef __ieee754_acos
 602 libm_alias_finite (__ieee754_acos, __acos)
 603 #endif