1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 /* 22 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23 */ 24 /* 25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved. 26 * Use is subject to license terms. 27 */ 28 29 #pragma weak __j0f = j0f 30 #pragma weak __j1f = j1f 31 #pragma weak __jnf = jnf 32 #pragma weak __y0f = y0f 33 #pragma weak __y1f = y1f 34 #pragma weak __ynf = ynf 35 36 #include "libm.h" 37 #include <float.h> 38 39 #if defined(__i386) && !defined(__amd64) 40 extern int __swapRP(int); 41 #endif 42 43 static const float 44 zerof = 0.0f, 45 onef = 1.0f; 46 47 static const double C[] = { 48 0.0, 49 -0.125, 50 0.25, 51 0.375, 52 0.5, 53 1.0, 54 2.0, 55 8.0, 56 0.5641895835477562869480794515607725858441, /* 1/sqrt(pi) */ 57 0.636619772367581343075535053490057448, /* 2/pi */ 58 1.0e9, 59 }; 60 61 #define zero C[0] 62 #define neighth C[1] 63 #define quarter C[2] 64 #define three8 C[3] 65 #define half C[4] 66 #define one C[5] 67 #define two C[6] 68 #define eight C[7] 69 #define isqrtpi C[8] 70 #define tpi C[9] 71 #define big C[10] 72 73 static const double Cj0y0[] = { 74 0.4861344183386052721391238447e5, /* pr */ 75 0.1377662549407112278133438945e6, 76 0.1222466364088289731869114004e6, 77 0.4107070084315176135583353374e5, 78 0.5026073801860637125889039915e4, 79 0.1783193659125479654541542419e3, 80 0.88010344055383421691677564e0, 81 0.4861344183386052721414037058e5, /* ps */ 82 0.1378196632630384670477582699e6, 83 0.1223967185341006542748936787e6, 84 0.4120150243795353639995862617e5, 85 0.5068271181053546392490184353e4, 86 0.1829817905472769960535671664e3, 87 1.0, 88 -0.1731210995701068539185611951e3, /* qr */ 89 -0.5522559165936166961235240613e3, 90 -0.5604935606637346590614529613e3, 91 -0.2200430300226009379477365011e3, 92 -0.323869355375648849771296746e2, 93 -0.14294979207907956223499258e1, 94 -0.834690374102384988158918e-2, 95 0.1107975037248683865326709645e5, /* qs */ 96 0.3544581680627082674651471873e5, 97 0.3619118937918394132179019059e5, 98 0.1439895563565398007471485822e5, 99 0.2190277023344363955930226234e4, 100 0.106695157020407986137501682e3, 101 1.0, 102 }; 103 104 #define pr Cj0y0 105 #define ps (Cj0y0+7) 106 #define qr (Cj0y0+14) 107 #define qs (Cj0y0+21) 108 109 static const double Cj0[] = { 110 -2.500000000000003622131880894830476755537e-0001, /* r0 */ 111 1.095597547334830263234433855932375353303e-0002, 112 -1.819734750463320921799187258987098087697e-0004, 113 9.977001946806131657544212501069893930846e-0007, 114 1.0, /* s0 */ 115 1.867609810662950169966782360588199673741e-0002, 116 1.590389206181565490878430827706972074208e-0004, 117 6.520867386742583632375520147714499522721e-0007, 118 9.999999999999999942156495584397047660949e-0001, /* r1 */ 119 -2.389887722731319130476839836908143731281e-0001, 120 1.293359476138939027791270393439493640570e-0002, 121 -2.770985642343140122168852400228563364082e-0004, 122 2.905241575772067678086738389169625218912e-0006, 123 -1.636846356264052597969042009265043251279e-0008, 124 5.072306160724884775085431059052611737827e-0011, 125 -8.187060730684066824228914775146536139112e-0014, 126 5.422219326959949863954297860723723423842e-0017, 127 1.0, /* s1 */ 128 1.101122772686807702762104741932076228349e-0002, 129 6.140169310641649223411427764669143978228e-0005, 130 2.292035877515152097976946119293215705250e-0007, 131 6.356910426504644334558832036362219583789e-0010, 132 1.366626326900219555045096999553948891401e-0012, 133 2.280399586866739522891837985560481180088e-0015, 134 2.801559820648939665270492520004836611187e-0018, 135 2.073101088320349159764410261466350732968e-0021, 136 }; 137 138 #define r0 Cj0 139 #define s0 (Cj0+4) 140 #define r1 (Cj0+8) 141 #define s1 (Cj0+17) 142 143 static const double Cy0[] = { 144 -7.380429510868722526754723020704317641941e-0002, /* u0 */ 145 1.772607102684869924301459663049874294814e-0001, 146 -1.524370666542713828604078090970799356306e-0002, 147 4.650819100693891757143771557629924591915e-0004, 148 -7.125768872339528975036316108718239946022e-0006, 149 6.411017001656104598327565004771515257146e-0008, 150 -3.694275157433032553021246812379258781665e-0010, 151 1.434364544206266624252820889648445263842e-0012, 152 -3.852064731859936455895036286874139896861e-0015, 153 7.182052899726138381739945881914874579696e-0018, 154 -9.060556574619677567323741194079797987200e-0021, 155 7.124435467408860515265552217131230511455e-0024, 156 -2.709726774636397615328813121715432044771e-0027, 157 1.0, /* v0 */ 158 4.678678931512549002587702477349214886475e-0003, 159 9.486828955529948534822800829497565178985e-0006, 160 1.001495929158861646659010844136682454906e-0008, 161 4.725338116256021660204443235685358593611e-0012, 162 }; 163 164 #define u0 Cy0 165 #define v0 (Cy0+13) 166 167 static const double Cj1y1[] = { 168 -0.4435757816794127857114720794e7, /* pr0 */ 169 -0.9942246505077641195658377899e7, 170 -0.6603373248364939109255245434e7, 171 -0.1523529351181137383255105722e7, 172 -0.1098240554345934672737413139e6, 173 -0.1611616644324610116477412898e4, 174 -0.4435757816794127856828016962e7, /* ps0 */ 175 -0.9934124389934585658967556309e7, 176 -0.6585339479723087072826915069e7, 177 -0.1511809506634160881644546358e7, 178 -0.1072638599110382011903063867e6, 179 -0.1455009440190496182453565068e4, 180 0.3322091340985722351859704442e5, /* qr0 */ 181 0.8514516067533570196555001171e5, 182 0.6617883658127083517939992166e5, 183 0.1849426287322386679652009819e5, 184 0.1706375429020768002061283546e4, 185 0.3526513384663603218592175580e2, 186 0.7087128194102874357377502472e6, /* qs0 */ 187 0.1819458042243997298924553839e7, 188 0.1419460669603720892855755253e7, 189 0.4002944358226697511708610813e6, 190 0.3789022974577220264142952256e5, 191 0.8638367769604990967475517183e3, 192 }; 193 194 #define pr0 Cj1y1 195 #define ps0 (Cj1y1+6) 196 #define qr0 (Cj1y1+12) 197 #define qs0 (Cj1y1+18) 198 199 static const double Cj1[] = { 200 -6.250000000000002203053200981413218949548e-0002, /* a0 */ 201 1.600998455640072901321605101981501263762e-0003, 202 -1.963888815948313758552511884390162864930e-0005, 203 8.263917341093549759781339713418201620998e-0008, 204 1.0e0, /* b0 */ 205 1.605069137643004242395356851797873766927e-0002, 206 1.149454623251299996428500249509098499383e-0004, 207 3.849701673735260970379681807910852327825e-0007, 208 4.999999999999999995517408894340485471724e-0001, 209 -6.003825028120475684835384519945468075423e-0002, 210 2.301719899263321828388344461995355419832e-0003, 211 -4.208494869238892934859525221654040304068e-0005, 212 4.377745135188837783031540029700282443388e-0007, 213 -2.854106755678624335145364226735677754179e-0009, 214 1.234002865443952024332943901323798413689e-0011, 215 -3.645498437039791058951273508838177134310e-0014, 216 7.404320596071797459925377103787837414422e-0017, 217 -1.009457448277522275262808398517024439084e-0019, 218 8.520158355824819796968771418801019930585e-0023, 219 -3.458159926081163274483854614601091361424e-0026, 220 1.0e0, /* b1 */ 221 4.923499437590484879081138588998986303306e-0003, 222 1.054389489212184156499666953501976688452e-0005, 223 1.180768373106166527048240364872043816050e-0008, 224 5.942665743476099355323245707680648588540e-0012, 225 }; 226 227 #define a0 Cj1 228 #define b0 (Cj1+4) 229 #define a1 (Cj1+8) 230 #define b1 (Cj1+20) 231 232 static const double Cy1[] = { 233 -1.960570906462389461018983259589655961560e-0001, /* c0 */ 234 4.931824118350661953459180060007970291139e-0002, 235 -1.626975871565393656845930125424683008677e-0003, 236 1.359657517926394132692884168082224258360e-0005, 237 1.0e0, /* d0 */ 238 2.565807214838390835108224713630901653793e-0002, 239 3.374175208978404268650522752520906231508e-0004, 240 2.840368571306070719539936935220728843177e-0006, 241 1.396387402048998277638900944415752207592e-0008, 242 -1.960570906462389473336339614647555351626e-0001, /* c1 */ 243 5.336268030335074494231369159933012844735e-0002, 244 -2.684137504382748094149184541866332033280e-0003, 245 5.737671618979185736981543498580051903060e-0005, 246 -6.642696350686335339171171785557663224892e-0007, 247 4.692417922568160354012347591960362101664e-0009, 248 -2.161728635907789319335231338621412258355e-0011, 249 6.727353419738316107197644431844194668702e-0014, 250 -1.427502986803861372125234355906790573422e-0016, 251 2.020392498726806769468143219616642940371e-0019, 252 -1.761371948595104156753045457888272716340e-0022, 253 7.352828391941157905175042420249225115816e-0026, 254 1.0e0, /* d1 */ 255 5.029187436727947764916247076102283399442e-0003, 256 1.102693095808242775074856548927801750627e-0005, 257 1.268035774543174837829534603830227216291e-0008, 258 6.579416271766610825192542295821308730206e-0012, 259 }; 260 261 #define c0 Cy1 262 #define d0 (Cy1+4) 263 #define c1 (Cy1+9) 264 #define d1 (Cy1+21) 265 266 267 /* core of j0f computation; assumes fx is finite */ 268 static double 269 __k_j0f(float fx) 270 { 271 double x, z, s, c, ss, cc, r, t, p0, q0; 272 int ix, i; 273 274 ix = *(int *)&fx & ~0x80000000; 275 x = fabs((double)fx); 276 if (ix > 0x41000000) { 277 /* x > 8; see comments in j0.c */ 278 s = sin(x); 279 c = cos(x); 280 if (signbit(s) != signbit(c)) { 281 ss = s - c; 282 cc = -cos(x + x) / ss; 283 } else { 284 cc = s + c; 285 ss = -cos(x + x) / cc; 286 } 287 if (ix > 0x501502f9) { 288 /* x > 1.0e10 */ 289 p0 = one; 290 q0 = neighth / x; 291 } else { 292 t = eight / x; 293 z = t * t; 294 p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] + 295 z * (pr[4] + z * (pr[5] + z * pr[6])))))) / 296 (ps[0] + z * (ps[1] + z * (ps[2] + z * (ps[3] + 297 z * (ps[4] + z * (ps[5] + z)))))); 298 q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] + 299 z * (qr[4] + z * (qr[5] + z * qr[6])))))) / 300 (qs[0] + z * (qs[1] + z * (qs[2] + z * (qs[3] + 301 z * (qs[4] + z * (qs[5] + z))))))) * t; 302 } 303 return (isqrtpi * (p0 * cc - q0 * ss) / sqrt(x)); 304 } 305 if (ix <= 0x3727c5ac) { 306 /* x <= 1.0e-5 */ 307 if (ix <= 0x219392ef) /* x <= 1.0e-18 */ 308 return (one - x); 309 return (one - x * x * quarter); 310 } 311 z = x * x; 312 if (ix <= 0x3fa3d70a) { 313 /* x <= 1.28 */ 314 r = r0[0] + z * (r0[1] + z * (r0[2] + z * r0[3])); 315 s = s0[0] + z * (s0[1] + z * (s0[2] + z * s0[3])); 316 return (one + z * (r / s)); 317 } 318 r = r1[8]; 319 s = s1[8]; 320 for (i = 7; i >= 0; i--) { 321 r = r * z + r1[i]; 322 s = s * z + s1[i]; 323 } 324 return (r / s); 325 } 326 327 float 328 j0f(float fx) 329 { 330 float f; 331 int ix; 332 #if defined(__i386) && !defined(__amd64) 333 int rp; 334 #endif 335 336 ix = *(int *)&fx & ~0x80000000; 337 if (ix >= 0x7f800000) { /* nan or inf */ 338 if (ix > 0x7f800000) 339 return (fx * fx); 340 return (zerof); 341 } 342 343 #if defined(__i386) && !defined(__amd64) 344 rp = __swapRP(fp_extended); 345 #endif 346 f = (float)__k_j0f(fx); 347 #if defined(__i386) && !defined(__amd64) 348 if (rp != fp_extended) 349 (void) __swapRP(rp); 350 #endif 351 return (f); 352 } 353 354 /* core of y0f computation; assumes fx is finite and positive */ 355 static double 356 __k_y0f(float fx) 357 { 358 double x, z, s, c, ss, cc, t, p0, q0, u, v; 359 int ix, i; 360 361 ix = *(int *)&fx; 362 x = (double)fx; 363 if (ix > 0x41000000) { 364 /* x > 8; see comments in j0.c */ 365 s = sin(x); 366 c = cos(x); 367 if (signbit(s) != signbit(c)) { 368 ss = s - c; 369 cc = -cos(x + x) / ss; 370 } else { 371 cc = s + c; 372 ss = -cos(x + x) / cc; 373 } 374 if (ix > 0x501502f9) { 375 /* x > 1.0e10 */ 376 p0 = one; 377 q0 = neighth / x; 378 } else { 379 t = eight / x; 380 z = t * t; 381 p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] + 382 z * (pr[4] + z * (pr[5] + z * pr[6])))))) / 383 (ps[0] + z * (ps[1] + z * (ps[2] + z * (ps[3] + 384 z * (ps[4] + z * (ps[5] + z)))))); 385 q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] + 386 z * (qr[4] + z * (qr[5] + z * qr[6])))))) / 387 (qs[0] + z * (qs[1] + z * (qs[2] + z * (qs[3] + 388 z * (qs[4] + z * (qs[5] + z))))))) * t; 389 } 390 return (isqrtpi * (p0 * ss + q0 * cc) / sqrt(x)); 391 } 392 if (ix <= 0x219392ef) /* x <= 1.0e-18 */ 393 return (u0[0] + tpi * log(x)); 394 z = x * x; 395 u = u0[12]; 396 for (i = 11; i >= 0; i--) 397 u = u * z + u0[i]; 398 v = v0[0] + z * (v0[1] + z * (v0[2] + z * (v0[3] + z * v0[4]))); 399 return (u / v + tpi * (__k_j0f(fx) * log(x))); 400 } 401 402 float 403 y0f(float fx) 404 { 405 float f; 406 int ix; 407 #if defined(__i386) && !defined(__amd64) 408 int rp; 409 #endif 410 411 ix = *(int *)&fx; 412 if ((ix & ~0x80000000) > 0x7f800000) /* nan */ 413 return (fx * fx); 414 if (ix <= 0) { /* zero or negative */ 415 if ((ix << 1) == 0) 416 return (-onef / zerof); 417 return (zerof / zerof); 418 } 419 if (ix == 0x7f800000) /* +inf */ 420 return (zerof); 421 422 #if defined(__i386) && !defined(__amd64) 423 rp = __swapRP(fp_extended); 424 #endif 425 f = (float)__k_y0f(fx); 426 #if defined(__i386) && !defined(__amd64) 427 if (rp != fp_extended) 428 (void) __swapRP(rp); 429 #endif 430 return (f); 431 } 432 433 /* core of j1f computation; assumes fx is finite */ 434 static double 435 __k_j1f(float fx) 436 { 437 double x, z, s, c, ss, cc, r, t, p1, q1; 438 int i, ix, sgn; 439 440 ix = *(int *)&fx; 441 sgn = (unsigned)ix >> 31; 442 ix &= ~0x80000000; 443 x = fabs((double)fx); 444 if (ix > 0x41000000) { 445 /* x > 8; see comments in j1.c */ 446 s = sin(x); 447 c = cos(x); 448 if (signbit(s) != signbit(c)) { 449 cc = s - c; 450 ss = cos(x + x) / cc; 451 } else { 452 ss = -s - c; 453 cc = cos(x + x) / ss; 454 } 455 if (ix > 0x501502f9) { 456 /* x > 1.0e10 */ 457 p1 = one; 458 q1 = three8 / x; 459 } else { 460 t = eight / x; 461 z = t * t; 462 p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z * 463 (pr0[3] + z * (pr0[4] + z * pr0[5]))))) / 464 (ps0[0] + z * (ps0[1] + z * (ps0[2] + z * 465 (ps0[3] + z * (ps0[4] + z * (ps0[5] + z)))))); 466 q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z * 467 (qr0[3] + z * (qr0[4] + z * qr0[5]))))) / 468 (qs0[0] + z * (qs0[1] + z * (qs0[2] + z * 469 (qs0[3] + z * (qs0[4] + z * (qs0[5] + z))))))) * t; 470 } 471 t = isqrtpi * (p1 * cc - q1 * ss) / sqrt(x); 472 return ((sgn)? -t : t); 473 } 474 if (ix <= 0x3727c5ac) { 475 /* x <= 1.0e-5 */ 476 if (ix <= 0x219392ef) /* x <= 1.0e-18 */ 477 t = half * x; 478 else 479 t = x * (half + neighth * x * x); 480 return ((sgn)? -t : t); 481 } 482 z = x * x; 483 if (ix < 0x3fa3d70a) { 484 /* x < 1.28 */ 485 r = a0[0] + z * (a0[1] + z * (a0[2] + z * a0[3])); 486 s = b0[0] + z * (b0[1] + z * (b0[2] + z * b0[3])); 487 t = x * half + x * (z * (r / s)); 488 } else { 489 r = a1[11]; 490 for (i = 10; i >= 0; i--) 491 r = r * z + a1[i]; 492 s = b1[0] + z * (b1[1] + z * (b1[2] + z * (b1[3] + z * b1[4]))); 493 t = x * (r / s); 494 } 495 return ((sgn)? -t : t); 496 } 497 498 float 499 j1f(float fx) 500 { 501 float f; 502 int ix; 503 #if defined(__i386) && !defined(__amd64) 504 int rp; 505 #endif 506 507 ix = *(int *)&fx & ~0x80000000; 508 if (ix >= 0x7f800000) /* nan or inf */ 509 return (onef / fx); 510 511 #if defined(__i386) && !defined(__amd64) 512 rp = __swapRP(fp_extended); 513 #endif 514 f = (float)__k_j1f(fx); 515 #if defined(__i386) && !defined(__amd64) 516 if (rp != fp_extended) 517 (void) __swapRP(rp); 518 #endif 519 return (f); 520 } 521 522 /* core of y1f computation; assumes fx is finite and positive */ 523 static double 524 __k_y1f(float fx) 525 { 526 double x, z, s, c, ss, cc, u, v, p1, q1, t; 527 int i, ix; 528 529 ix = *(int *)&fx; 530 x = (double)fx; 531 if (ix > 0x41000000) { 532 /* x > 8; see comments in j1.c */ 533 s = sin(x); 534 c = cos(x); 535 if (signbit(s) != signbit(c)) { 536 cc = s - c; 537 ss = cos(x + x) / cc; 538 } else { 539 ss = -s - c; 540 cc = cos(x + x) / ss; 541 } 542 if (ix > 0x501502f9) { 543 /* x > 1.0e10 */ 544 p1 = one; 545 q1 = three8 / x; 546 } else { 547 t = eight / x; 548 z = t * t; 549 p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z * 550 (pr0[3] + z * (pr0[4] + z * pr0[5]))))) / 551 (ps0[0] + z * (ps0[1] + z * (ps0[2] + z * 552 (ps0[3] + z * (ps0[4] + z * (ps0[5] + z)))))); 553 q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z * 554 (qr0[3] + z * (qr0[4] + z * qr0[5]))))) / 555 (qs0[0] + z * (qs0[1] + z * (qs0[2] + z * 556 (qs0[3] + z * (qs0[4] + z * (qs0[5] + z))))))) * t; 557 } 558 return (isqrtpi * (p1 * ss + q1 * cc) / sqrt(x)); 559 } 560 if (ix <= 0x219392ef) /* x <= 1.0e-18 */ 561 return (-tpi / x); 562 z = x * x; 563 if (ix < 0x3fa3d70a) { 564 /* x < 1.28 */ 565 u = c0[0] + z * (c0[1] + z * (c0[2] + z * c0[3])); 566 v = d0[0] + z * (d0[1] + z * (d0[2] + z * (d0[3] + z * d0[4]))); 567 } else { 568 u = c1[11]; 569 for (i = 10; i >= 0; i--) 570 u = u * z + c1[i]; 571 v = d1[0] + z * (d1[1] + z * (d1[2] + z * (d1[3] + z * d1[4]))); 572 } 573 return (x * (u / v) + tpi * (__k_j1f(fx) * log(x) - one / x)); 574 } 575 576 float 577 y1f(float fx) 578 { 579 float f; 580 int ix; 581 #if defined(__i386) && !defined(__amd64) 582 int rp; 583 #endif 584 585 ix = *(int *)&fx; 586 if ((ix & ~0x80000000) > 0x7f800000) /* nan */ 587 return (fx * fx); 588 if (ix <= 0) { /* zero or negative */ 589 if ((ix << 1) == 0) 590 return (-onef / zerof); 591 return (zerof / zerof); 592 } 593 if (ix == 0x7f800000) /* +inf */ 594 return (zerof); 595 596 #if defined(__i386) && !defined(__amd64) 597 rp = __swapRP(fp_extended); 598 #endif 599 f = (float)__k_y1f(fx); 600 #if defined(__i386) && !defined(__amd64) 601 if (rp != fp_extended) 602 (void) __swapRP(rp); 603 #endif 604 return (f); 605 } 606 607 float 608 jnf(int n, float fx) 609 { 610 double a, b, temp, x, z, w, t, q0, q1, h; 611 float f; 612 int i, ix, sgn, m, k; 613 #if defined(__i386) && !defined(__amd64) 614 int rp; 615 #endif 616 617 if (n < 0) { 618 n = -n; 619 fx = -fx; 620 } 621 if (n == 0) 622 return (j0f(fx)); 623 if (n == 1) 624 return (j1f(fx)); 625 626 ix = *(int *)&fx; 627 sgn = (n & 1)? ((unsigned)ix >> 31) : 0; 628 ix &= ~0x80000000; 629 if (ix >= 0x7f800000) { /* nan or inf */ 630 if (ix > 0x7f800000) 631 return (fx * fx); 632 return ((sgn)? -zerof : zerof); 633 } 634 if ((ix << 1) == 0) 635 return ((sgn)? -zerof : zerof); 636 637 #if defined(__i386) && !defined(__amd64) 638 rp = __swapRP(fp_extended); 639 #endif 640 fx = fabsf(fx); 641 x = (double)fx; 642 if ((double)n <= x) { 643 /* safe to use J(n+1,x) = 2n/x * J(n,x) - J(n-1,x) */ 644 a = __k_j0f(fx); 645 b = __k_j1f(fx); 646 for (i = 1; i < n; i++) { 647 temp = b; 648 b = b * ((double)(i + i) / x) - a; 649 a = temp; 650 } 651 f = (float)b; 652 #if defined(__i386) && !defined(__amd64) 653 if (rp != fp_extended) 654 (void) __swapRP(rp); 655 #endif 656 return ((sgn)? -f : f); 657 } 658 if (ix < 0x3089705f) { 659 /* x < 1.0e-9; use J(n,x) = 1/n! * (x / 2)^n */ 660 if (n > 6) 661 n = 6; /* result underflows to zero for n >= 6 */ 662 b = t = half * x; 663 a = one; 664 for (i = 2; i <= n; i++) { 665 b *= t; 666 a *= (double)i; 667 } 668 b /= a; 669 } else { 670 /* 671 * Use the backward recurrence: 672 * 673 * x x^2 x^2 674 * J(n,x)/J(n-1,x) = ---- - ------ - ------ ..... 675 * 2n 2(n+1) 2(n+2) 676 * 677 * Let w = 2n/x and h = 2/x. Then the above quotient 678 * is equal to the continued fraction: 679 * 1 680 * = ----------------------- 681 * 1 682 * w - ----------------- 683 * 1 684 * w+h - --------- 685 * w+2h - ... 686 * 687 * To determine how many terms are needed, run the 688 * recurrence 689 * 690 * Q(0) = w, 691 * Q(1) = w(w+h) - 1, 692 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2). 693 * 694 * Then when Q(k) > 1e4, k is large enough for single 695 * precision. 696 */ 697 /* XXX NOT DONE - rework this */ 698 w = (n + n) / x; 699 h = two / x; 700 q0 = w; 701 z = w + h; 702 q1 = w * z - one; 703 k = 1; 704 while (q1 < big) { 705 k++; 706 z += h; 707 temp = z * q1 - q0; 708 q0 = q1; 709 q1 = temp; 710 } 711 m = n + n; 712 t = zero; 713 for (i = (n + k) << 1; i >= m; i -= 2) 714 t = one / ((double)i / x - t); 715 a = t; 716 b = one; 717 /* 718 * estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) 719 * hence, if n*(log(2n/x)) > ... 720 * single 8.8722839355e+01 721 * double 7.09782712893383973096e+02 722 * then recurrent value may overflow and the result is 723 * likely underflow to zero 724 */ 725 temp = (double)n; 726 temp *= log((two / x) * temp); 727 if (temp < 7.09782712893383973096e+02) { 728 for (i = n - 1; i > 0; i--) { 729 temp = b; 730 b = b * ((double)(i + i) / x) - a; 731 a = temp; 732 } 733 } else { 734 for (i = n - 1; i > 0; i--) { 735 temp = b; 736 b = b * ((double)(i + i) / x) - a; 737 a = temp; 738 if (b > 1.0e100) { 739 a /= b; 740 t /= b; 741 b = one; 742 } 743 } 744 } 745 b = (t * __k_j0f(fx) / b); 746 } 747 f = (float)b; 748 #if defined(__i386) && !defined(__amd64) 749 if (rp != fp_extended) 750 (void) __swapRP(rp); 751 #endif 752 return ((sgn)? -f : f); 753 } 754 755 float 756 ynf(int n, float fx) 757 { 758 double a, b, temp, x; 759 float f; 760 int i, sign, ix; 761 #if defined(__i386) && !defined(__amd64) 762 int rp; 763 #endif 764 765 sign = 0; 766 if (n < 0) { 767 n = -n; 768 if (n & 1) 769 sign = 1; 770 } 771 if (n == 0) 772 return (y0f(fx)); 773 if (n == 1) 774 return ((sign)? -y1f(fx) : y1f(fx)); 775 776 ix = *(int *)&fx; 777 if ((ix & ~0x80000000) > 0x7f800000) /* nan */ 778 return (fx * fx); 779 if (ix <= 0) { /* zero or negative */ 780 if ((ix << 1) == 0) 781 return (-onef / zerof); 782 return (zerof / zerof); 783 } 784 if (ix == 0x7f800000) /* +inf */ 785 return (zerof); 786 787 #if defined(__i386) && !defined(__amd64) 788 rp = __swapRP(fp_extended); 789 #endif 790 a = __k_y0f(fx); 791 b = __k_y1f(fx); 792 x = (double)fx; 793 for (i = 1; i < n; i++) { 794 temp = b; 795 b *= (double)(i + i) / x; 796 if (b <= -DBL_MAX) 797 break; 798 b -= a; 799 a = temp; 800 } 801 f = (float)b; 802 #if defined(__i386) && !defined(__amd64) 803 if (rp != fp_extended) 804 (void) __swapRP(rp); 805 #endif 806 return ((sign)? -f : f); 807 } 808