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 45void __doasin(double x, double dx, double w[]); 46void __dubsin(double x, double dx, double v[]); 47void __dubcos(double x, double dx, double v[]); 48void __docos(double x, double dx, double v[]); 49double __sin32(double x, double res, double res1); 50double __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/***************************************************************************/ 56double __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 334double __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} 638