Print this page
11210 libm should be cstyle(1ONBLD) clean

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