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

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/m9x/tgamma.c
          +++ new/usr/src/lib/libm/common/m9x/tgamma.c
↓ open down ↓ 14 lines elided ↑ open up ↑
  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   23   * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24   24   */
       25 +
  25   26  /*
  26   27   * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27   28   * Use is subject to license terms.
  28   29   */
  29   30  
  30   31  #pragma weak __tgamma = tgamma
  31   32  
  32      -/* INDENT OFF */
       33 +/* BEGIN CSTYLED */
  33   34  /*
  34   35   * True gamma function
  35   36   * double tgamma(double x)
  36   37   *
  37   38   * Error:
  38   39   * ------
  39      - *      Less that one ulp for both positive and negative arguments.
       40 + *      Less that one ulp for both positive and negative arguments.
  40   41   *
  41   42   * Algorithm:
  42   43   * ---------
  43   44   *      A: For negative argument
  44   45   *              (1) gamma(-n or -inf) is NaN
  45   46   *              (2) Underflow Threshold
  46   47   *              (3) Reduction to gamma(1+x)
  47   48   *      B: For x between 1 and 2
  48      - *      C: For x between 0 and 1
       49 + *      C: For x between 0 and 1
  49   50   *      D: For x between 2 and 8
  50   51   *      E: Overflow thresold {see over.c}
  51   52   *      F: For overflow_threshold >= x >= 8
  52   53   *
  53   54   * Implementation details
  54   55   * -----------------------
  55   56   *                                                      -pi
  56   57   * (A) For negative argument, use gamma(-x) = ------------------------.
  57   58   *                                            (sin(pi*x)*gamma(1+x))
  58   59   *
↓ open down ↓ 27 lines elided ↑ open up ↑
  86   87   *      Since x = k+z,
  87   88   *                                                  k+1
  88   89   *              -sin(x*pi) = -sin(k*pi+z*pi) = (-1)   *sin(z*pi),
  89   90   *                               k+1
  90   91   *      we have -kpsin(x) = (-1)   * kpsin(z).  We can further
  91   92   *      reduce z to t by
  92   93   *         (I)   t = z       when 0.00000     <= z < 0.31830...
  93   94   *         (II)  t = 0.5-z   when 0.31830...  <= z < 0.681690...
  94   95   *         (III) t = 1-z     when 0.681690... <= z < 1.00000
  95   96   *      and correspondingly
  96      - *         (I)   kpsin(z) = kpsin(t)    ... 0<= z < 0.3184
  97      - *         (II)  kpsin(z) = kpcos(t)    ... |t|   < 0.182
  98      - *         (III) kpsin(z) = kpsin(t)    ... 0<= t < 0.3184
       97 + *         (I)   kpsin(z) = kpsin(t)    ... 0<= z < 0.3184
       98 + *         (II)  kpsin(z) = kpcos(t)    ... |t|   < 0.182
       99 + *         (III) kpsin(z) = kpsin(t)    ... 0<= t < 0.3184
  99  100   *
 100  101   *      Using a special Remez algorithm, we obtain the following polynomial
 101  102   *      approximation for kpsin(t) for 0<=t<0.3184:
 102  103   *
 103  104   *      Computation note: in simulating higher precision arithmetic, kcpsin
 104  105   *      return head = t and tail = ks[0]*t^3 + (...) to maintain extra bits.
 105  106   *
 106  107   *      Quad precision, remez error <= 2**(-129.74)
 107  108   *                                   3            5                   27
 108  109   *          kpsin(t) = t + ks[0] * t  + ks[1] * t  + ... + ks[12] * t
↓ open down ↓ 28 lines elided ↑ open up ↑
 137  138   *                                  3            5                  9
 138  139   *          kpsin(t) = t + ks[0] * t  + ks[1] * t  + ... + ks[3] * t
 139  140   *
 140  141   *       ks[0] =  -1.64493404985645811354476665052005342839447790544e+0000
 141  142   *       ks[1] =   8.11740794458351064092797249069438269367389272270e-0001
 142  143   *       ks[2] =  -1.90703144603551216933075809162889536878854055202e-0001
 143  144   *       ks[3] =   2.55742333994264563281155312271481108635575331201e-0002
 144  145   *
 145  146   *      Computation note: in simulating higher precision arithmetic, kcpsin
 146  147   *      return head = t and tail = kc[0]*t^3 + (...) to maintain extra bits
 147      - *      precision.
      148 + *      precision.
 148  149   *
 149  150   *      And for kpcos(t) for |t|< 0.183:
 150  151   *
 151  152   *      Quad precision, remez <= 2**(-122.48)
 152  153   *                                     2            4                  22
 153  154   *          kpcos(t) = 1/pi +  pi/2 * t  + kc[2] * t + ... + kc[11] * t
 154  155   *
 155  156   *       kc[2] =   1.29192819501249250731151312779548918765320728489e+0000
 156  157   *       kc[3] =  -4.25027339979557573976029596929319207009444090366e-0001
 157  158   *       kc[4] =   7.49080661650990096109672954618317623888421628613e-0002
↓ open down ↓ 27 lines elided ↑ open up ↑
 185  186   *      Computation note: in simulating higher precision arithmetic, kcpcos
 186  187   *      return head = 1/pi chopped, and tail = pi/2 *t^2 + (tail part of 1/pi
 187  188   *      + ...) to maintain extra bits precision. In particular, pi/2 * t^2
 188  189   *      is calculated with great care.
 189  190   *
 190  191   *      Thus, the computation of gamma(-x), x>0, is:
 191  192   *      Let k = int(x), z = x-k.
 192  193   *      For z in (I)
 193  194   *                                    k+1
 194  195   *                                (-1)
 195      - *              gamma(-x) = ------------------- ;
      196 + *              gamma(-x) = ------------------- ;
 196  197   *                          kpsin(z)*gamma(1+x)
 197  198   *
 198  199   *      otherwise, for z in (II),
 199  200   *                                      k+1
 200  201   *                                  (-1)
 201      - *              gamma(-x) = ----------------------- ;
      202 + *              gamma(-x) = ----------------------- ;
 202  203   *                          kpcos(0.5-z)*gamma(1+x)
 203  204   *
 204  205   *      otherwise, for z in (III),
 205  206   *                                      k+1
 206  207   *                                  (-1)
 207      - *              gamma(-x) = --------------------- .
      208 + *              gamma(-x) = --------------------- .
 208  209   *                          kpsin(1-z)*gamma(1+x)
 209  210   *
 210  211   *      Thus, the computation of gamma(-x) reduced to the computation of
 211  212   *      gamma(1+x) and kpsin(), kpcos().
 212  213   *
 213  214   * (B) For x between 1 and 2.  We break [1,2] into three parts:
 214  215   *      GT1 = [1.0000, 1.2845]
 215      - *      GT2 = [1.2844, 1.6374]
 216      - *      GT3 = [1.6373, 2.0000]
      216 + *      GT2 = [1.2844, 1.6374]
      217 + *      GT3 = [1.6373, 2.0000]
 217  218   *
 218  219   *    For x in GTi, i=1,2,3, let
 219      - *      z1  =  1.134861805732790769689793935774652917006
      220 + *      z1  =  1.134861805732790769689793935774652917006
 220  221   *      gz1 = gamma(z1)  =   0.9382046279096824494097535615803269576988
 221  222   *      tz1 = gamma'(z1) =  -0.3517214357852935791015625000000000000000
 222  223   *
 223  224   *      z2  =  1.461632144968362341262659542325721328468e+0000
 224  225   *      gz2 = gamma(z2)  = 0.8856031944108887002788159005825887332080
 225  226   *      tz2 = gamma'(z2) = 0.00
 226  227   *
 227  228   *      z3  =  1.819773101100500601787868704921606996312e+0000
 228  229   *      gz3 = gamma(z3)  = 0.9367814114636523216188468970808378497426
 229  230   *      tz3 = gamma'(z3) = 0.2805306315422058105468750000000000000000
↓ open down ↓ 83 lines elided ↑ open up ↑
 313  314   *       p2[5] =  -4.88158265593355093703112238534484636193260459574e-0002
 314  315   *
 315  316   *      i = 3
 316  317   *       p3[0] =   3.82409531118807759081121479786092134814808872880e-0001
 317  318   *       p3[1] =   2.65309888180188647956400403013495759365167853426e-0002
 318  319   *       p3[2] =   8.06815109775079171923561169415370309376296739835e-0002
 319  320   *       p3[3] =  -1.54821591666137613928840890835174351674007764799e-0002
 320  321   *       p3[4] =   1.76308239242717268530498313416899188157165183405e-0002
 321  322   *
 322  323   *    Coefficents: Double precision
 323      - *      i = 1:
      324 + *      i = 1:
 324  325   *       p1[0]   =   0.70908683619977797008004927192814648151397705078125000
 325  326   *       p1[1]   =   1.71987061393048558089579513384356441668351720061e-0001
 326  327   *       p1[2]   =  -3.19273345791990970293320316122813960527705450671e-0002
 327  328   *       p1[3]   =   8.36172645419110036267169600390549973563534476989e-0003
 328  329   *       p1[4]   =   1.13745336648572838333152213474277971244629758101e-0003
 329  330   *       q1[0]   =   1.0
 330  331   *       q1[1]   =   9.71980217826032937526460731778472389791321968082e-0001
 331  332   *       q1[2]   =  -7.43576743326756176594084137256042653497087666030e-0002
 332  333   *       q1[3]   =  -1.19345944932265559769719470515102012246995255372e-0001
 333  334   *       q1[4]   =   1.59913445751425002620935120470781382215050284762e-0002
 334  335   *       q1[5]   =   1.12601136853374984566572691306402321911547550783e-0003
 335      - *      i = 2:
      336 + *      i = 2:
 336  337   *       p2[0]   =   0.42848681585558601181418225678498856723308563232421875
 337  338   *       p2[1]   =   6.53596762668970816023718845105667418483122103629e-0002
 338  339   *       p2[2]   =  -6.97280829631212931321050770925128264272768936731e-0003
 339  340   *       p2[3]   =   6.46342359021981718947208605674813260166116632899e-0003
 340  341   *       q2[0]   =   1.0
 341  342   *       q2[1]   =   4.57572620560506047062553957454062012327519313936e-0001
 342  343   *       q2[2]   =  -2.52182594886075452859655003407796103083422572036e-0001
 343  344   *       q2[3]   =  -1.82970945407778594681348166040103197178711552827e-0002
 344  345   *       q2[4]   =   2.43574726993169566475227642128830141304953840502e-0002
 345  346   *       q2[5]   =  -5.20390406466942525358645957564897411258667085501e-0003
 346  347   *       q2[6]   =   4.79520251383279837635552431988023256031951133885e-0004
 347      - *      i = 3:
      348 + *      i = 3:
 348  349   *       p3[0]   =   0.382409479734567459008331979930517263710498809814453125
 349  350   *       p3[1]   =   1.42876048697668161599069814043449301572928034140e-0001
 350  351   *       p3[2]   =   3.42157571052250536817923866013561760785748899071e-0003
 351  352   *       p3[3]   =  -5.01542621710067521405087887856991700987709272937e-0004
 352  353   *       p3[4]   =   8.89285814866740910123834688163838287618332122670e-0004
 353  354   *       q3[0]   =   1.0
 354  355   *       q3[1]   =   3.04253086629444201002215640948957897906299633168e-0001
 355  356   *       q3[2]   =  -2.23162407379999477282555672834881213873185520006e-0001
 356  357   *       q3[3]   =  -1.05060867741952065921809811933670131427552903636e-0002
 357  358   *       q3[4]   =   1.70511763916186982473301861980856352005926669320e-0002
 358  359   *       q3[5]   =  -2.12950201683609187927899416700094630764182477464e-0003
 359  360   *
 360  361   *    Note that all pi0 are exact in double, which is obtained by a
 361  362   *    special Remez Algorithm.
 362  363   *
 363  364   *    Coefficents: Quad precision
 364      - *      i = 1:
      365 + *      i = 1:
 365  366   *       p1[0] =   0.709086836199777919037185741507610124611513720557
 366  367   *       p1[1] =   4.45754781206489035827915969367354835667391606951e-0001
 367  368   *       p1[2] =   3.21049298735832382311662273882632210062918153852e-0002
 368  369   *       p1[3] =  -5.71296796342106617651765245858289197369688864350e-0003
 369  370   *       p1[4] =   6.04666892891998977081619174969855831606965352773e-0003
 370  371   *       p1[5] =   8.99106186996888711939627812174765258822658645168e-0004
 371  372   *       p1[6] =  -6.96496846144407741431207008527018441810175568949e-0005
 372  373   *       p1[7] =   1.52597046118984020814225409300131445070213882429e-0005
 373  374   *       p1[8] =   5.68521076168495673844711465407432189190681541547e-0007
 374  375   *       p1[9] =   3.30749673519634895220582062520286565610418952979e-0008
 375  376   *       q1[0] =   1.0+0000
 376  377   *       q1[1] =   1.35806511721671070408570853537257079579490650668e+0000
 377  378   *       q1[2] =   2.97567810153429553405327140096063086994072952961e-0001
 378  379   *       q1[3] =  -1.52956835982588571502954372821681851681118097870e-0001
 379  380   *       q1[4] =  -2.88248519561420109768781615289082053597954521218e-0002
 380  381   *       q1[5] =   1.03475311719937405219789948456313936302378395955e-0002
 381  382   *       q1[6] =   4.12310203243891222368965360124391297374822742313e-0004
 382  383   *       q1[7] =  -3.12653708152290867248931925120380729518332507388e-0004
 383  384   *       q1[8] =   2.36672170850409745237358105667757760527014332458e-0005
 384  385   *
 385      - *      i = 2:
      386 + *      i = 2:
 386  387   *       p2[0] =   0.428486815855585429730209907810650616737756697477
 387  388   *       p2[1] =   2.63622124067885222919192651151581541943362617352e-0001
 388  389   *       p2[2] =   3.85520683670028865731877276741390421744971446855e-0002
 389  390   *       p2[3] =   3.05065978278128549958897133190295325258023525862e-0003
 390  391   *       p2[4] =   2.48232934951723128892080415054084339152450445081e-0003
 391  392   *       p2[5] =   3.67092777065632360693313762221411547741550105407e-0004
 392  393   *       p2[6] =   3.81228045616085789674530902563145250532194518946e-0006
 393  394   *       p2[7] =   4.61677225867087554059531455133839175822537617677e-0006
 394  395   *       p2[8] =   2.18209052385703200438239200991201916609364872993e-0007
 395  396   *       p2[9] =   1.00490538985245846460006244065624754421022542454e-0008
↓ open down ↓ 40 lines elided ↑ open up ↑
 436  437   *                1                       2
 437  438   *      Since  --------  ~  x + 0.577...*x  - ...,  we have, for small x,
 438  439   *              gamma(x)
 439  440   *           1                    1
 440  441   *      ----------- < gamma(x) < --- and
 441  442   *      x(1+0.578x)               x
 442  443   *              1                 1           1
 443  444   *        0 <  --- - gamma(x) <= ---  -  ----------- < 0.578
 444  445   *              x                 x      x(1+0.578x)
 445  446   *                                     1       1                        -P
 446      - *      The error is thus bounded by --- ulp(---) + 0.578. Since x <= 2   ,
      447 + *      The error is thus bounded by --- ulp(---) + 0.578. Since x <= 2   ,
 447  448   *                                     2       x
 448  449   *       1      P       1           P                                      1
 449  450   *      --- >= 2 , ulp(---) >= ulp(2  ) >= 2. Thus 0.578=0.289*2<=0.289ulp(-)
 450  451   *       x              x                                                  x
 451  452   *       Thus
 452  453   *                             1                                 1
 453  454   *              | gamma(x) - [---] rounded | <= (0.5+0.289)*ulp(---).
 454  455   *                             x                                 x
 455  456   *                         -P                              1
 456  457   *      Note that for x<= 2  , it is easy to see that ulp(---)=ulp(gamma(x))
 457  458   *                                                         x
 458  459   *                            n                             1
 459  460   *      except only when x = 2 , (n<= -53). In such cases, --- is exact
 460  461   *                                                          x
 461      - *      and therefore the error is bounded by
      462 + *      and therefore the error is bounded by
 462  463   *                         1
 463  464   *              0.298*ulp(---) = 0.298*2*ulp(gamma(x)) = 0.578ulp(gamma(x)).
 464  465   *                         x
 465  466   *      Thus we conclude that the error in gamma is less than 0.739 ulp.
 466  467   *
 467  468   *    (2)Otherwise, for x in GTi-1 (see B), let y = x-(zi-1). From (B) we obtain
 468  469   *                                                          gamma(1+x)
 469  470   *      gamma(1+x) = gy.h + gy.l,  then compute gamma(x) by -----------.
 470  471   *                                                               x
 471  472   *                                                          gy.h
↓ open down ↓ 63 lines elided ↑ open up ↑
 535  536   *      log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x))
 536  537   *                = L1 + L2 + L3,
 537  538   *    where
 538  539   *              L1(x) = (x-.5)*(log(x)-1),
 539  540   *              L2    = .5(log(2pi)-1) = 0.41893853....,
 540  541   *              L3(x) = (1/x)P(1/(x*x)),
 541  542   *
 542  543   *    The range of L1,L2, and L3 are as follows:
 543  544   *
 544  545   *      ------------------------------------------------------------------
 545      - *      Range(L1) =  (single) [8.09..,88.30..]   =[2** 3.01..,2**  6.46..]
      546 + *      Range(L1) =  (single) [8.09..,88.30..]   =[2** 3.01..,2**  6.46..]
 546  547   *                   (double) [8.09..,709.3..]   =[2** 3.01..,2**  9.47..]
 547  548   *                   (quad)   [8.09..,11356.10..]=[2** 3.01..,2** 13.47..]
 548      - *      Range(L2) = 0.41893853.....
      549 + *      Range(L2) = 0.41893853.....
 549  550   *      Range(L3) = [0.0104...., 0.00048....]    =[2**-6.58..,2**-11.02..]
 550  551   *      ------------------------------------------------------------------
 551  552   *
 552  553   *    Gamma(x) is then computed by exp(L1+L2+L3).
 553  554   *
 554  555   *    (2) Error analysis of (F):
 555  556   *    --------------------------
 556  557   *    The error in Gamma(x) depends on the error inherited in the computation
 557  558   *    of L= L1+L2+L3. Let L' be the computed value of L. The absolute error
 558  559   *    in L' is t = L-L'. Since exp(L') = exp(L-t) = exp(L)*exp(t) ~
↓ open down ↓ 46 lines elided ↑ open up ↑
 605  606   *
 606  607   *    where T1(n) = n*log(2)-1 and T2(j) = log(z(j)). Both T1 and T2 can be
 607  608   *    pre-calculated and be looked-up in a table. Note that 8 <= x < 1756
 608  609   *    implies 3<=n<=10 implies 1.079.. < T1(n) < 6.931.
 609  610   *
 610  611   *
 611  612   *                     y-z(i)          y       1+s
 612  613   *    For T3, let s = --------; then ----- =  ----- and
 613  614   *                     y+z(i)         z(i)     1-s
 614  615   *                1+s           2   3    2   5
 615      - *      T3 = log(-----) = 2s + --- s  + --- s  + ....
      616 + *      T3 = log(-----) = 2s + --- s  + --- s  + ....
 616  617   *                1-s           3        5
 617  618   *
 618  619   *    Suppose the first term 2s is compute in extra precision. The
 619  620   *    dominating error in T3 would then be the rounding error of the
 620  621   *    second term 2/3*s**3. To force the rounding bounded by
 621  622   *    the required accuracy, we have
 622  623   *        single:  |2/3*s**3| < 2**-11   == > |s|<0.09014...
 623  624   *        double:  |2/3*s**3| < 2**-14   == > |s|<0.04507...
 624  625   *        quad  :  |2/3*s**3| < 2**-18   == > |s|<0.01788... = 2**(-5.80..)
 625  626   *
↓ open down ↓ 41 lines elided ↑ open up ↑
 667  668   *    (3) Remez approximation for (T3(s)-2s)/s = T3[0]*s + T3[1]*s + ...:
 668  669   *       Single precision: 1 term (compute in double precision arithmetic)
 669  670   *          T3(s) = 2s + S1*s^3, S1 = 0.6666717231848518054693623697539230
 670  671   *          Remez error: |T3(s)/s - (2s+S1*s^3)| < 2**(-35.87)
 671  672   *       Double precision: 3 terms, Remez error is bounded by 2**(-72.40),
 672  673   *          see "tgamma_log"
 673  674   *       Quad precision: 7 terms, Remez error is bounded by 2**(-136.54),
 674  675   *          see "tgammal_log"
 675  676   *
 676  677   *   The computation of 0.5*(ln(2pi)-1):
 677      - *      0.5*(ln(2pi)-1) =  0.4189385332046727417803297364056176398614...
      678 + *      0.5*(ln(2pi)-1) =  0.4189385332046727417803297364056176398614...
 678  679   *      split 0.5*(ln(2pi)-1) to hln2pi_h + hln2pi_l, where hln2pi_h is the
 679  680   *      leading 21 bits of the constant.
 680  681   *          hln2pi_h= 0.4189383983612060546875
 681  682   *          hln2pi_l= 1.348434666870928297364056176398612173648e-07
 682  683   *
 683  684   *   The computation of 1/x*P(1/x^2) = log(G(x))-(x-.5)(ln(x)-1)-(.5ln(2pi)-1):
 684  685   *      Let s = 1/x <= 1/8 < 0.125. We have
 685  686   *      quad precision
 686  687   *          |GP(s) - s*P(s^2)| <= 2**(-120.6), where
 687  688   *                             3      5            39
↓ open down ↓ 41 lines elided ↑ open up ↑
 729  730   *        GP1 =  -2.77765545601667179767706600890361535225507762168e-0003
 730  731   *        GP2 =   7.77830853479775281781085278324621033523037489883e-0004
 731  732   *
 732  733   *
 733  734   *      Implementation note:
 734  735   *      z = (1/x), z2 = z*z, z4 = z2*z2;
 735  736   *      p = z*(GP0+z2*(GP1+....+z2*GP7))
 736  737   *        = z*(GP0+(z4*(GP2+z4*(GP4+z4*GP6))+z2*(GP1+z4*(GP3+z4*(GP5+z4*GP7)))))
 737  738   *
 738  739   *   Adding everything up:
 739      - *      t = rr.h*ww.h+hln2pi_h                  ... exact
      740 + *      t = rr.h*ww.h+hln2pi_h                  ... exact
 740  741   *      w = (hln2pi_l + ((x-0.5)*ww.l+rr.l*ww.h)) + p
 741  742   *
 742  743   *   Computing exp(t+w):
 743  744   *      s = t+w; write s = (n+j/32)*ln2+r, |r|<=(1/64)*ln2, then
 744  745   *      exp(s) = 2**n * (2**(j/32) + 2**(j/32)*expm1(r)), where
 745  746   *      expm1(r) = r + Et1*r^2 + Et2*r^3 + ... + Et5*r^6, and
 746  747   *      2**(j/32) is obtained by table look-up S[j]+S_trail[j].
 747  748   *      Remez error bound:
 748  749   *      |exp(r) - (1+r+Et1*r^2+...+Et5*r^6)| <= 2^(-63).
 749  750   */
      751 +/* END CSTYLED */
 750  752  
 751  753  #include "libm.h"
 752  754  
 753      -#define __HI(x) ((int *) &x)[HIWORD]
 754      -#define __LO(x) ((unsigned *) &x)[LOWORD]
      755 +#define __HI(x)         ((int *)&x)[HIWORD]
      756 +#define __LO(x)         ((unsigned *)&x)[LOWORD]
 755  757  
 756  758  struct Double {
 757  759          double h;
 758  760          double l;
 759  761  };
 760  762  
 761  763  /* Hex value of GP0 shoule be 3FB55555 55555555 */
 762  764  static const double c[] = {
 763  765          +1.0,
 764  766          +2.0,
 765  767          +0.5,
 766  768          +1.0e-300,
 767      -        +6.66666666666666740682e-01,                            /* A1=T3[0] */
 768      -        +3.99999999955626478023093908674902212920e-01,          /* A2=T3[1] */
 769      -        +2.85720221533145659809237398709372330980e-01,          /* A3=T3[2] */
 770      -        +0.0833333333333333287074040640618477,                  /* GP[0] */
      769 +        +6.66666666666666740682e-01,    /* A1=T3[0] */
      770 +        +3.99999999955626478023093908674902212920e-01,  /* A2=T3[1] */
      771 +        +2.85720221533145659809237398709372330980e-01,  /* A3=T3[2] */
      772 +        +0.0833333333333333287074040640618477,          /* GP[0] */
 771  773          -2.77777777776649355200565611114627670089130772843e-03,
 772  774          +7.93650787486083724805476194170211775784158551509e-04,
 773  775          -5.95236628558314928757811419580281294593903582971e-04,
 774  776          +8.41566473999853451983137162780427812781178932540e-04,
 775  777          -1.90424776670441373564512942038926168175921303212e-03,
 776  778          +5.84933161530949666312333949534482303007354299178e-03,
 777  779          -1.59453228931082030262124832506144392496561694550e-02,
 778  780          +4.18937683105468750000e-01,                            /* hln2pi_h */
 779  781          +8.50099203991780279640e-07,                            /* hln2pi_l */
 780  782          +4.18938533204672741744150788368695779923320328369e-01, /* hln2pi */
 781  783          +2.16608493865351192653e-02,                            /* ln2_32hi */
 782  784          +5.96317165397058656257e-12,                            /* ln2_32lo */
 783  785          +4.61662413084468283841e+01,                            /* invln2_32 */
 784  786          +5.0000000000000000000e-1,                              /* Et1 */
 785  787          +1.66666666665223585560605991943703896196054020060e-01, /* Et2 */
 786  788          +4.16666666665895103520154073534275286743788421687e-02, /* Et3 */
 787  789          +8.33336844093536520775865096538773197505523826029e-03, /* Et4 */
 788  790          +1.38889201930843436040204096950052984793587640227e-03, /* Et5 */
 789  791  };
 790  792  
 791      -#define one       c[0]
 792      -#define two       c[1]
 793      -#define half      c[2]
 794      -#define tiny      c[3]
 795      -#define A1        c[4]
 796      -#define A2        c[5]
 797      -#define A3        c[6]
 798      -#define GP0       c[7]
 799      -#define GP1       c[8]
 800      -#define GP2       c[9]
 801      -#define GP3       c[10]
 802      -#define GP4       c[11]
 803      -#define GP5       c[12]
 804      -#define GP6       c[13]
 805      -#define GP7       c[14]
 806      -#define hln2pi_h  c[15]
 807      -#define hln2pi_l  c[16]
 808      -#define hln2pi    c[17]
 809      -#define ln2_32hi  c[18]
 810      -#define ln2_32lo  c[19]
 811      -#define invln2_32 c[20]
 812      -#define Et1       c[21]
 813      -#define Et2       c[22]
 814      -#define Et3       c[23]
 815      -#define Et4       c[24]
 816      -#define Et5       c[25]
      793 +#define one             c[0]
      794 +#define two             c[1]
      795 +#define half            c[2]
      796 +#define tiny            c[3]
      797 +#define A1              c[4]
      798 +#define A2              c[5]
      799 +#define A3              c[6]
      800 +#define GP0             c[7]
      801 +#define GP1             c[8]
      802 +#define GP2             c[9]
      803 +#define GP3             c[10]
      804 +#define GP4             c[11]
      805 +#define GP5             c[12]
      806 +#define GP6             c[13]
      807 +#define GP7             c[14]
      808 +#define hln2pi_h        c[15]
      809 +#define hln2pi_l        c[16]
      810 +#define hln2pi          c[17]
      811 +#define ln2_32hi        c[18]
      812 +#define ln2_32lo        c[19]
      813 +#define invln2_32       c[20]
      814 +#define Et1             c[21]
      815 +#define Et2             c[22]
      816 +#define Et3             c[23]
      817 +#define Et4             c[24]
      818 +#define Et5             c[25]
 817  819  
 818  820  /*
 819  821   * double precision coefficients for computing log(x)-1 in tgamma.
 820  822   *  See "algorithm" for details
 821  823   *
 822  824   *  log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y,  1<=y<2,
 823  825   *  j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
 824  826   *       T1(n) = T1[2n,2n+1] = n*log(2)-1,
 825  827   *       T2(j) = T2[2j,2j+1] = log(z[j]),
 826  828   *       T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + T3[2]s^7
↓ open down ↓ 156 lines elided ↑ open up ↑
 983  985          +6.73422634601593017578e-01,    /* 0x3FE58CAD 0xA0000000 */
 984  986          +4.06105737027198329700e-08,    /* 0x3E65CD79 0x893092F2 */
 985  987          +6.81359171867370605469e-01,    /* 0x3FE5CDB1 0xC0000000 */
 986  988          +5.29405324634793230630e-08,    /* 0x3E6C6C17 0x648CF6E4 */
 987  989          +6.89233243465423583984e-01,    /* 0x3FE60E32 0xE0000000 */
 988  990          +3.77733853963405370102e-08,    /* 0x3E644788 0xD8CA7C89 */
 989  991  };
 990  992  
 991  993  /* S[j],S_trail[j] = 2**(j/32.) for the final computation of exp(t+w) */
 992  994  static const double S[] = {
 993      -        +1.00000000000000000000e+00,    /* 3FF0000000000000 */
 994      -        +1.02189714865411662714e+00,    /* 3FF059B0D3158574 */
 995      -        +1.04427378242741375480e+00,    /* 3FF0B5586CF9890F */
 996      -        +1.06714040067682369717e+00,    /* 3FF11301D0125B51 */
 997      -        +1.09050773266525768967e+00,    /* 3FF172B83C7D517B */
 998      -        +1.11438674259589243221e+00,    /* 3FF1D4873168B9AA */
 999      -        +1.13878863475669156458e+00,    /* 3FF2387A6E756238 */
1000      -        +1.16372485877757747552e+00,    /* 3FF29E9DF51FDEE1 */
1001      -        +1.18920711500272102690e+00,    /* 3FF306FE0A31B715 */
1002      -        +1.21524735998046895524e+00,    /* 3FF371A7373AA9CB */
1003      -        +1.24185781207348400201e+00,    /* 3FF3DEA64C123422 */
1004      -        +1.26905095719173321989e+00,    /* 3FF44E086061892D */
1005      -        +1.29683955465100964055e+00,    /* 3FF4BFDAD5362A27 */
1006      -        +1.32523664315974132322e+00,    /* 3FF5342B569D4F82 */
1007      -        +1.35425554693689265129e+00,    /* 3FF5AB07DD485429 */
1008      -        +1.38390988196383202258e+00,    /* 3FF6247EB03A5585 */
1009      -        +1.41421356237309514547e+00,    /* 3FF6A09E667F3BCD */
1010      -        +1.44518080697704665027e+00,    /* 3FF71F75E8EC5F74 */
1011      -        +1.47682614593949934623e+00,    /* 3FF7A11473EB0187 */
1012      -        +1.50916442759342284141e+00,    /* 3FF82589994CCE13 */
1013      -        +1.54221082540794074411e+00,    /* 3FF8ACE5422AA0DB */
1014      -        +1.57598084510788649659e+00,    /* 3FF93737B0CDC5E5 */
1015      -        +1.61049033194925428347e+00,    /* 3FF9C49182A3F090 */
1016      -        +1.64575547815396494578e+00,    /* 3FFA5503B23E255D */
1017      -        +1.68179283050742900407e+00,    /* 3FFAE89F995AD3AD */
1018      -        +1.71861929812247793414e+00,    /* 3FFB7F76F2FB5E47 */
1019      -        +1.75625216037329945351e+00,    /* 3FFC199BDD85529C */
1020      -        +1.79470907500310716820e+00,    /* 3FFCB720DCEF9069 */
1021      -        +1.83400808640934243066e+00,    /* 3FFD5818DCFBA487 */
1022      -        +1.87416763411029996256e+00,    /* 3FFDFC97337B9B5F */
1023      -        +1.91520656139714740007e+00,    /* 3FFEA4AFA2A490DA */
1024      -        +1.95714412417540017941e+00,    /* 3FFF50765B6E4540 */
      995 +        +1.00000000000000000000e+00, /* 3FF0000000000000 */
      996 +        +1.02189714865411662714e+00, /* 3FF059B0D3158574 */
      997 +        +1.04427378242741375480e+00, /* 3FF0B5586CF9890F */
      998 +        +1.06714040067682369717e+00, /* 3FF11301D0125B51 */
      999 +        +1.09050773266525768967e+00, /* 3FF172B83C7D517B */
     1000 +        +1.11438674259589243221e+00, /* 3FF1D4873168B9AA */
     1001 +        +1.13878863475669156458e+00, /* 3FF2387A6E756238 */
     1002 +        +1.16372485877757747552e+00, /* 3FF29E9DF51FDEE1 */
     1003 +        +1.18920711500272102690e+00, /* 3FF306FE0A31B715 */
     1004 +        +1.21524735998046895524e+00, /* 3FF371A7373AA9CB */
     1005 +        +1.24185781207348400201e+00, /* 3FF3DEA64C123422 */
     1006 +        +1.26905095719173321989e+00, /* 3FF44E086061892D */
     1007 +        +1.29683955465100964055e+00, /* 3FF4BFDAD5362A27 */
     1008 +        +1.32523664315974132322e+00, /* 3FF5342B569D4F82 */
     1009 +        +1.35425554693689265129e+00, /* 3FF5AB07DD485429 */
     1010 +        +1.38390988196383202258e+00, /* 3FF6247EB03A5585 */
     1011 +        +1.41421356237309514547e+00, /* 3FF6A09E667F3BCD */
     1012 +        +1.44518080697704665027e+00, /* 3FF71F75E8EC5F74 */
     1013 +        +1.47682614593949934623e+00, /* 3FF7A11473EB0187 */
     1014 +        +1.50916442759342284141e+00, /* 3FF82589994CCE13 */
     1015 +        +1.54221082540794074411e+00, /* 3FF8ACE5422AA0DB */
     1016 +        +1.57598084510788649659e+00, /* 3FF93737B0CDC5E5 */
     1017 +        +1.61049033194925428347e+00, /* 3FF9C49182A3F090 */
     1018 +        +1.64575547815396494578e+00, /* 3FFA5503B23E255D */
     1019 +        +1.68179283050742900407e+00, /* 3FFAE89F995AD3AD */
     1020 +        +1.71861929812247793414e+00, /* 3FFB7F76F2FB5E47 */
     1021 +        +1.75625216037329945351e+00, /* 3FFC199BDD85529C */
     1022 +        +1.79470907500310716820e+00, /* 3FFCB720DCEF9069 */
     1023 +        +1.83400808640934243066e+00, /* 3FFD5818DCFBA487 */
     1024 +        +1.87416763411029996256e+00, /* 3FFDFC97337B9B5F */
     1025 +        +1.91520656139714740007e+00, /* 3FFEA4AFA2A490DA */
     1026 +        +1.95714412417540017941e+00, /* 3FFF50765B6E4540 */
1025 1027  };
1026 1028  
1027 1029  static const double S_trail[] = {
1028 1030          +0.00000000000000000000e+00,
1029      -        +5.10922502897344389359e-17,    /* 3C8D73E2A475B465 */
1030      -        +8.55188970553796365958e-17,    /* 3C98A62E4ADC610A */
1031      -        -7.89985396684158212226e-17,    /* BC96C51039449B3A */
1032      -        -3.04678207981247114697e-17,    /* BC819041B9D78A76 */
1033      -        +1.04102784568455709549e-16,    /* 3C9E016E00A2643C */
1034      -        +8.91281267602540777782e-17,    /* 3C99B07EB6C70573 */
1035      -        +3.82920483692409349872e-17,    /* 3C8612E8AFAD1255 */
1036      -        +3.98201523146564611098e-17,    /* 3C86F46AD23182E4 */
1037      -        -7.71263069268148813091e-17,    /* BC963AEABF42EAE2 */
1038      -        +4.65802759183693679123e-17,    /* 3C8ADA0911F09EBC */
1039      -        +2.66793213134218609523e-18,    /* 3C489B7A04EF80D0 */
1040      -        +2.53825027948883149593e-17,    /* 3C7D4397AFEC42E2 */
1041      -        -2.85873121003886075697e-17,    /* BC807ABE1DB13CAC */
1042      -        +7.70094837980298946162e-17,    /* 3C96324C054647AD */
1043      -        -6.77051165879478628716e-17,    /* BC9383C17E40B497 */
1044      -        -9.66729331345291345105e-17,    /* BC9BDD3413B26456 */
1045      -        -3.02375813499398731940e-17,    /* BC816E4786887A99 */
1046      -        -3.48399455689279579579e-17,    /* BC841577EE04992F */
1047      -        -1.01645532775429503911e-16,    /* BC9D4C1DD41532D8 */
1048      -        +7.94983480969762085616e-17,    /* 3C96E9F156864B27 */
1049      -        -1.01369164712783039808e-17,    /* BC675FC781B57EBC */
1050      -        +2.47071925697978878522e-17,    /* 3C7C7C46B071F2BE */
1051      -        -1.01256799136747726038e-16,    /* BC9D2F6EDB8D41E1 */
1052      -        +8.19901002058149652013e-17,    /* 3C97A1CD345DCC81 */
1053      -        -1.85138041826311098821e-17,    /* BC75584F7E54AC3B */
1054      -        +2.96014069544887330703e-17,    /* 3C811065895048DD */
1055      -        +1.82274584279120867698e-17,    /* 3C7503CBD1E949DB */
1056      -        +3.28310722424562658722e-17,    /* 3C82ED02D75B3706 */
1057      -        -6.12276341300414256164e-17,    /* BC91A5CD4F184B5C */
1058      -        -1.06199460561959626376e-16,    /* BC9E9C23179C2893 */
1059      -        +8.96076779103666776760e-17,    /* 3C99D3E12DD8A18B */
     1031 +        +5.10922502897344389359e-17, /* 3C8D73E2A475B465 */
     1032 +        +8.55188970553796365958e-17, /* 3C98A62E4ADC610A */
     1033 +        -7.89985396684158212226e-17, /* BC96C51039449B3A */
     1034 +        -3.04678207981247114697e-17, /* BC819041B9D78A76 */
     1035 +        +1.04102784568455709549e-16, /* 3C9E016E00A2643C */
     1036 +        +8.91281267602540777782e-17, /* 3C99B07EB6C70573 */
     1037 +        +3.82920483692409349872e-17, /* 3C8612E8AFAD1255 */
     1038 +        +3.98201523146564611098e-17, /* 3C86F46AD23182E4 */
     1039 +        -7.71263069268148813091e-17, /* BC963AEABF42EAE2 */
     1040 +        +4.65802759183693679123e-17, /* 3C8ADA0911F09EBC */
     1041 +        +2.66793213134218609523e-18, /* 3C489B7A04EF80D0 */
     1042 +        +2.53825027948883149593e-17, /* 3C7D4397AFEC42E2 */
     1043 +        -2.85873121003886075697e-17, /* BC807ABE1DB13CAC */
     1044 +        +7.70094837980298946162e-17, /* 3C96324C054647AD */
     1045 +        -6.77051165879478628716e-17, /* BC9383C17E40B497 */
     1046 +        -9.66729331345291345105e-17, /* BC9BDD3413B26456 */
     1047 +        -3.02375813499398731940e-17, /* BC816E4786887A99 */
     1048 +        -3.48399455689279579579e-17, /* BC841577EE04992F */
     1049 +        -1.01645532775429503911e-16, /* BC9D4C1DD41532D8 */
     1050 +        +7.94983480969762085616e-17, /* 3C96E9F156864B27 */
     1051 +        -1.01369164712783039808e-17, /* BC675FC781B57EBC */
     1052 +        +2.47071925697978878522e-17, /* 3C7C7C46B071F2BE */
     1053 +        -1.01256799136747726038e-16, /* BC9D2F6EDB8D41E1 */
     1054 +        +8.19901002058149652013e-17, /* 3C97A1CD345DCC81 */
     1055 +        -1.85138041826311098821e-17, /* BC75584F7E54AC3B */
     1056 +        +2.96014069544887330703e-17, /* 3C811065895048DD */
     1057 +        +1.82274584279120867698e-17, /* 3C7503CBD1E949DB */
     1058 +        +3.28310722424562658722e-17, /* 3C82ED02D75B3706 */
     1059 +        -6.12276341300414256164e-17, /* BC91A5CD4F184B5C */
     1060 +        -1.06199460561959626376e-16, /* BC9E9C23179C2893 */
     1061 +        +8.96076779103666776760e-17, /* 3C99D3E12DD8A18B */
1060 1062  };
1061 1063  
1062 1064  /* Primary interval GTi() */
1063 1065  static const double cr[] = {
1064 1066  /* p1, q1 */
1065 1067          +0.70908683619977797008004927192814648151397705078125000,
1066 1068          +1.71987061393048558089579513384356441668351720061e-0001,
1067 1069          -3.19273345791990970293320316122813960527705450671e-0002,
1068 1070          +8.36172645419110036267169600390549973563534476989e-0003,
1069 1071          +1.13745336648572838333152213474277971244629758101e-0003,
↓ open down ↓ 22 lines elided ↑ open up ↑
1092 1094          -5.01542621710067521405087887856991700987709272937e-0004,
1093 1095          +8.89285814866740910123834688163838287618332122670e-0004,
1094 1096          +1.0,
1095 1097          +3.04253086629444201002215640948957897906299633168e-0001,
1096 1098          -2.23162407379999477282555672834881213873185520006e-0001,
1097 1099          -1.05060867741952065921809811933670131427552903636e-0002,
1098 1100          +1.70511763916186982473301861980856352005926669320e-0002,
1099 1101          -2.12950201683609187927899416700094630764182477464e-0003,
1100 1102  };
1101 1103  
1102      -#define P10   cr[0]
1103      -#define P11   cr[1]
1104      -#define P12   cr[2]
1105      -#define P13   cr[3]
1106      -#define P14   cr[4]
1107      -#define Q10   cr[5]
1108      -#define Q11   cr[6]
1109      -#define Q12   cr[7]
1110      -#define Q13   cr[8]
1111      -#define Q14   cr[9]
1112      -#define Q15   cr[10]
1113      -#define P20   cr[11]
1114      -#define P21   cr[12]
1115      -#define P22   cr[13]
1116      -#define P23   cr[14]
1117      -#define Q20   cr[15]
1118      -#define Q21   cr[16]
1119      -#define Q22   cr[17]
1120      -#define Q23   cr[18]
1121      -#define Q24   cr[19]
1122      -#define Q25   cr[20]
1123      -#define Q26   cr[21]
1124      -#define P30   cr[22]
1125      -#define P31   cr[23]
1126      -#define P32   cr[24]
1127      -#define P33   cr[25]
1128      -#define P34   cr[26]
1129      -#define Q30   cr[27]
1130      -#define Q31   cr[28]
1131      -#define Q32   cr[29]
1132      -#define Q33   cr[30]
1133      -#define Q34   cr[31]
1134      -#define Q35   cr[32]
     1104 +#define P10             cr[0]
     1105 +#define P11             cr[1]
     1106 +#define P12             cr[2]
     1107 +#define P13             cr[3]
     1108 +#define P14             cr[4]
     1109 +#define Q10             cr[5]
     1110 +#define Q11             cr[6]
     1111 +#define Q12             cr[7]
     1112 +#define Q13             cr[8]
     1113 +#define Q14             cr[9]
     1114 +#define Q15             cr[10]
     1115 +#define P20             cr[11]
     1116 +#define P21             cr[12]
     1117 +#define P22             cr[13]
     1118 +#define P23             cr[14]
     1119 +#define Q20             cr[15]
     1120 +#define Q21             cr[16]
     1121 +#define Q22             cr[17]
     1122 +#define Q23             cr[18]
     1123 +#define Q24             cr[19]
     1124 +#define Q25             cr[20]
     1125 +#define Q26             cr[21]
     1126 +#define P30             cr[22]
     1127 +#define P31             cr[23]
     1128 +#define P32             cr[24]
     1129 +#define P33             cr[25]
     1130 +#define P34             cr[26]
     1131 +#define Q30             cr[27]
     1132 +#define Q31             cr[28]
     1133 +#define Q32             cr[29]
     1134 +#define Q33             cr[30]
     1135 +#define Q34             cr[31]
     1136 +#define Q35             cr[32]
1135 1137  
1136      -static const double
1137      -        GZ1_h = +0.938204627909682398190,
     1138 +static const double GZ1_h = +0.938204627909682398190,
1138 1139          GZ1_l = +5.121952600248205157935e-17,
1139 1140          GZ2_h = +0.885603194410888749921,
1140 1141          GZ2_l = -4.964236872556339810692e-17,
1141 1142          GZ3_h = +0.936781411463652347038,
1142 1143          GZ3_l = -2.541923110834479415023e-17,
1143 1144          TZ1 = -0.3517214357852935791015625,
1144 1145          TZ3 = +0.280530631542205810546875;
1145      -/* INDENT ON */
1146 1146  
1147      -/* compute gamma(y=yh+yl) for y in GT1 = [1.0000, 1.2845] */
1148      -/* assume yh got 20 significant bits */
     1147 +/*
     1148 + * compute gamma(y=yh+yl) for y in GT1 = [1.0000, 1.2845]
     1149 + * assume yh got 20 significant bits
     1150 + */
1149 1151  static struct Double
1150      -GT1(double yh, double yl) {
     1152 +GT1(double yh, double yl)
     1153 +{
1151 1154          double t3, t4, y, z;
1152 1155          struct Double r;
1153 1156  
1154 1157          y = yh + yl;
1155 1158          z = y * y;
1156      -        t3 = (z * (P10 + y * ((P11 + y * P12) + z * (P13 + y * P14)))) /
1157      -                (Q10 + y * ((Q11 + y * Q12) + z * ((Q13 + Q14 * y) + z * Q15)));
     1159 +        t3 = (z * (P10 + y * ((P11 + y * P12) + z * (P13 + y * P14)))) / (Q10 +
     1160 +            y * ((Q11 + y * Q12) + z * ((Q13 + Q14 * y) + z * Q15)));
1158 1161          t3 += (TZ1 * yl + GZ1_l);
1159 1162          t4 = TZ1 * yh;
1160      -        r.h = (double) ((float) (t4 + GZ1_h + t3));
     1163 +        r.h = (double)((float)(t4 + GZ1_h + t3));
1161 1164          t3 += (t4 - (r.h - GZ1_h));
1162 1165          r.l = t3;
1163 1166          return (r);
1164 1167  }
1165 1168  
1166      -/* compute gamma(y=yh+yl) for y in GT2 = [1.2844, 1.6374] */
1167      -/* assume yh got 20 significant bits */
     1169 +/*
     1170 + * compute gamma(y=yh+yl) for y in GT2 = [1.2844, 1.6374]
     1171 + * assume yh got 20 significant bits
     1172 + */
1168 1173  static struct Double
1169      -GT2(double yh, double yl) {
     1174 +GT2(double yh, double yl)
     1175 +{
1170 1176          double t3, y, z;
1171 1177          struct Double r;
1172 1178  
1173 1179          y = yh + yl;
1174 1180          z = y * y;
1175      -        t3 = (z * (P20 + y * P21 + z * (P22 + y * P23))) /
1176      -                (Q20 + (y * ((Q21 + Q22 * y) + z * Q23) +
1177      -                (z * z) * ((Q24 + Q25 * y) + z * Q26))) + GZ2_l;
1178      -        r.h = (double) ((float) (GZ2_h + t3));
     1181 +        t3 = (z * (P20 + y * P21 + z * (P22 + y * P23))) / (Q20 + (y * ((Q21 +
     1182 +            Q22 * y) + z * Q23) + (z * z) * ((Q24 + Q25 * y) + z * Q26))) +
     1183 +            GZ2_l;
     1184 +        r.h = (double)((float)(GZ2_h + t3));
1179 1185          r.l = t3 - (r.h - GZ2_h);
1180 1186          return (r);
1181 1187  }
1182 1188  
1183      -/* compute gamma(y=yh+yl) for y in GT3 = [1.6373, 2.0000] */
1184      -/* assume yh got 20 significant bits */
     1189 +/*
     1190 + * compute gamma(y=yh+yl) for y in GT3 = [1.6373, 2.0000]
     1191 + * assume yh got 20 significant bits
     1192 + */
1185 1193  static struct Double
1186      -GT3(double yh, double yl) {
     1194 +GT3(double yh, double yl)
     1195 +{
1187 1196          double t3, t4, y, z;
1188 1197          struct Double r;
1189 1198  
1190 1199          y = yh + yl;
1191 1200          z = y * y;
1192      -        t3 = (z * (P30 + y * ((P31 + y * P32) + z * (P33 + y * P34)))) /
1193      -                (Q30 + y * ((Q31 + y * Q32) + z * ((Q33 + Q34 * y) + z * Q35)));
     1201 +        t3 = (z * (P30 + y * ((P31 + y * P32) + z * (P33 + y * P34)))) / (Q30 +
     1202 +            y * ((Q31 + y * Q32) + z * ((Q33 + Q34 * y) + z * Q35)));
1194 1203          t3 += (TZ3 * yl + GZ3_l);
1195 1204          t4 = TZ3 * yh;
1196      -        r.h = (double) ((float) (t4 + GZ3_h + t3));
     1205 +        r.h = (double)((float)(t4 + GZ3_h + t3));
1197 1206          t3 += (t4 - (r.h - GZ3_h));
1198 1207          r.l = t3;
1199 1208          return (r);
1200 1209  }
1201 1210  
1202      -/* INDENT OFF */
     1211 +
1203 1212  /*
1204 1213   * return tgamma(x) scaled by 2**-m for 8<x<=171.62... using Stirling's formula
1205 1214   *     log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x))
1206 1215   *                = L1 + L2 + L3,
1207 1216   */
1208      -/* INDENT ON */
1209 1217  static struct Double
1210      -large_gam(double x, int *m) {
1211      -        double z, t1, t2, t3, z2, t5, w, y, u, r, z4, v, t24 = 16777216.0,
1212      -                p24 = 1.0 / 16777216.0;
     1218 +large_gam(double x, int *m)
     1219 +{
     1220 +        double z, t1, t2, t3, z2, t5, w, y, u, r, z4, v, t24 = 16777216.0, p24 =
     1221 +            1.0 / 16777216.0;
1213 1222          int n2, j2, k, ix, j;
1214 1223          unsigned lx;
1215 1224          struct Double zz;
1216 1225          double u2, ss_h, ss_l, r_h, w_h, w_l, t4;
1217 1226  
1218      -/* INDENT OFF */
     1227 +
1219 1228  /*
1220 1229   * compute ss = ss.h+ss.l = log(x)-1 (see tgamma_log.h for details)
1221 1230   *
1222 1231   *  log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y,  1<=y<2,
1223 1232   *  j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
1224 1233   *       T1(n) = T1[2n,2n+1] = n*log(2)-1,
1225 1234   *       T2(j) = T2[2j,2j+1] = log(z[j]),
1226 1235   *       T3(s) = 2s + A1[0]s^3 + A2[1]s^5 + A3[2]s^7
1227 1236   *  Note
1228 1237   *  (1) the leading entries are truncated to 24 binary point.
↓ open down ↓ 3 lines elided ↑ open up ↑
1232 1241   *               T1(n):     |_________|___________________|
1233 1242   *                             _______ ______________________
1234 1243   *               T2(j):       |_______|______________________|
1235 1244   *                                ____ _______________________
1236 1245   *               2s:             |____|_______________________|
1237 1246   *                                    __________________________
1238 1247   *          +    T3(s)-2s:           |__________________________|
1239 1248   *                       -------------------------------------------
1240 1249   *                          [leading] + [Trailing]
1241 1250   */
1242      -/* INDENT ON */
1243 1251          ix = __HI(x);
1244 1252          lx = __LO(x);
1245      -        n2 = (ix >> 20) - 0x3ff;        /* exponent of x, range:3-7 */
1246      -        n2 += n2;                       /* 2n */
     1253 +        n2 = (ix >> 20) - 0x3ff;                /* exponent of x, range:3-7 */
     1254 +        n2 += n2;                               /* 2n */
1247 1255          ix = (ix & 0x000fffff) | 0x3ff00000;    /* y = scale x to [1,2] */
1248 1256          __HI(y) = ix;
1249 1257          __LO(y) = lx;
1250 1258          __HI(z) = (ix & 0xffffc000) | 0x2000;   /* z[j]=1+j/64+1/128 */
1251 1259          __LO(z) = 0;
1252      -        j2 = (ix >> 13) & 0x7e; /* 2j */
     1260 +        j2 = (ix >> 13) & 0x7e;                 /* 2j */
1253 1261          t1 = y + z;
1254 1262          t2 = y - z;
1255 1263          r = one / t1;
1256      -        t1 = (double) ((float) t1);
1257      -        u = r * t2;             /* u = (y-z)/(y+z) */
     1264 +        t1 = (double)((float)t1);
     1265 +        u = r * t2;                     /* u = (y-z)/(y+z) */
1258 1266          t4 = T2[j2 + 1] + T1[n2 + 1];
1259 1267          z2 = u * u;
1260 1268          k = __HI(u) & 0x7fffffff;
1261 1269          t3 = T2[j2] + T1[n2];
     1270 +
1262 1271          if ((k >> 20) < 0x3ec) {        /* |u|<2**-19 */
1263 1272                  t2 = t4 + u * ((two + z2 * A1) + (z2 * z2) * (A2 + z2 * A3));
1264 1273          } else {
1265 1274                  t5 = t4 + u * (z2 * A1 + (z2 * z2) * (A2 + z2 * A3));
1266 1275                  u2 = u + u;
1267      -                v = (double) ((int) (u2 * t24)) * p24;
     1276 +                v = (double)((int)(u2 * t24)) * p24;
1268 1277                  t2 = t5 + r * ((two * t2 - v * t1) - v * (y - (t1 - z)));
1269 1278                  t3 += v;
1270 1279          }
1271      -        ss_h = (double) ((float) (t2 + t3));
     1280 +
     1281 +        ss_h = (double)((float)(t2 + t3));
1272 1282          ss_l = t2 - (ss_h - t3);
1273 1283  
1274 1284          /*
1275 1285           * compute ww = (x-.5)*(log(x)-1) + .5*(log(2pi)-1) + 1/x*(P(1/x^2)))
1276 1286           * where ss = log(x) - 1 in already in extra precision
1277 1287           */
1278 1288          z = one / x;
1279 1289          r = x - half;
1280      -        r_h = (double) ((float) r);
     1290 +        r_h = (double)((float)r);
1281 1291          w_h = r_h * ss_h + hln2pi_h;
1282 1292          z2 = z * z;
1283 1293          w = (r - r_h) * ss_h + r * ss_l;
1284 1294          z4 = z2 * z2;
1285 1295          t1 = z2 * (GP1 + z4 * (GP3 + z4 * (GP5 + z4 * GP7)));
1286 1296          t2 = z4 * (GP2 + z4 * (GP4 + z4 * GP6));
1287 1297          t1 += t2;
1288 1298          w += hln2pi_l;
1289 1299          w_l = z * (GP0 + t1) + w;
1290      -        k = (int) ((w_h + w_l) * invln2_32 + half);
     1300 +        k = (int)((w_h + w_l) * invln2_32 + half);
1291 1301  
1292 1302          /* compute the exponential of w_h+w_l */
1293 1303          j = k & 0x1f;
1294 1304          *m = (k >> 5);
1295      -        t3 = (double) k;
     1305 +        t3 = (double)k;
1296 1306  
1297 1307          /* perform w - k*ln2_32 (represent as w_h - w_l) */
1298 1308          t1 = w_h - t3 * ln2_32hi;
1299 1309          t2 = t3 * ln2_32lo;
1300 1310          w = w_l - t2;
1301 1311          w_h = t1 + w_l;
1302 1312          w_l = t2 - (w_l - (w_h - t1));
1303 1313  
1304 1314          /* compute exp(w_h+w_l) */
1305 1315          z = w_h - w_l;
1306 1316          z2 = z * z;
1307 1317          t1 = z2 * (Et1 + z2 * (Et3 + z2 * Et5));
1308 1318          t2 = z2 * (Et2 + z2 * Et4);
1309 1319          t3 = w_h - (w_l - (t1 + z * t2));
1310 1320          zz.l = S_trail[j] * (one + t3) + S[j] * t3;
1311 1321          zz.h = S[j];
1312 1322          return (zz);
1313 1323  }
1314 1324  
1315      -/* INDENT OFF */
     1325 +
1316 1326  /*
1317 1327   * kpsin(x)= sin(pi*x)/pi
1318 1328   *                 3        5        7        9        11        13        15
1319 1329   *      = x+ks[0]*x +ks[1]*x +ks[2]*x +ks[3]*x +ks[4]*x  +ks[5]*x  +ks[6]*x
1320 1330   */
1321 1331  static const double ks[] = {
1322 1332          -1.64493406684822640606569,
1323 1333          +8.11742425283341655883668741874008920850698590621e-0001,
1324 1334          -1.90751824120862873825597279118304943994042258291e-0001,
1325 1335          +2.61478477632554278317289628332654539353521911570e-0002,
1326 1336          -2.34607978510202710377617190278735525354347705866e-0003,
1327 1337          +1.48413292290051695897242899977121846763824221705e-0004,
1328 1338          -6.87730769637543488108688726777687262485357072242e-0006,
1329 1339  };
1330      -/* INDENT ON */
1331 1340  
1332 1341  /* assume x is not tiny and positive */
1333 1342  static struct Double
1334      -kpsin(double x) {
     1343 +kpsin(double x)
     1344 +{
1335 1345          double z, t1, t2, t3, t4;
1336 1346          struct Double xx;
1337 1347  
1338 1348          z = x * x;
1339 1349          xx.h = x;
1340 1350          t1 = z * x;
1341 1351          t2 = z * z;
1342 1352          t4 = t1 * ks[0];
1343      -        t3 = (t1 * z) * ((ks[1] + z * ks[2] + t2 * ks[3]) + (z * t2) *
1344      -                (ks[4] + z * ks[5] + t2 * ks[6]));
     1353 +        t3 = (t1 * z) * ((ks[1] + z * ks[2] + t2 * ks[3]) + (z * t2) * (ks[4] +
     1354 +            z * ks[5] + t2 * ks[6]));
1345 1355          xx.l = t4 + t3;
1346 1356          return (xx);
1347 1357  }
1348 1358  
1349      -/* INDENT OFF */
     1359 +
1350 1360  /*
1351 1361   * kpcos(x)= cos(pi*x)/pi
1352 1362   *                     2        4        6        8        10        12
1353 1363   *      = 1/pi +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x +kc[4]*x  +kc[5]*x
1354 1364   */
1355 1365  
1356 1366  static const double one_pi_h = 0.318309886183790635705292970,
1357      -                one_pi_l = 3.583247455607534006714276420e-17;
     1367 +        one_pi_l = 3.583247455607534006714276420e-17;
1358 1368  static const double npi_2_h = -1.5625,
1359      -                npi_2_l = -0.00829632679489661923132169163975055099555883223;
     1369 +        npi_2_l = -0.00829632679489661923132169163975055099555883223;
1360 1370  static const double kc[] = {
1361 1371          -1.57079632679489661923132169163975055099555883223e+0000,
1362 1372          +1.29192819501230224953283586722575766189551966008e+0000,
1363 1373          -4.25027339940149518500158850753393173519732149213e-0001,
1364 1374          +7.49080625187015312373925142219429422375556727752e-0002,
1365 1375          -8.21442040906099210866977352284054849051348692715e-0003,
1366 1376          +6.10411356829515414575566564733632532333904115968e-0004,
1367 1377  };
1368      -/* INDENT ON */
1369 1378  
1370 1379  /* assume x is not tiny and positive */
1371 1380  static struct Double
1372      -kpcos(double x) {
     1381 +kpcos(double x)
     1382 +{
1373 1383          double z, t1, t2, t3, t4, x4, x8;
1374 1384          struct Double xx;
1375 1385  
1376 1386          z = x * x;
1377 1387          xx.h = one_pi_h;
1378      -        t1 = (double) ((float) x);
     1388 +        t1 = (double)((float)x);
1379 1389          x4 = z * z;
1380 1390          t2 = npi_2_l * z + npi_2_h * (x + t1) * (x - t1);
1381      -        t3 = one_pi_l + x4 * ((kc[1] + z * kc[2]) + x4 * (kc[3] + z *
1382      -                kc[4] + x4 * kc[5]));
1383      -        t4 = t1 * t1;   /* 48 bits mantissa */
     1391 +        t3 = one_pi_l + x4 * ((kc[1] + z * kc[2]) + x4 * (kc[3] + z * kc[4] +
     1392 +            x4 * kc[5]));
     1393 +        t4 = t1 * t1;           /* 48 bits mantissa */
1384 1394          x8 = t2 + t3;
1385      -        t4 *= npi_2_h;  /* npi_2_h is 5 bits const. The product is exact */
1386      -        xx.l = x8 + t4; /* that will minimized the rounding error in xx.l */
     1395 +        t4 *= npi_2_h;    /* npi_2_h is 5 bits const. The product is exact */
     1396 +        xx.l = x8 + t4;   /* that will minimized the rounding error in xx.l */
1387 1397          return (xx);
1388 1398  }
1389 1399  
1390      -/* INDENT OFF */
1391 1400  static const double
1392      -        /* 0.134861805732790769689793935774652917006 */
1393      -        t0z1   =  0.1348618057327907737708,
     1401 +/* 0.134861805732790769689793935774652917006 */
     1402 +        t0z1 = 0.1348618057327907737708,
1394 1403          t0z1_l = -4.0810077708578299022531e-18,
1395      -        /* 0.461632144968362341262659542325721328468 */
1396      -        t0z2   =  0.4616321449683623567850,
     1404 +/* 0.461632144968362341262659542325721328468 */
     1405 +        t0z2 = 0.4616321449683623567850,
1397 1406          t0z2_l = -1.5522348162858676890521e-17,
1398      -        /* 0.819773101100500601787868704921606996312 */
1399      -        t0z3   =  0.8197731011005006118708,
     1407 +/* 0.819773101100500601787868704921606996312 */
     1408 +        t0z3 = 0.8197731011005006118708,
1400 1409          t0z3_l = -1.0082945122487103498325e-17;
1401      -        /* 1.134861805732790769689793935774652917006 */
1402      -/* INDENT ON */
     1410 +
     1411 +/*
     1412 + * 1.134861805732790769689793935774652917006
     1413 + */
1403 1414  
1404 1415  /* gamma(x+i) for 0 <= x < 1  */
1405 1416  static struct Double
1406      -gam_n(int i, double x) {
1407      -        struct Double rr = {0.0L, 0.0L}, yy;
     1417 +gam_n(int i, double x)
     1418 +{
     1419 +        struct Double rr = { 0.0L, 0.0L }, yy;
1408 1420          double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl;
1409 1421  
1410 1422          /* compute yy = gamma(x+1) */
1411 1423          if (x > 0.2845) {
1412 1424                  if (x > 0.6374) {
1413 1425                          r1 = x - t0z3;
1414      -                        r2 = (double) ((float) (r1 - t0z3_l));
     1426 +                        r2 = (double)((float)(r1 - t0z3_l));
1415 1427                          t2 = r1 - r2;
1416 1428                          yy = GT3(r2, t2 - t0z3_l);
1417 1429                  } else {
1418 1430                          r1 = x - t0z2;
1419      -                        r2 = (double) ((float) (r1 - t0z2_l));
     1431 +                        r2 = (double)((float)(r1 - t0z2_l));
1420 1432                          t2 = r1 - r2;
1421 1433                          yy = GT2(r2, t2 - t0z2_l);
1422 1434                  }
1423 1435          } else {
1424 1436                  r1 = x - t0z1;
1425      -                r2 = (double) ((float) (r1 - t0z1_l));
     1437 +                r2 = (double)((float)(r1 - t0z1_l));
1426 1438                  t2 = r1 - r2;
1427 1439                  yy = GT1(r2, t2 - t0z1_l);
1428 1440          }
1429 1441  
1430 1442          /* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */
1431 1443          switch (i) {
1432      -        case 0:         /* yy/x */
     1444 +        case 0:                                 /* yy/x */
1433 1445                  r1 = one / x;
1434      -                xh = (double) ((float) x);      /* x is not tiny */
1435      -                rr.h = (double) ((float) ((yy.h + yy.l) * r1));
1436      -                rr.l = r1 * (yy.h - rr.h * xh) -
1437      -                        ((r1 * rr.h) * (x - xh) - r1 * yy.l);
     1446 +                xh = (double)((float)x);        /* x is not tiny */
     1447 +                rr.h = (double)((float)((yy.h + yy.l) * r1));
     1448 +                rr.l = r1 * (yy.h - rr.h * xh) - ((r1 * rr.h) * (x - xh) - r1 *
     1449 +                    yy.l);
1438 1450                  break;
1439      -        case 1:         /* yy */
     1451 +        case 1:                         /* yy */
1440 1452                  rr.h = yy.h;
1441 1453                  rr.l = yy.l;
1442 1454                  break;
1443      -        case 2:         /* (x+1)*yy */
1444      -                z = x + one;    /* may not be exact */
1445      -                zh = (double) ((float) z);
     1455 +        case 2:                         /* (x+1)*yy */
     1456 +                z = x + one;            /* may not be exact */
     1457 +                zh = (double)((float)z);
1446 1458                  rr.h = zh * yy.h;
1447 1459                  rr.l = z * yy.l + (x - (zh - one)) * yy.h;
1448 1460                  break;
1449      -        case 3:         /* (x+2)*(x+1)*yy */
     1461 +        case 3:                         /* (x+2)*(x+1)*yy */
1450 1462                  z1 = x + one;
1451 1463                  z2 = x + 2.0;
1452 1464                  z = z1 * z2;
1453      -                xh = (double) ((float) z);
1454      -                zh = (double) ((float) z1);
     1465 +                xh = (double)((float)z);
     1466 +                zh = (double)((float)z1);
1455 1467                  xl = (x - (zh - one)) * (z2 + zh) - (xh - zh * (zh + one));
1456 1468                  rr.h = xh * yy.h;
1457 1469                  rr.l = z * yy.l + xl * yy.h;
1458 1470                  break;
1459 1471  
1460      -        case 4:         /* (x+1)*(x+3)*(x+2)*yy */
     1472 +        case 4:                         /* (x+1)*(x+3)*(x+2)*yy */
1461 1473                  z1 = x + 2.0;
1462 1474                  z2 = (x + one) * (x + 3.0);
1463 1475                  zh = z1;
1464 1476                  __LO(zh) = 0;
1465 1477                  __HI(zh) &= 0xfffffff8; /* zh 18 bits mantissa */
1466 1478                  zl = x - (zh - 2.0);
1467 1479                  z = z1 * z2;
1468      -                xh = (double) ((float) z);
     1480 +                xh = (double)((float)z);
1469 1481                  xl = zl * (z2 + zh * (z1 + zh)) - (xh - zh * (zh * zh - one));
1470 1482                  rr.h = xh * yy.h;
1471 1483                  rr.l = z * yy.l + xl * yy.h;
1472 1484                  break;
1473      -        case 5:         /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */
     1485 +        case 5:                         /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */
1474 1486                  z1 = x + 2.0;
1475 1487                  z2 = x + 3.0;
1476 1488                  z = z1 * z2;
1477      -                zh = (double) ((float) z1);
1478      -                yh = (double) ((float) z);
     1489 +                zh = (double)((float)z1);
     1490 +                yh = (double)((float)z);
1479 1491                  yl = (x - (zh - 2.0)) * (z2 + zh) - (yh - zh * (zh + one));
1480 1492                  z2 = z - 2.0;
1481 1493                  z *= z2;
1482      -                xh = (double) ((float) z);
     1494 +                xh = (double)((float)z);
1483 1495                  xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0));
1484 1496                  rr.h = xh * yy.h;
1485 1497                  rr.l = z * yy.l + xl * yy.h;
1486 1498                  break;
1487      -        case 6:         /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
     1499 +        case 6:                         /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
1488 1500                  z1 = x + 2.0;
1489 1501                  z2 = x + 3.0;
1490 1502                  z = z1 * z2;
1491      -                zh = (double) ((float) z1);
1492      -                yh = (double) ((float) z);
     1503 +                zh = (double)((float)z1);
     1504 +                yh = (double)((float)z);
1493 1505                  z1 = x - (zh - 2.0);
1494 1506                  yl = z1 * (z2 + zh) - (yh - zh * (zh + one));
1495 1507                  z2 = z - 2.0;
1496 1508                  x5 = x + 5.0;
1497 1509                  z *= z2;
1498      -                xh = (double) ((float) z);
     1510 +                xh = (double)((float)z);
1499 1511                  zh += 3.0;
1500 1512                  xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0));
1501      -                                                /* xh+xl=(x+1)*...*(x+4) */
1502      -                /* wh+wl=(x+5)*yy */
1503      -                wh = (double) ((float) (x5 * (yy.h + yy.l)));
     1513 +
     1514 +                /*
     1515 +                 * xh+xl=(x+1)*...*(x+4)
     1516 +                 * wh+wl=(x+5)*yy
     1517 +                 */
     1518 +                wh = (double)((float)(x5 * (yy.h + yy.l)));
1504 1519                  wl = (z1 * yy.h + x5 * yy.l) - (wh - zh * yy.h);
1505 1520                  rr.h = wh * xh;
1506 1521                  rr.l = z * wl + xl * wh;
1507 1522                  break;
1508      -        case 7:         /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
     1523 +        case 7:                 /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
1509 1524                  z1 = x + 3.0;
1510 1525                  z2 = x + 4.0;
1511 1526                  z = z2 * z1;
1512      -                zh = (double) ((float) z1);
1513      -                yh = (double) ((float) z);      /* yh+yl = (x+3)(x+4) */
     1527 +                zh = (double)((float)z1);
     1528 +                yh = (double)((float)z);        /* yh+yl = (x+3)(x+4) */
1514 1529                  yl = (x - (zh - 3.0)) * (z2 + zh) - (yh - (zh * (zh + one)));
1515 1530                  z1 = x + 6.0;
1516      -                z2 = z - 2.0;   /* z2 = (x+2)*(x+5) */
     1531 +                z2 = z - 2.0;                   /* z2 = (x+2)*(x+5) */
1517 1532                  z *= z2;
1518      -                xh = (double) ((float) z);
     1533 +                xh = (double)((float)z);
1519 1534                  xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0));
1520      -                                                /* xh+xl=(x+2)*...*(x+5) */
1521      -                /* wh+wl=(x+1)(x+6)*yy */
1522      -                z2 -= 4.0;      /* z2 = (x+1)(x+6) */
1523      -                wh = (double) ((float) (z2 * (yy.h + yy.l)));
     1535 +
     1536 +                /*
     1537 +                 * xh+xl=(x+2)*...*(x+5)
     1538 +                 * wh+wl=(x+1)(x+6)*yy
     1539 +                 */
     1540 +                z2 -= 4.0;              /* z2 = (x+1)(x+6) */
     1541 +                wh = (double)((float)(z2 * (yy.h + yy.l)));
1524 1542                  wl = (z2 * yy.l + yl * yy.h) - (wh - (yh - 6.0) * yy.h);
1525 1543                  rr.h = wh * xh;
1526 1544                  rr.l = z * wl + xl * wh;
1527 1545          }
     1546 +
1528 1547          return (rr);
1529 1548  }
1530 1549  
1531 1550  double
1532      -tgamma(double x) {
     1551 +tgamma(double x)
     1552 +{
1533 1553          struct Double ss, ww;
1534 1554          double t, t1, t2, t3, t4, t5, w, y, z, z1, z2, z3, z5;
1535 1555          int i, j, k, m, ix, hx, xk;
1536 1556          unsigned lx;
1537 1557  
1538 1558          hx = __HI(x);
1539 1559          lx = __LO(x);
1540 1560          ix = hx & 0x7fffffff;
1541 1561          y = x;
1542 1562  
1543 1563          if (ix < 0x3ca00000)
1544 1564                  return (one / x);       /* |x| < 2**-53 */
     1565 +
1545 1566          if (ix >= 0x7ff00000)
1546      -                        /* +Inf -> +Inf, -Inf or NaN -> NaN */
1547      -                return (x * ((hx < 0)? 0.0 : x));
1548      -        if (hx > 0x406573fa ||  /* x > 171.62... overflow to +inf */
     1567 +                /* +Inf -> +Inf, -Inf or NaN -> NaN */
     1568 +                return (x * ((hx < 0) ? 0.0 : x));
     1569 +
     1570 +        if (hx > 0x406573fa ||          /* x > 171.62... overflow to +inf */
1549 1571              (hx == 0x406573fa && lx > 0xE561F647)) {
1550 1572                  z = x / tiny;
1551 1573                  return (z * z);
1552 1574          }
1553      -        if (hx >= 0x40200000) { /* x >= 8 */
     1575 +
     1576 +        if (hx >= 0x40200000) {         /* x >= 8 */
1554 1577                  ww = large_gam(x, &m);
1555 1578                  w = ww.h + ww.l;
1556 1579                  __HI(w) += m << 20;
1557 1580                  return (w);
1558 1581          }
1559      -        if (hx > 0) {           /* 0 < x < 8 */
1560      -                i = (int) x;
1561      -                ww = gam_n(i, x - (double) i);
     1582 +
     1583 +        if (hx > 0) {                   /* 0 < x < 8 */
     1584 +                i = (int)x;
     1585 +                ww = gam_n(i, x - (double)i);
1562 1586                  return (ww.h + ww.l);
1563 1587          }
1564 1588  
1565      -        /* negative x */
1566      -        /* INDENT OFF */
     1589 +        /*
     1590 +         * negative x
     1591 +         */
     1592 +
1567 1593          /*
1568 1594           * compute: xk =
1569 1595           *      -2 ... x is an even int (-inf is even)
1570 1596           *      -1 ... x is an odd int
1571 1597           *      +0 ... x is not an int but chopped to an even int
1572 1598           *      +1 ... x is not an int but chopped to an odd int
1573 1599           */
1574      -        /* INDENT ON */
1575 1600          xk = 0;
     1601 +
1576 1602          if (ix >= 0x43300000) {
1577 1603                  if (ix >= 0x43400000)
1578 1604                          xk = -2;
1579 1605                  else
1580 1606                          xk = -2 + (lx & 1);
1581 1607          } else if (ix >= 0x3ff00000) {
1582 1608                  k = (ix >> 20) - 0x3ff;
     1609 +
1583 1610                  if (k > 20) {
1584 1611                          j = lx >> (52 - k);
     1612 +
1585 1613                          if ((j << (52 - k)) == lx)
1586 1614                                  xk = -2 + (j & 1);
1587 1615                          else
1588 1616                                  xk = j & 1;
1589 1617                  } else {
1590 1618                          j = ix >> (20 - k);
     1619 +
1591 1620                          if ((j << (20 - k)) == ix && lx == 0)
1592 1621                                  xk = -2 + (j & 1);
1593 1622                          else
1594 1623                                  xk = j & 1;
1595 1624                  }
1596 1625          }
     1626 +
1597 1627          if (xk < 0)
1598 1628                  /* ideally gamma(-n)= (-1)**(n+1) * inf, but c99 expect NaN */
1599      -                return ((x - x) / (x - x));             /* 0/0 = NaN */
1600      -
     1629 +                return ((x - x) / (x - x));     /* 0/0 = NaN */
1601 1630  
1602 1631          /* negative underflow thresold */
1603 1632          if (ix > 0x4066e000 || (ix == 0x4066e000 && lx > 11)) {
1604 1633                  /* x < -183.0 - 11ulp */
1605 1634                  z = tiny / x;
     1635 +
1606 1636                  if (xk == 1)
1607 1637                          z = -z;
     1638 +
1608 1639                  return (z * tiny);
1609 1640          }
1610 1641  
1611 1642          /* now compute gamma(x) by  -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x */
1612 1643  
1613 1644          /*
1614 1645           * First compute ss = -sin(pi*y)/pi , so that
1615 1646           * gamma(x) = 1/(ss*gamma(1+y))
1616 1647           */
1617 1648          y = -x;
1618      -        j = (int) y;
1619      -        z = y - (double) j;
1620      -        if (z > 0.3183098861837906715377675)
     1649 +        j = (int)y;
     1650 +        z = y - (double)j;
     1651 +
     1652 +        if (z > 0.3183098861837906715377675) {
1621 1653                  if (z > 0.6816901138162093284622325)
1622 1654                          ss = kpsin(one - z);
1623 1655                  else
1624 1656                          ss = kpcos(0.5 - z);
1625      -        else
     1657 +        } else {
1626 1658                  ss = kpsin(z);
     1659 +        }
     1660 +
1627 1661          if (xk == 0) {
1628 1662                  ss.h = -ss.h;
1629 1663                  ss.l = -ss.l;
1630 1664          }
1631 1665  
1632 1666          /* Then compute ww = gamma(1+y), note that result scale to 2**m */
1633 1667          m = 0;
     1668 +
1634 1669          if (j < 7) {
1635 1670                  ww = gam_n(j + 1, z);
1636 1671          } else {
1637 1672                  w = y + one;
     1673 +
1638 1674                  if ((lx & 1) == 0) {    /* y+1 exact (note that y<184) */
1639 1675                          ww = large_gam(w, &m);
1640 1676                  } else {
1641 1677                          t = w - one;
     1678 +
1642 1679                          if (t == y) {   /* y+one exact */
1643 1680                                  ww = large_gam(w, &m);
1644 1681                          } else {        /* use y*gamma(y) */
1645 1682                                  if (j == 7)
1646 1683                                          ww = gam_n(j, z);
1647 1684                                  else
1648 1685                                          ww = large_gam(y, &m);
     1686 +
1649 1687                                  t4 = ww.h + ww.l;
1650      -                                t1 = (double) ((float) y);
1651      -                                t2 = (double) ((float) t4);
1652      -                                                /* t4 will not be too large */
     1688 +                                t1 = (double)((float)y);
     1689 +                                t2 = (double)((float)t4);
     1690 +                                /* t4 will not be too large */
1653 1691                                  ww.l = y * (ww.l - (t2 - ww.h)) + (y - t1) * t2;
1654 1692                                  ww.h = t1 * t2;
1655 1693                          }
1656 1694                  }
1657 1695          }
1658 1696  
1659 1697          /* compute 1/(ss*ww) */
1660 1698          t3 = ss.h + ss.l;
1661 1699          t4 = ww.h + ww.l;
1662      -        t1 = (double) ((float) t3);
1663      -        t2 = (double) ((float) t4);
     1700 +        t1 = (double)((float)t3);
     1701 +        t2 = (double)((float)t4);
1664 1702          z1 = ss.l - (t1 - ss.h);        /* (t1,z1) = ss */
1665 1703          z2 = ww.l - (t2 - ww.h);        /* (t2,z2) = ww */
1666 1704          t3 = t3 * t4;                   /* t3 = ss*ww */
1667 1705          z3 = one / t3;                  /* z3 = 1/(ss*ww) */
1668 1706          t5 = t1 * t2;
1669 1707          z5 = z1 * t4 + t1 * z2;         /* (t5,z5) = ss*ww */
1670      -        t1 = (double) ((float) t3);     /* (t1,z1) = ss*ww */
     1708 +        t1 = (double)((float)t3);       /* (t1,z1) = ss*ww */
1671 1709          z1 = z5 - (t1 - t5);
1672      -        t2 = (double) ((float) z3);     /* leading 1/(ss*ww) */
     1710 +        t2 = (double)((float)z3);       /* leading 1/(ss*ww) */
1673 1711          z2 = z3 * (t2 * z1 - (one - t2 * t1));
1674 1712          z = t2 - z2;
1675 1713  
1676 1714          /* check whether z*2**-m underflow */
1677 1715          if (m != 0) {
1678 1716                  hx = __HI(z);
1679 1717                  i = hx & 0x80000000;
1680 1718                  ix = hx ^ i;
1681 1719                  j = ix >> 20;
     1720 +
1682 1721                  if (j > m) {
1683 1722                          ix -= m << 20;
1684 1723                          __HI(z) = ix ^ i;
1685 1724                  } else if ((m - j) > 52) {
1686 1725                          /* underflow */
1687 1726                          if (xk == 0)
1688 1727                                  z = -tiny * tiny;
1689 1728                          else
1690 1729                                  z = tiny * tiny;
1691 1730                  } else {
1692 1731                          /* subnormal */
1693 1732                          m -= 60;
1694 1733                          t = one;
1695 1734                          __HI(t) -= 60 << 20;
1696 1735                          ix -= m << 20;
1697 1736                          __HI(z) = ix ^ i;
1698 1737                          z *= t;
1699 1738                  }
1700 1739          }
     1740 +
1701 1741          return (z);
1702 1742  }
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX