sysdeps/ieee754/dbl-64/e_asin.c

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