usr/src/lib/libm/common/R/besself.c

   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 }