1 /* ef_j0.c -- float version of e_j0.c. 2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 3 */ 4 5 /* 6 * ==================================================== 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 8 * 9 * Developed at SunPro, a Sun Microsystems, Inc. business. 10 * Permission to use, copy, modify, and distribute this 11 * software is freely granted, provided that this notice 12 * is preserved. 13 * ==================================================== 14 */ 15 16 #include "fdlibm.h" 17 18 #ifdef __STDC__ 19 static float pzerof(float), qzerof(float); 20 #else 21 static float pzerof(), qzerof(); 22 #endif 23 24 #ifdef __STDC__ 25 static const float 26 #else 27 static float 28 #endif 29 huge = 1e30, 30 one = 1.0, 31 invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */ 32 tpi = 6.3661974669e-01, /* 0x3f22f983 */ 33 /* R0/S0 on [0, 2.00] */ 34 R02 = 1.5625000000e-02, /* 0x3c800000 */ 35 R03 = -1.8997929874e-04, /* 0xb947352e */ 36 R04 = 1.8295404516e-06, /* 0x35f58e88 */ 37 R05 = -4.6183270541e-09, /* 0xb19eaf3c */ 38 S01 = 1.5619102865e-02, /* 0x3c7fe744 */ 39 S02 = 1.1692678527e-04, /* 0x38f53697 */ 40 S03 = 5.1354652442e-07, /* 0x3509daa6 */ 41 S04 = 1.1661400734e-09; /* 0x30a045e8 */ 42 43 #ifdef __STDC__ 44 static const float zero = 0.0; 45 #else 46 static float zero = 0.0; 47 #endif 48 49 #ifdef __STDC__ __ieee754_j0f(float x)50 float __ieee754_j0f(float x) 51 #else 52 float __ieee754_j0f(x) 53 float x; 54 #endif 55 { 56 float z, s,c,ss,cc,r,u,v; 57 __int32_t hx,ix; 58 59 GET_FLOAT_WORD(hx,x); 60 ix = hx&0x7fffffff; 61 if(!FLT_UWORD_IS_FINITE(ix)) return one/(x*x); 62 x = fabsf(x); 63 if(ix >= 0x40000000) { /* |x| >= 2.0 */ 64 s = sinf(x); 65 c = cosf(x); 66 ss = s-c; 67 cc = s+c; 68 if(ix<=FLT_UWORD_HALF_MAX) { /* make sure x+x not overflow */ 69 z = -cosf(x+x); 70 if ((s*c)<zero) cc = z/ss; 71 else ss = z/cc; 72 } 73 /* 74 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) 75 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) 76 */ 77 if(ix>0x80000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(x); 78 else { 79 u = pzerof(x); v = qzerof(x); 80 z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrtf(x); 81 } 82 return z; 83 } 84 if(ix<0x39000000) { /* |x| < 2**-13 */ 85 if(huge+x>one) { /* raise inexact if x != 0 */ 86 if(ix<0x32000000) return one; /* |x|<2**-27 */ 87 else return one - (float)0.25*x*x; 88 } 89 } 90 z = x*x; 91 r = z*(R02+z*(R03+z*(R04+z*R05))); 92 s = one+z*(S01+z*(S02+z*(S03+z*S04))); 93 if(ix < 0x3F800000) { /* |x| < 1.00 */ 94 return one + z*((float)-0.25+(r/s)); 95 } else { 96 u = (float)0.5*x; 97 return((one+u)*(one-u)+z*(r/s)); 98 } 99 } 100 101 #ifdef __STDC__ 102 static const float 103 #else 104 static float 105 #endif 106 u00 = -7.3804296553e-02, /* 0xbd9726b5 */ 107 u01 = 1.7666645348e-01, /* 0x3e34e80d */ 108 u02 = -1.3818567619e-02, /* 0xbc626746 */ 109 u03 = 3.4745343146e-04, /* 0x39b62a69 */ 110 u04 = -3.8140706238e-06, /* 0xb67ff53c */ 111 u05 = 1.9559013964e-08, /* 0x32a802ba */ 112 u06 = -3.9820518410e-11, /* 0xae2f21eb */ 113 v01 = 1.2730483897e-02, /* 0x3c509385 */ 114 v02 = 7.6006865129e-05, /* 0x389f65e0 */ 115 v03 = 2.5915085189e-07, /* 0x348b216c */ 116 v04 = 4.4111031494e-10; /* 0x2ff280c2 */ 117 118 #ifdef __STDC__ __ieee754_y0f(float x)119 float __ieee754_y0f(float x) 120 #else 121 float __ieee754_y0f(x) 122 float x; 123 #endif 124 { 125 float z, s,c,ss,cc,u,v; 126 __int32_t hx,ix; 127 128 GET_FLOAT_WORD(hx,x); 129 ix = 0x7fffffff&hx; 130 /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ 131 if(!FLT_UWORD_IS_FINITE(ix)) return one/(x+x*x); 132 if(FLT_UWORD_IS_ZERO(ix)) return -one/zero; 133 if(hx<0) return zero/zero; 134 if(ix >= 0x40000000) { /* |x| >= 2.0 */ 135 /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) 136 * where x0 = x-pi/4 137 * Better formula: 138 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) 139 * = 1/sqrt(2) * (sin(x) + cos(x)) 140 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) 141 * = 1/sqrt(2) * (sin(x) - cos(x)) 142 * To avoid cancellation, use 143 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) 144 * to compute the worse one. 145 */ 146 s = sinf(x); 147 c = cosf(x); 148 ss = s-c; 149 cc = s+c; 150 /* 151 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) 152 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) 153 */ 154 if(ix<=FLT_UWORD_HALF_MAX) { /* make sure x+x not overflow */ 155 z = -cosf(x+x); 156 if ((s*c)<zero) cc = z/ss; 157 else ss = z/cc; 158 } 159 if(ix>0x80000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x); 160 else { 161 u = pzerof(x); v = qzerof(x); 162 z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x); 163 } 164 return z; 165 } 166 if(ix<=0x32000000) { /* x < 2**-27 */ 167 return(u00 + tpi*__ieee754_logf(x)); 168 } 169 z = x*x; 170 u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); 171 v = one+z*(v01+z*(v02+z*(v03+z*v04))); 172 return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x))); 173 } 174 175 /* The asymptotic expansions of pzero is 176 * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. 177 * For x >= 2, We approximate pzero by 178 * pzero(x) = 1 + (R/S) 179 * where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10 180 * S = 1 + pS0*s^2 + ... + pS4*s^10 181 * and 182 * | pzero(x)-1-R/S | <= 2 ** ( -60.26) 183 */ 184 #ifdef __STDC__ 185 static const float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 186 #else 187 static float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 188 #endif 189 0.0000000000e+00, /* 0x00000000 */ 190 -7.0312500000e-02, /* 0xbd900000 */ 191 -8.0816707611e+00, /* 0xc1014e86 */ 192 -2.5706311035e+02, /* 0xc3808814 */ 193 -2.4852163086e+03, /* 0xc51b5376 */ 194 -5.2530439453e+03, /* 0xc5a4285a */ 195 }; 196 #ifdef __STDC__ 197 static const float pS8[5] = { 198 #else 199 static float pS8[5] = { 200 #endif 201 1.1653436279e+02, /* 0x42e91198 */ 202 3.8337448730e+03, /* 0x456f9beb */ 203 4.0597855469e+04, /* 0x471e95db */ 204 1.1675296875e+05, /* 0x47e4087c */ 205 4.7627726562e+04, /* 0x473a0bba */ 206 }; 207 #ifdef __STDC__ 208 static const float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 209 #else 210 static float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 211 #endif 212 -1.1412546255e-11, /* 0xad48c58a */ 213 -7.0312492549e-02, /* 0xbd8fffff */ 214 -4.1596107483e+00, /* 0xc0851b88 */ 215 -6.7674766541e+01, /* 0xc287597b */ 216 -3.3123129272e+02, /* 0xc3a59d9b */ 217 -3.4643338013e+02, /* 0xc3ad3779 */ 218 }; 219 #ifdef __STDC__ 220 static const float pS5[5] = { 221 #else 222 static float pS5[5] = { 223 #endif 224 6.0753936768e+01, /* 0x42730408 */ 225 1.0512523193e+03, /* 0x44836813 */ 226 5.9789707031e+03, /* 0x45bad7c4 */ 227 9.6254453125e+03, /* 0x461665c8 */ 228 2.4060581055e+03, /* 0x451660ee */ 229 }; 230 231 #ifdef __STDC__ 232 static const float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 233 #else 234 static float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 235 #endif 236 -2.5470459075e-09, /* 0xb12f081b */ 237 -7.0311963558e-02, /* 0xbd8fffb8 */ 238 -2.4090321064e+00, /* 0xc01a2d95 */ 239 -2.1965976715e+01, /* 0xc1afba52 */ 240 -5.8079170227e+01, /* 0xc2685112 */ 241 -3.1447946548e+01, /* 0xc1fb9565 */ 242 }; 243 #ifdef __STDC__ 244 static const float pS3[5] = { 245 #else 246 static float pS3[5] = { 247 #endif 248 3.5856033325e+01, /* 0x420f6c94 */ 249 3.6151397705e+02, /* 0x43b4c1ca */ 250 1.1936077881e+03, /* 0x44953373 */ 251 1.1279968262e+03, /* 0x448cffe6 */ 252 1.7358093262e+02, /* 0x432d94b8 */ 253 }; 254 255 #ifdef __STDC__ 256 static const float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 257 #else 258 static float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 259 #endif 260 -8.8753431271e-08, /* 0xb3be98b7 */ 261 -7.0303097367e-02, /* 0xbd8ffb12 */ 262 -1.4507384300e+00, /* 0xbfb9b1cc */ 263 -7.6356959343e+00, /* 0xc0f4579f */ 264 -1.1193166733e+01, /* 0xc1331736 */ 265 -3.2336456776e+00, /* 0xc04ef40d */ 266 }; 267 #ifdef __STDC__ 268 static const float pS2[5] = { 269 #else 270 static float pS2[5] = { 271 #endif 272 2.2220300674e+01, /* 0x41b1c32d */ 273 1.3620678711e+02, /* 0x430834f0 */ 274 2.7047027588e+02, /* 0x43873c32 */ 275 1.5387539673e+02, /* 0x4319e01a */ 276 1.4657617569e+01, /* 0x416a859a */ 277 }; 278 279 #ifdef __STDC__ pzerof(float x)280 static float pzerof(float x) 281 #else 282 static float pzerof(x) 283 float x; 284 #endif 285 { 286 #ifdef __STDC__ 287 const float *p,*q; 288 #else 289 float *p,*q; 290 #endif 291 float z,r,s; 292 __int32_t ix; 293 GET_FLOAT_WORD(ix,x); 294 ix &= 0x7fffffff; 295 if(ix>=0x41000000) {p = pR8; q= pS8;} 296 else if(ix>=0x40f71c58){p = pR5; q= pS5;} 297 else if(ix>=0x4036db68){p = pR3; q= pS3;} 298 else {p = pR2; q= pS2;} 299 z = one/(x*x); 300 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 301 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); 302 return one+ r/s; 303 } 304 305 306 /* For x >= 8, the asymptotic expansions of qzero is 307 * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. 308 * We approximate qzero by 309 * qzero(x) = s*(-1.25 + (R/S)) 310 * where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10 311 * S = 1 + qS0*s^2 + ... + qS5*s^12 312 * and 313 * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) 314 */ 315 #ifdef __STDC__ 316 static const float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 317 #else 318 static float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 319 #endif 320 0.0000000000e+00, /* 0x00000000 */ 321 7.3242187500e-02, /* 0x3d960000 */ 322 1.1768206596e+01, /* 0x413c4a93 */ 323 5.5767340088e+02, /* 0x440b6b19 */ 324 8.8591972656e+03, /* 0x460a6cca */ 325 3.7014625000e+04, /* 0x471096a0 */ 326 }; 327 #ifdef __STDC__ 328 static const float qS8[6] = { 329 #else 330 static float qS8[6] = { 331 #endif 332 1.6377603149e+02, /* 0x4323c6aa */ 333 8.0983447266e+03, /* 0x45fd12c2 */ 334 1.4253829688e+05, /* 0x480b3293 */ 335 8.0330925000e+05, /* 0x49441ed4 */ 336 8.4050156250e+05, /* 0x494d3359 */ 337 -3.4389928125e+05, /* 0xc8a7eb69 */ 338 }; 339 340 #ifdef __STDC__ 341 static const float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 342 #else 343 static float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 344 #endif 345 1.8408595828e-11, /* 0x2da1ec79 */ 346 7.3242180049e-02, /* 0x3d95ffff */ 347 5.8356351852e+00, /* 0x40babd86 */ 348 1.3511157227e+02, /* 0x43071c90 */ 349 1.0272437744e+03, /* 0x448067cd */ 350 1.9899779053e+03, /* 0x44f8bf4b */ 351 }; 352 #ifdef __STDC__ 353 static const float qS5[6] = { 354 #else 355 static float qS5[6] = { 356 #endif 357 8.2776611328e+01, /* 0x42a58da0 */ 358 2.0778142090e+03, /* 0x4501dd07 */ 359 1.8847289062e+04, /* 0x46933e94 */ 360 5.6751113281e+04, /* 0x475daf1d */ 361 3.5976753906e+04, /* 0x470c88c1 */ 362 -5.3543427734e+03, /* 0xc5a752be */ 363 }; 364 365 #ifdef __STDC__ 366 static const float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 367 #else 368 static float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 369 #endif 370 4.3774099900e-09, /* 0x3196681b */ 371 7.3241114616e-02, /* 0x3d95ff70 */ 372 3.3442313671e+00, /* 0x405607e3 */ 373 4.2621845245e+01, /* 0x422a7cc5 */ 374 1.7080809021e+02, /* 0x432acedf */ 375 1.6673394775e+02, /* 0x4326bbe4 */ 376 }; 377 #ifdef __STDC__ 378 static const float qS3[6] = { 379 #else 380 static float qS3[6] = { 381 #endif 382 4.8758872986e+01, /* 0x42430916 */ 383 7.0968920898e+02, /* 0x44316c1c */ 384 3.7041481934e+03, /* 0x4567825f */ 385 6.4604252930e+03, /* 0x45c9e367 */ 386 2.5163337402e+03, /* 0x451d4557 */ 387 -1.4924745178e+02, /* 0xc3153f59 */ 388 }; 389 390 #ifdef __STDC__ 391 static const float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 392 #else 393 static float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 394 #endif 395 1.5044444979e-07, /* 0x342189db */ 396 7.3223426938e-02, /* 0x3d95f62a */ 397 1.9981917143e+00, /* 0x3fffc4bf */ 398 1.4495602608e+01, /* 0x4167edfd */ 399 3.1666231155e+01, /* 0x41fd5471 */ 400 1.6252708435e+01, /* 0x4182058c */ 401 }; 402 #ifdef __STDC__ 403 static const float qS2[6] = { 404 #else 405 static float qS2[6] = { 406 #endif 407 3.0365585327e+01, /* 0x41f2ecb8 */ 408 2.6934811401e+02, /* 0x4386ac8f */ 409 8.4478375244e+02, /* 0x44533229 */ 410 8.8293585205e+02, /* 0x445cbbe5 */ 411 2.1266638184e+02, /* 0x4354aa98 */ 412 -5.3109550476e+00, /* 0xc0a9f358 */ 413 }; 414 415 #ifdef __STDC__ qzerof(float x)416 static float qzerof(float x) 417 #else 418 static float qzerof(x) 419 float x; 420 #endif 421 { 422 #ifdef __STDC__ 423 const float *p,*q; 424 #else 425 float *p,*q; 426 #endif 427 float s,r,z; 428 __int32_t ix; 429 GET_FLOAT_WORD(ix,x); 430 ix &= 0x7fffffff; 431 if(ix>=0x41000000) {p = qR8; q= qS8;} 432 else if(ix>=0x40f71c58){p = qR5; q= qS5;} 433 else if(ix>=0x4036db68){p = qR3; q= qS3;} 434 else {p = qR2; q= qS2;} 435 z = one/(x*x); 436 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 437 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); 438 return (-(float).125 + r/s)/x; 439 } 440