libgcc-math/dbl-64/e_asin.c

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