Print this page
11210 libm should be cstyle(1ONBLD) clean
   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]


  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 }
   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 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 
  26 /*
  27  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __j0f = j0f
  32 #pragma weak __j1f = j1f
  33 #pragma weak __jnf = jnf
  34 #pragma weak __y0f = y0f
  35 #pragma weak __y1f = y1f
  36 #pragma weak __ynf = ynf
  37 
  38 #include "libm.h"
  39 #include <float.h>
  40 
  41 #if defined(__i386) && !defined(__amd64)
  42 extern int __swapRP(int);
  43 #endif
  44 
  45 static const float zerof = 0.0f, 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]


  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 /* core of j0f computation; assumes fx is finite */
 267 static double
 268 __k_j0f(float fx)
 269 {
 270         double x, z, s, c, ss, cc, r, t, p0, q0;
 271         int ix, i;
 272 
 273         ix = *(int *)&fx & ~0x80000000;
 274         x = fabs((double)fx);
 275 
 276         if (ix > 0x41000000) {
 277                 /* x > 8; see comments in j0.c */
 278                 s = sin(x);
 279                 c = cos(x);
 280 
 281                 if (signbit(s) != signbit(c)) {
 282                         ss = s - c;
 283                         cc = -cos(x + x) / ss;
 284                 } else {
 285                         cc = s + c;
 286                         ss = -cos(x + x) / cc;
 287                 }
 288 
 289                 if (ix > 0x501502f9) {
 290                         /* x > 1.0e10 */
 291                         p0 = one;
 292                         q0 = neighth / x;
 293                 } else {
 294                         t = eight / x;
 295                         z = t * t;
 296                         p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] + z *
 297                             (pr[4] + z * (pr[5] + z * pr[6])))))) / (ps[0] + z *
 298                             (ps[1] + z * (ps[2] + z * (ps[3] + z * (ps[4] + z *
 299                             (ps[5] + z))))));
 300                         q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
 301                             z * (qr[4] + z * (qr[5] + z * qr[6])))))) / (qs[0] +
 302                             z * (qs[1] + z * (qs[2] + z * (qs[3] + z * (qs[4] +
 303                             z * (qs[5] + z))))))) * t;
 304                 }
 305 
 306                 return (isqrtpi * (p0 * cc - q0 * ss) / sqrt(x));
 307         }
 308 
 309         if (ix <= 0x3727c5ac) {
 310                 /* x <= 1.0e-5 */
 311                 if (ix <= 0x219392ef)        /* x <= 1.0e-18 */
 312                         return (one - x);
 313 
 314                 return (one - x * x * quarter);
 315         }
 316 
 317         z = x * x;
 318 
 319         if (ix <= 0x3fa3d70a) {
 320                 /* x <= 1.28 */
 321                 r = r0[0] + z * (r0[1] + z * (r0[2] + z * r0[3]));
 322                 s = s0[0] + z * (s0[1] + z * (s0[2] + z * s0[3]));
 323                 return (one + z * (r / s));
 324         }
 325 
 326         r = r1[8];
 327         s = s1[8];
 328 
 329         for (i = 7; i >= 0; i--) {
 330                 r = r * z + r1[i];
 331                 s = s * z + s1[i];
 332         }
 333 
 334         return (r / s);
 335 }
 336 
 337 float
 338 j0f(float fx)
 339 {
 340         float f;
 341         int ix;
 342 
 343 #if defined(__i386) && !defined(__amd64)
 344         int rp;
 345 #endif
 346 
 347         ix = *(int *)&fx & ~0x80000000;
 348 
 349         if (ix >= 0x7f800000) {              /* nan or inf */
 350                 if (ix > 0x7f800000)
 351                         return (fx * fx);
 352 
 353                 return (zerof);
 354         }
 355 
 356 #if defined(__i386) && !defined(__amd64)
 357         rp = __swapRP(fp_extended);
 358 #endif
 359         f = (float)__k_j0f(fx);
 360 #if defined(__i386) && !defined(__amd64)
 361         if (rp != fp_extended)
 362                 (void) __swapRP(rp);
 363 #endif
 364         return (f);
 365 }
 366 
 367 /* core of y0f computation; assumes fx is finite and positive */
 368 static double
 369 __k_y0f(float fx)
 370 {
 371         double x, z, s, c, ss, cc, t, p0, q0, u, v;
 372         int ix, i;
 373 
 374         ix = *(int *)&fx;
 375         x = (double)fx;
 376 
 377         if (ix > 0x41000000) {
 378                 /* x > 8; see comments in j0.c */
 379                 s = sin(x);
 380                 c = cos(x);
 381 
 382                 if (signbit(s) != signbit(c)) {
 383                         ss = s - c;
 384                         cc = -cos(x + x) / ss;
 385                 } else {
 386                         cc = s + c;
 387                         ss = -cos(x + x) / cc;
 388                 }
 389 
 390                 if (ix > 0x501502f9) {
 391                         /* x > 1.0e10 */
 392                         p0 = one;
 393                         q0 = neighth / x;
 394                 } else {
 395                         t = eight / x;
 396                         z = t * t;
 397                         p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] + z *
 398                             (pr[4] + z * (pr[5] + z * pr[6])))))) / (ps[0] + z *
 399                             (ps[1] + z * (ps[2] + z * (ps[3] + z * (ps[4] + z *
 400                             (ps[5] + z))))));
 401                         q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
 402                             z * (qr[4] + z * (qr[5] + z * qr[6])))))) / (qs[0] +
 403                             z * (qs[1] + z * (qs[2] + z * (qs[3] + z * (qs[4] +
 404                             z * (qs[5] + z))))))) * t;
 405                 }
 406 
 407                 return (isqrtpi * (p0 * ss + q0 * cc) / sqrt(x));
 408         }
 409 
 410         if (ix <= 0x219392ef)                /* x <= 1.0e-18 */
 411                 return (u0[0] + tpi * log(x));
 412 
 413         z = x * x;
 414         u = u0[12];
 415 
 416         for (i = 11; i >= 0; i--)
 417                 u = u * z + u0[i];
 418 
 419         v = v0[0] + z * (v0[1] + z * (v0[2] + z * (v0[3] + z * v0[4])));
 420         return (u / v + tpi * (__k_j0f(fx) * log(x)));
 421 }
 422 
 423 float
 424 y0f(float fx)
 425 {
 426         float f;
 427         int ix;
 428 
 429 #if defined(__i386) && !defined(__amd64)
 430         int rp;
 431 #endif
 432 
 433         ix = *(int *)&fx;
 434 
 435         if ((ix & ~0x80000000) > 0x7f800000)     /* nan */
 436                 return (fx * fx);
 437 
 438         if (ix <= 0) {                               /* zero or negative */
 439                 if ((ix << 1) == 0)
 440                         return (-onef / zerof);
 441 
 442                 return (zerof / zerof);
 443         }
 444 
 445         if (ix == 0x7f800000)           /* +inf */
 446                 return (zerof);
 447 
 448 #if defined(__i386) && !defined(__amd64)
 449         rp = __swapRP(fp_extended);
 450 #endif
 451         f = (float)__k_y0f(fx);
 452 #if defined(__i386) && !defined(__amd64)
 453         if (rp != fp_extended)
 454                 (void) __swapRP(rp);
 455 #endif
 456         return (f);
 457 }
 458 
 459 /* core of j1f computation; assumes fx is finite */
 460 static double
 461 __k_j1f(float fx)
 462 {
 463         double x, z, s, c, ss, cc, r, t, p1, q1;
 464         int i, ix, sgn;
 465 
 466         ix = *(int *)&fx;
 467         sgn = (unsigned)ix >> 31;
 468         ix &= ~0x80000000;
 469         x = fabs((double)fx);
 470 
 471         if (ix > 0x41000000) {
 472                 /* x > 8; see comments in j1.c */
 473                 s = sin(x);
 474                 c = cos(x);
 475 
 476                 if (signbit(s) != signbit(c)) {
 477                         cc = s - c;
 478                         ss = cos(x + x) / cc;
 479                 } else {
 480                         ss = -s - c;
 481                         cc = cos(x + x) / ss;
 482                 }
 483 
 484                 if (ix > 0x501502f9) {
 485                         /* x > 1.0e10 */
 486                         p1 = one;
 487                         q1 = three8 / x;
 488                 } else {
 489                         t = eight / x;
 490                         z = t * t;
 491                         p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z * (pr0[3] +
 492                             z * (pr0[4] + z * pr0[5]))))) / (ps0[0] + z *
 493                             (ps0[1] + z * (ps0[2] + z * (ps0[3] + z * (ps0[4] +
 494                             z * (ps0[5] + z))))));
 495                         q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
 496                             (qr0[3] + z * (qr0[4] + z * qr0[5]))))) / (qs0[0] +
 497                             z * (qs0[1] + z * (qs0[2] + z * (qs0[3] + z *
 498                             (qs0[4] + z * (qs0[5] + z))))))) * t;
 499                 }
 500 
 501                 t = isqrtpi * (p1 * cc - q1 * ss) / sqrt(x);
 502                 return ((sgn) ? -t : t);
 503         }
 504 
 505         if (ix <= 0x3727c5ac) {
 506                 /* x <= 1.0e-5 */
 507                 if (ix <= 0x219392ef)        /* x <= 1.0e-18 */
 508                         t = half * x;
 509                 else
 510                         t = x * (half + neighth * x * x);
 511 
 512                 return ((sgn) ? -t : t);
 513         }
 514 
 515         z = x * x;
 516 
 517         if (ix < 0x3fa3d70a) {
 518                 /* x < 1.28 */
 519                 r = a0[0] + z * (a0[1] + z * (a0[2] + z * a0[3]));
 520                 s = b0[0] + z * (b0[1] + z * (b0[2] + z * b0[3]));
 521                 t = x * half + x * (z * (r / s));
 522         } else {
 523                 r = a1[11];
 524 
 525                 for (i = 10; i >= 0; i--)
 526                         r = r * z + a1[i];
 527 
 528                 s = b1[0] + z * (b1[1] + z * (b1[2] + z * (b1[3] + z * b1[4])));
 529                 t = x * (r / s);
 530         }
 531 
 532         return ((sgn) ? -t : t);
 533 }
 534 
 535 float
 536 j1f(float fx)
 537 {
 538         float f;
 539         int ix;
 540 
 541 #if defined(__i386) && !defined(__amd64)
 542         int rp;
 543 #endif
 544 
 545         ix = *(int *)&fx & ~0x80000000;
 546 
 547         if (ix >= 0x7f800000)                /* nan or inf */
 548                 return (onef / fx);
 549 
 550 #if defined(__i386) && !defined(__amd64)
 551         rp = __swapRP(fp_extended);
 552 #endif
 553         f = (float)__k_j1f(fx);
 554 #if defined(__i386) && !defined(__amd64)
 555         if (rp != fp_extended)
 556                 (void) __swapRP(rp);
 557 #endif
 558         return (f);
 559 }
 560 
 561 /* core of y1f computation; assumes fx is finite and positive */
 562 static double
 563 __k_y1f(float fx)
 564 {
 565         double x, z, s, c, ss, cc, u, v, p1, q1, t;
 566         int i, ix;
 567 
 568         ix = *(int *)&fx;
 569         x = (double)fx;
 570 
 571         if (ix > 0x41000000) {
 572                 /* x > 8; see comments in j1.c */
 573                 s = sin(x);
 574                 c = cos(x);
 575 
 576                 if (signbit(s) != signbit(c)) {
 577                         cc = s - c;
 578                         ss = cos(x + x) / cc;
 579                 } else {
 580                         ss = -s - c;
 581                         cc = cos(x + x) / ss;
 582                 }
 583 
 584                 if (ix > 0x501502f9) {
 585                         /* x > 1.0e10 */
 586                         p1 = one;
 587                         q1 = three8 / x;
 588                 } else {
 589                         t = eight / x;
 590                         z = t * t;
 591                         p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z * (pr0[3] +
 592                             z * (pr0[4] + z * pr0[5]))))) / (ps0[0] + z *
 593                             (ps0[1] + z * (ps0[2] + z * (ps0[3] + z * (ps0[4] +
 594                             z * (ps0[5] + z))))));
 595                         q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
 596                             (qr0[3] + z * (qr0[4] + z * qr0[5]))))) / (qs0[0] +
 597                             z * (qs0[1] + z * (qs0[2] + z * (qs0[3] + z *
 598                             (qs0[4] + z * (qs0[5] + z))))))) * t;
 599                 }
 600 
 601                 return (isqrtpi * (p1 * ss + q1 * cc) / sqrt(x));
 602         }
 603 
 604         if (ix <= 0x219392ef)                /* x <= 1.0e-18 */
 605                 return (-tpi / x);
 606 
 607         z = x * x;
 608 
 609         if (ix < 0x3fa3d70a) {
 610                 /* x < 1.28 */
 611                 u = c0[0] + z * (c0[1] + z * (c0[2] + z * c0[3]));
 612                 v = d0[0] + z * (d0[1] + z * (d0[2] + z * (d0[3] + z * d0[4])));
 613         } else {
 614                 u = c1[11];
 615 
 616                 for (i = 10; i >= 0; i--)
 617                         u = u * z + c1[i];
 618 
 619                 v = d1[0] + z * (d1[1] + z * (d1[2] + z * (d1[3] + z * d1[4])));
 620         }
 621 
 622         return (x * (u / v) + tpi * (__k_j1f(fx) * log(x) - one / x));
 623 }
 624 
 625 float
 626 y1f(float fx)
 627 {
 628         float f;
 629         int ix;
 630 
 631 #if defined(__i386) && !defined(__amd64)
 632         int rp;
 633 #endif
 634 
 635         ix = *(int *)&fx;
 636 
 637         if ((ix & ~0x80000000) > 0x7f800000)     /* nan */
 638                 return (fx * fx);
 639 
 640         if (ix <= 0) {                               /* zero or negative */
 641                 if ((ix << 1) == 0)
 642                         return (-onef / zerof);
 643 
 644                 return (zerof / zerof);
 645         }
 646 
 647         if (ix == 0x7f800000)           /* +inf */
 648                 return (zerof);
 649 
 650 #if defined(__i386) && !defined(__amd64)
 651         rp = __swapRP(fp_extended);
 652 #endif
 653         f = (float)__k_y1f(fx);
 654 #if defined(__i386) && !defined(__amd64)
 655         if (rp != fp_extended)
 656                 (void) __swapRP(rp);
 657 #endif
 658         return (f);
 659 }
 660 
 661 float
 662 jnf(int n, float fx)
 663 {
 664         double a, b, temp, x, z, w, t, q0, q1, h;
 665         float f;
 666         int i, ix, sgn, m, k;
 667 
 668 #if defined(__i386) && !defined(__amd64)
 669         int rp;
 670 #endif
 671 
 672         if (n < 0) {
 673                 n = -n;
 674                 fx = -fx;
 675         }
 676 
 677         if (n == 0)
 678                 return (j0f(fx));
 679 
 680         if (n == 1)
 681                 return (j1f(fx));
 682 
 683         ix = *(int *)&fx;
 684         sgn = (n & 1) ? ((unsigned)ix >> 31) : 0;
 685         ix &= ~0x80000000;
 686 
 687         if (ix >= 0x7f800000) {              /* nan or inf */
 688                 if (ix > 0x7f800000)
 689                         return (fx * fx);
 690 
 691                 return ((sgn) ? -zerof : zerof);
 692         }
 693 
 694         if ((ix << 1) == 0)
 695                 return ((sgn) ? -zerof : zerof);
 696 
 697 #if defined(__i386) && !defined(__amd64)
 698         rp = __swapRP(fp_extended);
 699 #endif
 700         fx = fabsf(fx);
 701         x = (double)fx;
 702 
 703         if ((double)n <= x) {
 704                 /* safe to use J(n+1,x) = 2n/x * J(n,x) - J(n-1,x) */
 705                 a = __k_j0f(fx);
 706                 b = __k_j1f(fx);
 707 
 708                 for (i = 1; i < n; i++) {
 709                         temp = b;
 710                         b = b * ((double)(i + i) / x) - a;
 711                         a = temp;
 712                 }
 713 
 714                 f = (float)b;
 715 #if defined(__i386) && !defined(__amd64)
 716                 if (rp != fp_extended)
 717                         (void) __swapRP(rp);
 718 #endif
 719                 return ((sgn) ? -f : f);
 720         }
 721 
 722         if (ix < 0x3089705f) {
 723                 /* x < 1.0e-9; use J(n,x) = 1/n! * (x / 2)^n */
 724                 if (n > 6)
 725                         n = 6;  /* result underflows to zero for n >= 6 */
 726 
 727                 b = t = half * x;
 728                 a = one;
 729 
 730                 for (i = 2; i <= n; i++) {
 731                         b *= t;
 732                         a *= (double)i;
 733                 }
 734 
 735                 b /= a;
 736         } else {
 737                 /* BEGIN CSTYLED */
 738                 /*
 739                  * Use the backward recurrence:
 740                  *
 741                  *                      x      x^2      x^2
 742                  *  J(n,x)/J(n-1,x) =  ---- - ------ - ------   .....
 743                  *                      2n    2(n+1)   2(n+2)
 744                  *
 745                  * Let w = 2n/x and h = 2/x.  Then the above quotient
 746                  * is equal to the continued fraction:
 747                  *                   1
 748                  *      = -----------------------
 749                  *                      1
 750                  *         w - -----------------
 751                  *                        1
 752                  *              w+h - ---------
 753                  *                      w+2h - ...
 754                  *
 755                  * To determine how many terms are needed, run the
 756                  * recurrence
 757                  *
 758                  *      Q(0) = w,
 759                  *      Q(1) = w(w+h) - 1,
 760                  *      Q(k) = (w+k*h)*Q(k-1) - Q(k-2).
 761                  *
 762                  * Then when Q(k) > 1e4, k is large enough for single
 763                  * precision.
 764                  */
 765                 /* END CSTYLED  */
 766 /* XXX NOT DONE - rework this */
 767                 w = (n + n) / x;
 768                 h = two / x;
 769                 q0 = w;
 770                 z = w + h;
 771                 q1 = w * z - one;
 772                 k = 1;
 773 
 774                 while (q1 < big) {
 775                         k++;
 776                         z += h;
 777                         temp = z * q1 - q0;
 778                         q0 = q1;
 779                         q1 = temp;
 780                 }
 781 
 782                 m = n + n;
 783                 t = zero;
 784 
 785                 for (i = (n + k) << 1; i >= m; i -= 2)
 786                         t = one / ((double)i / x - t);
 787 
 788                 a = t;
 789                 b = one;
 790 
 791                 /*
 792                  * estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
 793                  * hence, if n*(log(2n/x)) > ...
 794                  *      single 8.8722839355e+01
 795                  *      double 7.09782712893383973096e+02
 796                  *      then recurrent value may overflow and the result is
 797                  *      likely underflow to zero
 798                  */
 799                 temp = (double)n;
 800                 temp *= log((two / x) * temp);
 801 
 802                 if (temp < 7.09782712893383973096e+02) {
 803                         for (i = n - 1; i > 0; i--) {
 804                                 temp = b;
 805                                 b = b * ((double)(i + i) / x) - a;
 806                                 a = temp;
 807                         }
 808                 } else {
 809                         for (i = n - 1; i > 0; i--) {
 810                                 temp = b;
 811                                 b = b * ((double)(i + i) / x) - a;
 812                                 a = temp;
 813 
 814                                 if (b > 1.0e100) {
 815                                         a /= b;
 816                                         t /= b;
 817                                         b = one;
 818                                 }
 819                         }
 820                 }
 821 
 822                 b = (t * __k_j0f(fx) / b);
 823         }
 824 
 825         f = (float)b;
 826 #if defined(__i386) && !defined(__amd64)
 827         if (rp != fp_extended)
 828                 (void) __swapRP(rp);
 829 #endif
 830         return ((sgn) ? -f : f);
 831 }
 832 
 833 float
 834 ynf(int n, float fx)
 835 {
 836         double a, b, temp, x;
 837         float f;
 838         int i, sign, ix;
 839 
 840 #if defined(__i386) && !defined(__amd64)
 841         int rp;
 842 #endif
 843 
 844         sign = 0;
 845 
 846         if (n < 0) {
 847                 n = -n;
 848 
 849                 if (n & 1)
 850                         sign = 1;
 851         }
 852 
 853         if (n == 0)
 854                 return (y0f(fx));
 855 
 856         if (n == 1)
 857                 return ((sign) ? -y1f(fx) : y1f(fx));
 858 
 859         ix = *(int *)&fx;
 860 
 861         if ((ix & ~0x80000000) > 0x7f800000)     /* nan */
 862                 return (fx * fx);
 863 
 864         if (ix <= 0) {                               /* zero or negative */
 865                 if ((ix << 1) == 0)
 866                         return (-onef / zerof);
 867 
 868                 return (zerof / zerof);
 869         }
 870 
 871         if (ix == 0x7f800000)           /* +inf */
 872                 return (zerof);
 873 
 874 #if defined(__i386) && !defined(__amd64)
 875         rp = __swapRP(fp_extended);
 876 #endif
 877         a = __k_y0f(fx);
 878         b = __k_y1f(fx);
 879         x = (double)fx;
 880 
 881         for (i = 1; i < n; i++) {
 882                 temp = b;
 883                 b *= (double)(i + i) / x;
 884 
 885                 if (b <= -DBL_MAX)
 886                         break;
 887 
 888                 b -= a;
 889                 a = temp;
 890         }
 891 
 892         f = (float)b;
 893 #if defined(__i386) && !defined(__amd64)
 894         if (rp != fp_extended)
 895                 (void) __swapRP(rp);
 896 #endif
 897         return ((sign) ? -f : f);
 898 }