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

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/m9x/tgammal.c
          +++ new/usr/src/lib/libm/common/m9x/tgammal.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 __tgammal = tgammal
  31   32  
  32   33  #include "libm.h"
  33   34  #include <sys/isa_defs.h>
  34   35  
  35   36  #if defined(_BIG_ENDIAN)
  36      -#define H0_WORD(x)      ((unsigned *) &x)[0]
  37      -#define H3_WORD(x)      ((unsigned *) &x)[3]
  38      -#define CHOPPED(x)      (long double) ((double) (x))
       37 +#define H0_WORD(x)              ((unsigned *)&x)[0]
       38 +#define H3_WORD(x)              ((unsigned *)&x)[3]
       39 +#define CHOPPED(x)              (long double)((double)(x))
  39   40  #else
  40      -#define H0_WORD(x)      ((((int *) &x)[2] << 16) | \
  41      -                        (0x0000ffff & (((unsigned *) &x)[1] >> 15)))
  42      -#define H3_WORD(x)      ((unsigned *) &x)[0]
  43      -#define CHOPPED(x)      (long double) ((float) (x))
       41 +#define H0_WORD(x)              ((((int *)&x)[2] << 16) | (0x0000ffff & \
       42 +        (((unsigned *)&x)[1] >> 15)))
       43 +#define H3_WORD(x)              ((unsigned *)&x)[0]
       44 +#define CHOPPED(x)              (long double)((float)(x))
  44   45  #endif
  45   46  
  46   47  struct LDouble {
  47   48          long double h, l;
  48   49  };
  49   50  
  50      -/* INDENT OFF */
  51      -/* Primary interval GTi() */
       51 +/*
       52 + * Primary interval GTi()
       53 + */
  52   54  static const long double P1[] = {
  53   55          +0.709086836199777919037185741507610124611513720557L,
  54   56          +4.45754781206489035827915969367354835667391606951e-0001L,
  55   57          +3.21049298735832382311662273882632210062918153852e-0002L,
  56   58          -5.71296796342106617651765245858289197369688864350e-0003L,
  57   59          +6.04666892891998977081619174969855831606965352773e-0003L,
  58   60          +8.99106186996888711939627812174765258822658645168e-0004L,
  59   61          -6.96496846144407741431207008527018441810175568949e-0005L,
  60   62          +1.52597046118984020814225409300131445070213882429e-0005L,
  61   63          +5.68521076168495673844711465407432189190681541547e-0007L,
  62   64          +3.30749673519634895220582062520286565610418952979e-0008L,
  63   65  };
       66 +
  64   67  static const long double Q1[] = {
  65      -        +1.0+0000L,
       68 +        +1.0 + 0000L,
  66   69          +1.35806511721671070408570853537257079579490650668e+0000L,
  67   70          +2.97567810153429553405327140096063086994072952961e-0001L,
  68   71          -1.52956835982588571502954372821681851681118097870e-0001L,
  69   72          -2.88248519561420109768781615289082053597954521218e-0002L,
  70   73          +1.03475311719937405219789948456313936302378395955e-0002L,
  71   74          +4.12310203243891222368965360124391297374822742313e-0004L,
  72   75          -3.12653708152290867248931925120380729518332507388e-0004L,
  73   76          +2.36672170850409745237358105667757760527014332458e-0005L,
  74   77  };
       78 +
  75   79  static const long double P2[] = {
  76   80          +0.428486815855585429730209907810650135255270600668084114L,
  77   81          +2.62768479103809762805691743305424077975230551176e-0001L,
  78   82          +3.81187532685392297608310837995193946591425896150e-0002L,
  79   83          +3.00063075891811043820666846129131255948527925381e-0003L,
  80   84          +2.47315407812279164228398470797498649142513408654e-0003L,
  81   85          +3.62838199917848372586173483147214880464782938664e-0004L,
  82   86          +3.43991105975492623982725644046473030098172692423e-0006L,
  83   87          +4.56902151569603272237014240794257659159045432895e-0006L,
  84   88          +2.13734755837595695602045100675540011352948958453e-0007L,
  85   89          +9.74123440547918230781670266967882492234877125358e-0009L,
  86   90  };
       91 +
  87   92  static const long double Q2[] = {
  88   93          +1.0L,
  89   94          +9.18284118632506842664645516830761489700556179701e-0001L,
  90   95          -6.41430858837830766045202076965923776189154874947e-0003L,
  91   96          -1.24400885809771073213345747437964149775410921376e-0001L,
  92   97          +4.69803798146251757538856567522481979624746875964e-0003L,
  93   98          +7.18309447069495315914284705109868696262662082731e-0003L,
  94   99          -8.75812626987894695112722600697653425786166399105e-0004L,
  95  100          -1.23539972377769277995959339188431498626674835169e-0004L,
  96  101          +3.10019017590151598732360097849672925448587547746e-0005L,
  97  102          -1.77260223349332617658921874288026777465782364070e-0006L,
  98  103  };
      104 +
  99  105  static const long double P3[] = {
 100  106          +0.3824094797345675048502747661075355640070439388902L,
 101  107          +3.42198093076618495415854906335908427159833377774e-0001L,
 102  108          +9.63828189500585568303961406863153237440702754858e-0002L,
 103  109          +8.76069421042696384852462044188520252156846768667e-0003L,
 104  110          +1.86477890389161491224872014149309015261897537488e-0003L,
 105  111          +8.16871354540309895879974742853701311541286944191e-0004L,
 106  112          +6.83783483674600322518695090864659381650125625216e-0005L,
 107  113          -1.10168269719261574708565935172719209272190828456e-0006L,
 108  114          +9.66243228508380420159234853278906717065629721016e-0007L,
 109  115          +2.31858885579177250541163820671121664974334728142e-0008L,
 110  116  };
      117 +
 111  118  static const long double Q3[] = {
 112      -        +1.0L,
 113      -        +8.25479821168813634632437430090376252512793067339e-0001L,
      119 +        +1.0L, +8.25479821168813634632437430090376252512793067339e-0001L,
 114  120          -1.62251363073937769739639623669295110346015576320e-0002L,
 115  121          -1.10621286905916732758745130629426559691187579852e-0001L,
 116  122          +3.48309693970985612644446415789230015515365291459e-0003L,
 117  123          +6.73553737487488333032431261131289672347043401328e-0003L,
 118  124          -7.63222008393372630162743587811004613050245128051e-0004L,
 119  125          -1.35792670669190631476784768961953711773073251336e-0004L,
 120  126          +3.19610150954223587006220730065608156460205690618e-0005L,
 121  127          -1.82096553862822346610109522015129585693354348322e-0006L,
 122  128  };
 123  129  
 124  130  static const long double
 125  131  #if defined(__x86)
 126      -GZ1_h   =  0.938204627909682449364570100414084663498215377L,
 127      -GZ1_l   =  4.518346116624229420055327632718530617227944106e-20L,
 128      -GZ2_h   =  0.885603194410888700264725126309883762587560340L,
 129      -GZ2_l   =  1.409077427270497062039119290776508217077297169e-20L,
 130      -GZ3_h   =  0.936781411463652321613537060640553022494714241L,
 131      -GZ3_l   =  5.309836440284827247897772963887219035221996813e-21L,
      132 +        GZ1_h = 0.938204627909682449364570100414084663498215377L,
      133 +        GZ1_l = 4.518346116624229420055327632718530617227944106e-20L,
      134 +        GZ2_h = 0.885603194410888700264725126309883762587560340L,
      135 +        GZ2_l = 1.409077427270497062039119290776508217077297169e-20L,
      136 +        GZ3_h = 0.936781411463652321613537060640553022494714241L,
      137 +        GZ3_l = 5.309836440284827247897772963887219035221996813e-21L,
 132  138  #else
 133      -GZ1_h   =  0.938204627909682449409753561580326910854647031L,
 134      -GZ1_l   =  4.684412162199460089642452580902345976446297037e-35L,
 135      -GZ2_h   =  0.885603194410888700278815900582588658192658794L,
 136      -GZ2_l   =  7.501529273890253789219935569758713534641074860e-35L,
 137      -GZ3_h   =  0.936781411463652321618846897080837818855399840L,
 138      -GZ3_l   =  3.088721217404784363585591914529361687403776917e-35L,
      139 +        GZ1_h = 0.938204627909682449409753561580326910854647031L,
      140 +        GZ1_l = 4.684412162199460089642452580902345976446297037e-35L,
      141 +        GZ2_h = 0.885603194410888700278815900582588658192658794L,
      142 +        GZ2_l = 7.501529273890253789219935569758713534641074860e-35L,
      143 +        GZ3_h = 0.936781411463652321618846897080837818855399840L,
      144 +        GZ3_l = 3.088721217404784363585591914529361687403776917e-35L,
 139  145  #endif
 140      -TZ1     = -0.3517214357852935791015625L,
 141      -TZ3     =  0.280530631542205810546875L;
 142      -/* INDENT ON */
      146 +        TZ1 = -0.3517214357852935791015625L,
      147 +        TZ3 = 0.280530631542205810546875L;
      148 +
 143  149  
 144      -/* INDENT OFF */
 145  150  /*
 146  151   * compute gamma(y=yh+yl) for y in GT1 = [1.0000, 1.2845]
 147  152   * ...assume yh got 53 or 24(i386) significant bits
 148  153   */
 149      -/* INDENT ON */
 150  154  static struct LDouble
 151      -GT1(long double yh, long double yl) {
      155 +GT1(long double yh, long double yl)
      156 +{
 152  157          long double t3, t4, y;
 153  158          int i;
 154  159          struct LDouble r;
 155  160  
 156  161          y = yh + yl;
      162 +
 157  163          for (t4 = Q1[8], t3 = P1[8] + y * P1[9], i = 7; i >= 0; i--) {
 158  164                  t4 = t4 * y + Q1[i];
 159  165                  t3 = t3 * y + P1[i];
 160  166          }
      167 +
 161  168          t3 = (y * y) * t3 / t4;
 162  169          t3 += (TZ1 * yl + GZ1_l);
 163  170          t4 = TZ1 * yh;
 164  171          r.h = CHOPPED((t4 + GZ1_h + t3));
 165  172          t3 += (t4 - (r.h - GZ1_h));
 166  173          r.l = t3;
 167  174          return (r);
 168  175  }
 169  176  
 170      -/* INDENT OFF */
      177 +
 171  178  /*
 172  179   * compute gamma(y=yh+yl) for y in GT2 = [1.2844, 1.6374]
 173  180   * ...assume yh got 53 significant bits
 174  181   */
 175      -/* INDENT ON */
 176  182  static struct LDouble
 177      -GT2(long double yh, long double yl) {
      183 +GT2(long double yh, long double yl)
      184 +{
 178  185          long double t3, t4, y;
 179  186          int i;
 180  187          struct LDouble r;
 181  188  
 182  189          y = yh + yl;
      190 +
 183  191          for (t4 = Q2[9], t3 = P2[9], i = 8; i >= 0; i--) {
 184  192                  t4 = t4 * y + Q2[i];
 185  193                  t3 = t3 * y + P2[i];
 186  194          }
      195 +
 187  196          t3 = GZ2_l + (y * y) * t3 / t4;
 188  197          r.h = CHOPPED((GZ2_h + t3));
 189  198          r.l = t3 - (r.h - GZ2_h);
 190  199          return (r);
 191  200  }
 192  201  
 193      -/* INDENT OFF */
      202 +
 194  203  /*
 195  204   * compute gamma(y=yh+yl) for y in GT3 = [1.6373, 2.0000]
 196  205   * ...assume yh got 53 significant bits
 197  206   */
 198      -/* INDENT ON */
 199  207  static struct LDouble
 200      -GT3(long double yh, long double yl) {
      208 +GT3(long double yh, long double yl)
      209 +{
 201  210          long double t3, t4, y;
 202  211          int i;
 203  212          struct LDouble r;
 204  213  
 205  214          y = yh + yl;
      215 +
 206  216          for (t4 = Q3[9], t3 = P3[9], i = 8; i >= 0; i--) {
 207  217                  t4 = t4 * y + Q3[i];
 208  218                  t3 = t3 * y + P3[i];
 209  219          }
      220 +
 210  221          t3 = (y * y) * t3 / t4;
 211  222          t3 += (TZ3 * yl + GZ3_l);
 212  223          t4 = TZ3 * yh;
 213  224          r.h = CHOPPED((t4 + GZ3_h + t3));
 214  225          t3 += (t4 - (r.h - GZ3_h));
 215  226          r.l = t3;
 216  227          return (r);
 217  228  }
 218  229  
 219      -/* INDENT OFF */
 220      -/* Hex value of GP[0] shoule be 3FB55555 55555555 */
      230 +/*
      231 + * Hex value of GP[0] shoule be 3FB55555 55555555
      232 + */
 221  233  static const long double GP[] = {
 222  234          +0.083333333333333333333333333333333172839171301L,
 223  235          -2.77777777777777777777777777492501211999399424104e-0003L,
 224  236          +7.93650793650793650793635650541638236350020883243e-0004L,
 225  237          -5.95238095238095238057299772679324503339241961704e-0004L,
 226  238          +8.41750841750841696138422987977683524926142600321e-0004L,
 227  239          -1.91752691752686682825032547823699662178842123308e-0003L,
 228  240          +6.41025641022403480921891559356473451161279359322e-0003L,
 229  241          -2.95506535798414019189819587455577003732808185071e-0002L,
 230  242          +1.79644367229970031486079180060923073476568732136e-0001L,
 231  243          -1.39243086487274662174562872567057200255649290646e+0000L,
 232  244          +1.34025874044417962188677816477842265259608269775e+0001L,
 233  245          -1.56803713480127469414495545399982508700748274318e+0002L,
 234  246          +2.18739841656201561694927630335099313968924493891e+0003L,
 235  247          -3.55249848644100338419187038090925410976237921269e+0004L,
 236  248          +6.43464880437835286216768959439484376449179576452e+0005L,
 237  249          -1.20459154385577014992600342782821389605893904624e+0007L,
 238  250          +2.09263249637351298563934942349749718491071093210e+0008L,
 239  251          -2.96247483183169219343745316433899599834685703457e+0009L,
 240  252          +2.88984933605896033154727626086506756972327292981e+0010L,
 241      -        -1.40960434146030007732838382416230610302678063984e+0011L,      /* 19 */
      253 +        -1.40960434146030007732838382416230610302678063984e+0011L, /* 19 */
 242  254  };
 243  255  
 244  256  static const long double T3[] = {
 245      -        +0.666666666666666666666666666666666634567834260213L,   /* T3[0] */
 246      -        +0.400000000000000000000000000040853636176634934140L,   /* T3[1] */
 247      -        +0.285714285714285714285696975252753987869020263448L,   /* T3[2] */
 248      -        +0.222222222222222225593221101192317258554772129875L,   /* T3[3] */
 249      -        +0.181818181817850192105847183461778186703779262916L,   /* T3[4] */
 250      -        +0.153846169861348633757101285952333369222567014596L,   /* T3[5] */
 251      -        +0.133033462889260193922261296772841229985047571265L,   /* T3[6] */
      257 +        +0.666666666666666666666666666666666634567834260213L, /* T3[0] */
      258 +        +0.400000000000000000000000000040853636176634934140L, /* T3[1] */
      259 +        +0.285714285714285714285696975252753987869020263448L, /* T3[2] */
      260 +        +0.222222222222222225593221101192317258554772129875L, /* T3[3] */
      261 +        +0.181818181817850192105847183461778186703779262916L, /* T3[4] */
      262 +        +0.153846169861348633757101285952333369222567014596L, /* T3[5] */
      263 +        +0.133033462889260193922261296772841229985047571265L, /* T3[6] */
 252  264  };
 253      -
      265 +/* BEGIN CSTYLED */
 254  266  static const long double c[] = {
 255      -0.0L,
 256      -1.0L,
 257      -2.0L,
 258      -0.5L,
 259      -1.0e-4930L,                                                     /* tiny */
 260      -4.18937683105468750000e-01L,                                    /* hln2pim1_h */
 261      -8.50099203991780329736405617639861397473637783412817152e-07L,   /* hln2pim1_l */
 262      -0.418938533204672741780329736405617639861397473637783412817152L, /* hln2pim1 */
 263      -2.16608493865351192653179168701171875e-02L,                     /* ln2_32hi */
 264      -5.96317165397058692545083025235937919875797669127130e-12L,      /* ln2_32lo */
 265      -46.16624130844682903551758979206054839765267053289554989233L,   /* invln2_32 */
      267 +        0.0L,
      268 +        1.0L,
      269 +        2.0L,
      270 +        0.5L,
      271 +        1.0e-4930L,                  /* tiny */
      272 +        4.18937683105468750000e-01L, /* hln2pim1_h */
      273 +        8.50099203991780329736405617639861397473637783412817152e-07L, /* hln2pim1_l */
      274 +        0.418938533204672741780329736405617639861397473637783412817152L, /* hln2pim1 */
      275 +        2.16608493865351192653179168701171875e-02L, /* ln2_32hi */
      276 +        5.96317165397058692545083025235937919875797669127130e-12L, /* ln2_32lo */
      277 +        46.16624130844682903551758979206054839765267053289554989233L, /* invln2_32 */
 266  278  #if defined(__x86)
 267      -1.7555483429044629170023839037639845628291e+03L,                /* overflow */
      279 +        1.7555483429044629170023839037639845628291e+03L, /* overflow */
 268  280  #else
 269      -1.7555483429044629170038892160702032034177e+03L,                /* overflow */
      281 +        1.7555483429044629170038892160702032034177e+03L, /* overflow */
 270  282  #endif
 271  283  };
      284 +/* END CSTYLED */
 272  285  
 273      -#define zero            c[0]
 274      -#define one             c[1]
 275      -#define two             c[2]
 276      -#define half            c[3]
 277      -#define tiny            c[4]
 278      -#define hln2pim1_h      c[5]
 279      -#define hln2pim1_l      c[6]
 280      -#define hln2pim1        c[7]
 281      -#define ln2_32hi        c[8]
 282      -#define ln2_32lo        c[9]
 283      -#define invln2_32       c[10]
 284      -#define overflow        c[11]
      286 +#define zero                    c[0]
      287 +#define one                     c[1]
      288 +#define two                     c[2]
      289 +#define half                    c[3]
      290 +#define tiny                    c[4]
      291 +#define hln2pim1_h              c[5]
      292 +#define hln2pim1_l              c[6]
      293 +#define hln2pim1                c[7]
      294 +#define ln2_32hi                c[8]
      295 +#define ln2_32lo                c[9]
      296 +#define invln2_32               c[10]
      297 +#define overflow                c[11]
 285  298  
 286  299  /*
 287  300   * |exp(r) - (1+r+Et0*r^2+...+Et10*r^12)| <= 2^(-128.88) for |r|<=ln2/64
 288  301   */
 289  302  static const long double Et[] = {
 290  303          +5.0000000000000000000e-1L,
 291  304          +1.66666666666666666666666666666828835166292152466e-0001L,
 292  305          +4.16666666666666666666666666666693398646592712189e-0002L,
 293  306          +8.33333333333333333333331748774512601775591115951e-0003L,
 294  307          +1.38888888888888888888888845356011511394764753997e-0003L,
↓ open down ↓ 12 lines elided ↑ open up ↑
 307  320   *  log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y,  1<=y<2,
 308  321   *  j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
 309  322   *       T1(n) = T1[2n,2n+1] = n*log(2)-1,
 310  323   *       T2(j) = T2[2j,2j+1] = log(z[j]),
 311  324   *       T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + T3[2]s^7 + ... + T3[6]s^15
 312  325   *  Note
 313  326   *  (1) the leading entries are truncated to 24 binary point.
 314  327   *  (2) Remez error for T3(s) is bounded by 2**(-136.54)
 315  328   */
 316  329  static const long double T1[] = {
 317      --1.000000000000000000000000000000000000000000e+00L,
      330 +        -1.000000000000000000000000000000000000000000e+00L,
 318  331          +0.000000000000000000000000000000000000000000e+00L,
 319      --3.068528175354003906250000000000000000000000e-01L,
 320      --1.904654299957767878541823431924500011926579e-09L,
      332 +        -3.068528175354003906250000000000000000000000e-01L,
      333 +        -1.904654299957767878541823431924500011926579e-09L,
 321  334          +3.862943053245544433593750000000000000000000e-01L,
 322  335          +5.579533617547508924291635313615100141107647e-08L,
 323  336          +1.079441487789154052734375000000000000000000e+00L,
 324  337          +5.389068187551732136437452970422650211661470e-08L,
 325  338          +1.772588670253753662109375000000000000000000e+00L,
 326  339          +5.198602757555955348583270627230200282215294e-08L,
 327  340          +2.465735852718353271484375000000000000000000e+00L,
 328  341          +5.008137327560178560729088284037750352769117e-08L,
 329  342          +3.158883035182952880859375000000000000000000e+00L,
 330  343          +4.817671897564401772874905940845299849351090e-08L,
↓ open down ↓ 206 lines elided ↑ open up ↑
 537  550          +1.68179283050742908606225095246642969e+00L,
 538  551          +1.71861929812247791562934437645631244e+00L,
 539  552          +1.75625216037329948311216061937531314e+00L,
 540  553          +1.79470907500310718642770324212778174e+00L,
 541  554          +1.83400808640934246348708318958828892e+00L,
 542  555          +1.87416763411029990132999894995444645e+00L,
 543  556          +1.91520656139714729387261127029583086e+00L,
 544  557          +1.95714412417540026901832225162687149e+00L,
 545  558  #endif
 546  559  };
      560 +
 547  561  static const long double S_trail[] = {
 548  562  #if defined(__x86)
 549  563          +0.0000000000000000000000000e+00L,
 550  564          +2.6327965667180882569382524e-20L,
 551  565          +8.3765863521895191129661899e-20L,
 552  566          +3.9798705777454504249209575e-20L,
 553  567          +1.0668046596651558640993042e-19L,
 554  568          +1.9376009847285360448117114e-20L,
 555  569          +6.7081819456112953751277576e-21L,
 556  570          +1.9711680502629186462729727e-20L,
↓ open down ↓ 17 lines elided ↑ open up ↑
 574  588          +2.6052966871016580981769728e-20L,
 575  589          +2.6876456344632553875309579e-21L,
 576  590          +1.2861930155613700201703279e-20L,
 577  591          +8.8166633394037485606572294e-20L,
 578  592          +2.9788615389580190940837037e-20L,
 579  593          +5.2352341619805098677422139e-20L,
 580  594          +5.2578463064010463732242363e-20L,
 581  595  #else
 582  596          +0.00000000000000000000000000000000000e+00L,
 583  597          +1.80506787420330954745573333054573786e-35L,
 584      --9.37452029228042742195756741973083214e-35L,
 585      --1.59696844729275877071290963023149997e-35L,
      598 +        -9.37452029228042742195756741973083214e-35L,
      599 +        -1.59696844729275877071290963023149997e-35L,
 586  600          +9.11249341012502297851168610167248666e-35L,
 587      --6.50422820697854828723037477525938871e-35L,
 588      --8.14846884452585113732569176748815532e-35L,
 589      --5.06621457672180031337233074514290335e-35L,
 590      --1.35983097468881697374987563824591912e-35L,
      601 +        -6.50422820697854828723037477525938871e-35L,
      602 +        -8.14846884452585113732569176748815532e-35L,
      603 +        -5.06621457672180031337233074514290335e-35L,
      604 +        -1.35983097468881697374987563824591912e-35L,
 591  605          +9.49742763556319647030771056643324660e-35L,
 592      --3.28317052317699860161506596533391526e-36L,
 593      --5.01723570938719041029018653045842895e-35L,
 594      --2.39147479768910917162283430160264014e-35L,
 595      --8.35057135763390881529889073794408385e-36L,
      606 +        -3.28317052317699860161506596533391526e-36L,
      607 +        -5.01723570938719041029018653045842895e-35L,
      608 +        -2.39147479768910917162283430160264014e-35L,
      609 +        -8.35057135763390881529889073794408385e-36L,
 596  610          +7.03675688907326504242173719067187644e-35L,
 597      --5.18248485306464645753689301856695619e-35L,
      611 +        -5.18248485306464645753689301856695619e-35L,
 598  612          +9.42224254862183206569211673639406488e-35L,
 599      --3.96750082539886230916730613021641828e-35L,
      613 +        -3.96750082539886230916730613021641828e-35L,
 600  614          +7.14352899156330061452327361509276724e-35L,
 601  615          +1.15987125286798512424651783410044433e-35L,
 602  616          +4.69693347835811549530973921320187447e-35L,
 603      --3.38651317599500471079924198499981917e-35L,
 604      --8.58731877429824706886865593510387445e-35L,
 605      --9.60595154874935050318549936224606909e-35L,
      617 +        -3.38651317599500471079924198499981917e-35L,
      618 +        -8.58731877429824706886865593510387445e-35L,
      619 +        -9.60595154874935050318549936224606909e-35L,
 606  620          +9.60973393212801278450755869714178581e-35L,
 607  621          +6.37839792144002843924476144978084855e-35L,
 608  622          +7.79243078569586424945646112516927770e-35L,
 609  623          +7.36133776758845652413193083663393220e-35L,
 610      --6.47299514791334723003521457561217053e-35L,
      624 +        -6.47299514791334723003521457561217053e-35L,
 611  625          +8.58747441795369869427879806229522962e-35L,
 612  626          +2.37181542282517483569165122830269098e-35L,
 613      --3.02689168209611877300459737342190031e-37L,
      627 +        -3.02689168209611877300459737342190031e-37L,
 614  628  #endif
 615  629  };
 616      -/* INDENT ON */
 617  630  
 618      -/* INDENT OFF */
      631 +
 619  632  /*
 620  633   * return tgamma(x) scaled by 2**-m for 8<x<=171.62... using Stirling's formula
 621  634   *     log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x))
 622  635   *                = L1 + L2 + L3,
 623  636   */
 624      -/* INDENT ON */
 625  637  static struct LDouble
 626      -large_gam(long double x, int *m) {
      638 +large_gam(long double x, int *m)
      639 +{
 627  640          long double z, t1, t2, t3, z2, t5, w, y, u, r, v;
 628  641          long double t24 = 16777216.0L, p24 = 1.0L / 16777216.0L;
 629  642          int n2, j2, k, ix, j, i;
 630  643          struct LDouble zz;
 631  644          long double u2, ss_h, ss_l, r_h, w_h, w_l, t4;
 632  645  
 633      -/* INDENT OFF */
      646 +/* BEGIN CSTYLED */
 634  647  /*
 635  648   * compute ss = ss.h+ss.l = log(x)-1 (see tgamma_log.h for details)
 636  649   *
 637  650   *  log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y,  1<=y<2,
 638  651   *  j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
 639  652   *       T1(n) = T1[2n,2n+1] = n*log(2)-1,
 640  653   *       T2(j) = T2[2j,2j+1] = log(z[j]),
 641  654   *       T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + ... + T3[6]s^15
 642  655   *  Note
 643  656   *  (1) the leading entries are truncated to 24 binary point.
↓ open down ↓ 3 lines elided ↑ open up ↑
 647  660   *               T1(n):     |_________|___________________|
 648  661   *                             _______ ______________________
 649  662   *               T2(j):       |_______|______________________|
 650  663   *                                ____ _______________________
 651  664   *               2s:             |____|_______________________|
 652  665   *                                    __________________________
 653  666   *          +    T3(s)-2s:           |__________________________|
 654  667   *                       -------------------------------------------
 655  668   *                          [leading] + [Trailing]
 656  669   */
 657      -        /* INDENT ON */
      670 +        /* END CSTYLED */
 658  671          ix = H0_WORD(x);
 659      -        n2 = (ix >> 16) - 0x3fff;       /* exponent of x, range:3-10 */
 660      -        y = scalbnl(x, -n2);    /* y = scale x to [1,2] */
 661      -        n2 += n2;               /* 2n */
 662      -        j = (ix >> 10) & 0x3f;  /* j */
 663      -        z = 1.0078125L + (long double) j * 0.015625L;   /* z[j]=1+j/64+1/128 */
      672 +        n2 = (ix >> 16) - 0x3fff; /* exponent of x, range:3-10 */
      673 +        y = scalbnl(x, -n2);      /* y = scale x to [1,2] */
      674 +        n2 += n2;                 /* 2n */
      675 +        j = (ix >> 10) & 0x3f;    /* j */
      676 +        z = 1.0078125L + (long double)j * 0.015625L;    /* z[j]=1+j/64+1/128 */
 664  677          j2 = j + j;
 665  678          t1 = y + z;
 666  679          t2 = y - z;
 667  680          r = one / t1;
 668      -        u = r * t2;             /* u = (y-z)/(y+z) */
      681 +        u = r * t2;                     /* u = (y-z)/(y+z) */
 669  682          t1 = CHOPPED(t1);
 670  683          t4 = T2[j2 + 1] + T1[n2 + 1];
 671  684          z2 = u * u;
 672  685          k = H0_WORD(u) & 0x7fffffff;
 673  686          t3 = T2[j2] + T1[n2];
      687 +
 674  688          for (t5 = T3[6], i = 5; i >= 0; i--)
 675  689                  t5 = z2 * t5 + T3[i];
      690 +
 676  691          if ((k >> 16) < 0x3fec) {       /* |u|<2**-19 */
 677  692                  t2 = t4 + u * (two + z2 * t5);
 678  693          } else {
 679  694                  t5 = t4 + (u * z2) * t5;
 680  695                  u2 = u + u;
 681      -                v = (long double) ((int) (u2 * t24)) * p24;
      696 +                v = (long double)((int)(u2 * t24)) * p24;
 682  697                  t2 = t5 + r * ((two * t2 - v * t1) - v * (y - (t1 - z)));
 683  698                  t3 += v;
 684  699          }
      700 +
 685  701          ss_h = CHOPPED((t2 + t3));
 686  702          ss_l = t2 - (ss_h - t3);
 687      -/* INDENT OFF */
      703 +
 688  704  /*
 689  705   * compute ww = (x-.5)*(log(x)-1) + .5*(log(2pi)-1) + 1/x*(P(1/x^2)))
 690  706   * where ss = log(x) - 1 in already in extra precision
 691  707   */
 692      -        /* INDENT ON */
 693  708          z = one / x;
 694  709          r = x - half;
 695  710          r_h = CHOPPED((r));
 696  711          w_h = r_h * ss_h + hln2pim1_h;
 697  712          z2 = z * z;
 698  713          w = (r - r_h) * ss_h + r * ss_l;
 699  714          t1 = GP[19];
      715 +
 700  716          for (i = 18; i > 0; i--)
 701  717                  t1 = z2 * t1 + GP[i];
      718 +
 702  719          w += hln2pim1_l;
 703  720          w_l = z * (GP[0] + z2 * t1) + w;
 704      -        k = (int) ((w_h + w_l) * invln2_32 + half);
      721 +        k = (int)((w_h + w_l) * invln2_32 + half);
 705  722  
 706  723          /* compute the exponential of w_h+w_l */
 707  724  
 708  725          j = k & 0x1f;
 709  726          *m = k >> 5;
 710      -        t3 = (long double) k;
      727 +        t3 = (long double)k;
 711  728  
 712  729          /* perform w - k*ln2_32 (represent as w_h - w_l) */
 713  730          t1 = w_h - t3 * ln2_32hi;
 714  731          t2 = t3 * ln2_32lo;
 715  732          w = t2 - w_l;
 716  733          w_h = t1 - w;
 717  734          w_l = w - (t1 - w_h);
 718  735  
 719  736          /* compute exp(w_h-w_l) */
 720  737          z = w_h - w_l;
      738 +
 721  739          for (t1 = Et[10], i = 9; i >= 0; i--)
 722  740                  t1 = z * t1 + Et[i];
      741 +
 723  742          t3 = w_h - (w_l - (z * z) * t1);        /* t3 = expm1(z) */
 724  743          zz.l = S_trail[j] * (one + t3) + S[j] * t3;
 725  744          zz.h = S[j];
 726  745          return (zz);
 727  746  }
 728  747  
 729      -/* INDENT OFF */
      748 +
 730  749  /*
 731  750   * kpsin(x)= sin(pi*x)/pi
 732  751   *                 3        5        7        9        11                27
 733  752   *      = x+ks[0]*x +ks[1]*x +ks[2]*x +ks[3]*x +ks[4]*x  + ... + ks[12]*x
 734  753   */
 735  754  static const long double ks[] = {
 736  755          -1.64493406684822643647241516664602518705158902870e+0000L,
 737  756          +8.11742425283353643637002772405874238094995726160e-0001L,
 738  757          -1.90751824122084213696472111835337366232282723933e-0001L,
 739  758          +2.61478478176548005046532613563241288115395517084e-0002L,
 740  759          -2.34608103545582363750893072647117829448016479971e-0003L,
 741  760          +1.48428793031071003684606647212534027556262040158e-0004L,
 742  761          -6.97587366165638046518462722252768122615952898698e-0006L,
 743  762          +2.53121740413702536928659271747187500934840057929e-0007L,
 744  763          -7.30471182221385990397683641695766121301933621956e-0009L,
 745  764          +1.71653847451163495739958249695549313987973589884e-0010L,
 746  765          -3.34813314714560776122245796929054813458341420565e-0012L,
 747  766          +5.50724992262622033449487808306969135431411753047e-0014L,
 748  767          -7.67678132753577998601234393215802221104236979928e-0016L,
 749  768  };
 750      -/* INDENT ON */
 751  769  
 752  770  /*
 753  771   * assume x is not tiny and positive
 754  772   */
 755  773  static struct LDouble
 756      -kpsin(long double x) {
      774 +kpsin(long double x)
      775 +{
 757  776          long double z, t1, t2;
 758  777          struct LDouble xx;
 759  778          int i;
 760  779  
 761  780          z = x * x;
 762  781          xx.h = x;
      782 +
 763  783          for (t2 = ks[12], i = 11; i > 0; i--)
 764  784                  t2 = z * t2 + ks[i];
      785 +
 765  786          t1 = z * x;
 766  787          t2 *= z * t1;
 767  788          xx.l = t1 * ks[0] + t2;
 768  789          return (xx);
 769  790  }
 770  791  
 771      -/* INDENT OFF */
      792 +
 772  793  /*
 773  794   * kpcos(x)= cos(pi*x)/pi
 774  795   *                     2        4        6        8        10        12
 775  796   *      = 1/pi +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x +kc[4]*x  +kc[5]*x
 776  797   *
 777  798   *                     2        4        6        8        10            22
 778  799   *      = 1/pi - pi/2*x +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x  +...+kc[9]*x
 779  800   *
 780  801   * -pi/2*x*x = (npi_2_h + npi_2_l) * (x_f+x_l)*(x_f+x_l)
 781  802   *         =  npi_2_h*(x_f+x_l)*(x_f+x_l) + npi_2_l*x*x
↓ open down ↓ 10 lines elided ↑ open up ↑
 792  813   *  -1.5893254773528196691639751442098584699687552910487472296153e-8
 793  814   * 1/pi(in hex) =
 794  815   *  .517CC1B727220A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D5EF5DE2B
 795  816   * will be splitted into:
 796  817   *  one_pi_h = 1/pi chopped to 48 bits = .517CC1B727220000000000...  and
 797  818   *  one_pi_l = .0000000000000A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D5EF5DE2B
 798  819   */
 799  820  
 800  821  static const long double
 801  822  #if defined(__x86)
 802      -one_pi_h = 0.3183098861481994390487670898437500L,       /* 31 bits */
 803      -one_pi_l = 3.559123248900043690127872406891929148e-11L,
      823 +        one_pi_h = 0.3183098861481994390487670898437500L,       /* 31 bits */
      824 +        one_pi_l = 3.559123248900043690127872406891929148e-11L,
 804  825  #else
 805      -one_pi_h = 0.31830988618379052468299050815403461456298828125L,
 806      -one_pi_l = 1.46854777018590994109505931010230912897495334688117e-16L,
      826 +        one_pi_h = 0.31830988618379052468299050815403461456298828125L,
      827 +        one_pi_l = 1.46854777018590994109505931010230912897495334688117e-16L,
 807  828  #endif
 808      -npi_2_h = -1.570796310901641845703125000000000L,
 809      -npi_2_l = -1.5893254773528196691639751442098584699687552910e-8L;
 810      -
      829 +        npi_2_h = -1.570796310901641845703125000000000L,
      830 +        npi_2_l = -1.5893254773528196691639751442098584699687552910e-8L;
 811  831  static const long double kc[] = {
 812  832          +1.29192819501249250731151312779548918765320728489e+0000L,
 813  833          -4.25027339979557573976029596929319207009444090366e-0001L,
 814  834          +7.49080661650990096109672954618317623888421628613e-0002L,
 815  835          -8.21458866111282287985539464173976555436050215120e-0003L,
 816  836          +6.14202578809529228503205255165761204750211603402e-0004L,
 817  837          -3.33073432691149607007217330302595267179545908740e-0005L,
 818  838          +1.36970959047832085796809745461530865597993680204e-0006L,
 819  839          -4.41780774262583514450246512727201806217271097336e-0008L,
 820  840          +1.14741409212381858820016567664488123478660705759e-0009L,
 821  841          -2.44261236114707374558437500654381006300502749632e-0011L,
 822  842  };
 823      -/* INDENT ON */
 824  843  
 825  844  /*
 826  845   * assume x is not tiny and positive
 827  846   */
 828  847  static struct LDouble
 829      -kpcos(long double x) {
      848 +kpcos(long double x)
      849 +{
 830  850          long double z, t1, t2, t3, t4, x4, x8;
 831  851          int i;
 832  852          struct LDouble xx;
 833  853  
 834  854          z = x * x;
 835  855          xx.h = one_pi_h;
 836      -        t1 = (long double) ((float) x);
      856 +        t1 = (long double)((float)x);
 837  857          x4 = z * z;
 838  858          t2 = npi_2_l * z + npi_2_h * (x + t1) * (x - t1);
      859 +
 839  860          for (i = 8, t3 = kc[9]; i >= 0; i--)
 840  861                  t3 = z * t3 + kc[i];
      862 +
 841  863          t3 = one_pi_l + x4 * t3;
 842  864          t4 = t1 * t1 * npi_2_h;
 843  865          x8 = t2 + t3;
 844  866          xx.l = x8 + t4;
 845  867          return (xx);
 846  868  }
 847  869  
 848      -/* INDENT OFF */
 849  870  static const long double
 850      -        /* 0.13486180573279076968979393577465291700642511139552429398233 */
      871 +/* 0.13486180573279076968979393577465291700642511139552429398233 */
 851  872  #if defined(__x86)
 852      -t0z1   =  0.1348618057327907696779385054997035808810L,
 853      -t0z1_l =  1.1855430274949336125392717150257379614654e-20L,
      873 +        t0z1 = 0.1348618057327907696779385054997035808810L,
      874 +        t0z1_l = 1.1855430274949336125392717150257379614654e-20L,
 854  875  #else
 855      -t0z1   =  0.1348618057327907696897939357746529168654L,
 856      -t0z1_l =  1.4102088588676879418739164486159514674310e-37L,
      876 +        t0z1 = 0.1348618057327907696897939357746529168654L,
      877 +        t0z1_l = 1.4102088588676879418739164486159514674310e-37L,
 857  878  #endif
 858      -        /* 0.46163214496836234126265954232572132846819620400644635129599 */
      879 +/* 0.46163214496836234126265954232572132846819620400644635129599 */
 859  880  #if defined(__x86)
 860      -t0z2   =  0.4616321449683623412538115843295472018326L,
 861      -t0z2_l =  8.84795799617412663558532305039261747030640e-21L,
      881 +        t0z2 = 0.4616321449683623412538115843295472018326L,
      882 +        t0z2_l = 8.84795799617412663558532305039261747030640e-21L,
 862  883  #else
 863      -t0z2   =  0.46163214496836234126265954232572132343318L,
 864      -t0z2_l =  5.03501162329616380465302666480916271611101e-36L,
      884 +        t0z2 = 0.46163214496836234126265954232572132343318L,
      885 +        t0z2_l = 5.03501162329616380465302666480916271611101e-36L,
 865  886  #endif
 866      -        /* 0.81977310110050060178786870492160699631174407846245179119586 */
      887 +/* 0.81977310110050060178786870492160699631174407846245179119586 */
 867  888  #if defined(__x86)
 868      -t0z3   =  0.81977310110050060178773362329351925836817L,
 869      -t0z3_l =  1.350816280877379435658077052534574556256230e-22L
      889 +        t0z3 = 0.81977310110050060178773362329351925836817L,
      890 +        t0z3_l = 1.350816280877379435658077052534574556256230e-22L
 870  891  #else
 871      -t0z3   =  0.8197731011005006017878687049216069516957449L,
 872      -t0z3_l =  4.461599916947014419045492615933551648857380e-35L
      892 +        t0z3 = 0.8197731011005006017878687049216069516957449L,
      893 +        t0z3_l = 4.461599916947014419045492615933551648857380e-35L
 873  894  #endif
 874  895  ;
 875      -/* INDENT ON */
 876  896  
 877  897  /*
 878  898   * gamma(x+i) for 0 <= x < 1
 879  899   */
 880  900  static struct LDouble
 881      -gam_n(int i, long double x) {
 882      -        struct LDouble rr = {0.0L, 0.0L}, yy;
      901 +gam_n(int i, long double x)
      902 +{
      903 +        struct LDouble rr = { 0.0L, 0.0L }, yy;
 883  904          long double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl;
 884  905  
 885  906          /* compute yy = gamma(x+1) */
 886  907          if (x > 0.2845L) {
 887  908                  if (x > 0.6374L) {
 888  909                          r1 = x - t0z3;
 889  910                          r2 = CHOPPED((r1 - t0z3_l));
 890  911                          t2 = r1 - r2;
 891  912                          yy = GT3(r2, t2 - t0z3_l);
 892  913                  } else {
↓ open down ↓ 1 lines elided ↑ open up ↑
 894  915                          r2 = CHOPPED((r1 - t0z2_l));
 895  916                          t2 = r1 - r2;
 896  917                          yy = GT2(r2, t2 - t0z2_l);
 897  918                  }
 898  919          } else {
 899  920                  r1 = x - t0z1;
 900  921                  r2 = CHOPPED((r1 - t0z1_l));
 901  922                  t2 = r1 - r2;
 902  923                  yy = GT1(r2, t2 - t0z1_l);
 903  924          }
      925 +
 904  926          /* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */
 905  927          switch (i) {
 906      -        case 0:         /* yy/x */
      928 +        case 0:                         /* yy/x */
 907  929                  r1 = one / x;
 908  930                  xh = CHOPPED((x));      /* x is not tiny */
 909  931                  rr.h = CHOPPED(((yy.h + yy.l) * r1));
 910      -                rr.l = r1 * (yy.h - rr.h * xh) - ((r1 * rr.h) * (x - xh) -
 911      -                        r1 * yy.l);
      932 +                rr.l = r1 * (yy.h - rr.h * xh) - ((r1 * rr.h) * (x - xh) - r1 *
      933 +                    yy.l);
 912  934                  break;
 913      -        case 1:         /* yy */
      935 +        case 1:                         /* yy */
 914  936                  rr.h = yy.h;
 915  937                  rr.l = yy.l;
 916  938                  break;
 917      -        case 2:         /* (x+1)*yy */
 918      -                z = x + one;    /* may not be exact */
      939 +        case 2:                         /* (x+1)*yy */
      940 +                z = x + one;            /* may not be exact */
 919  941                  zh = CHOPPED((z));
 920  942                  rr.h = zh * yy.h;
 921  943                  rr.l = z * yy.l + (x - (zh - one)) * yy.h;
 922  944                  break;
 923      -        case 3:         /* (x+2)*(x+1)*yy */
      945 +        case 3:                         /* (x+2)*(x+1)*yy */
 924  946                  z1 = x + one;
 925  947                  z2 = x + 2.0L;
 926  948                  z = z1 * z2;
 927  949                  xh = CHOPPED((z));
 928  950                  zh = CHOPPED((z1));
 929  951                  xl = (x - (zh - one)) * (z2 + zh) - (xh - zh * (zh + one));
 930  952  
 931  953                  rr.h = xh * yy.h;
 932  954                  rr.l = z * yy.l + xl * yy.h;
 933  955                  break;
 934  956  
 935      -        case 4:         /* (x+1)*(x+3)*(x+2)*yy */
      957 +        case 4:                         /* (x+1)*(x+3)*(x+2)*yy */
 936  958                  z1 = x + 2.0L;
 937  959                  z2 = (x + one) * (x + 3.0L);
 938  960                  zh = CHOPPED(z1);
 939  961                  zl = x - (zh - 2.0L);
 940  962                  xh = CHOPPED(z2);
 941  963                  xl = zl * (zh + z1) - (xh - (zh * zh - one));
 942  964  
 943  965                  /* wh+wl=(x+2)*yy */
 944  966                  wh = CHOPPED((z1 * (yy.h + yy.l)));
 945  967                  wl = (zl * yy.h + z1 * yy.l) - (wh - zh * yy.h);
 946  968  
 947  969                  rr.h = xh * wh;
 948  970                  rr.l = z2 * wl + xl * wh;
 949  971  
 950  972                  break;
 951      -        case 5:         /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */
      973 +        case 5:                         /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */
 952  974                  z1 = x + 2.0L;
 953  975                  z2 = x + 3.0L;
 954  976                  z = z1 * z2;
 955  977                  zh = CHOPPED((z1));
 956  978                  yh = CHOPPED((z));
 957  979                  yl = (x - (zh - 2.0L)) * (z2 + zh) - (yh - zh * (zh + one));
 958  980                  z2 = z - 2.0L;
 959  981                  z *= z2;
 960  982                  xh = CHOPPED((z));
 961  983                  xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L));
 962  984                  rr.h = xh * yy.h;
 963  985                  rr.l = z * yy.l + xl * yy.h;
 964  986                  break;
 965      -        case 6:         /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
      987 +        case 6:                         /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
 966  988                  z1 = x + 2.0L;
 967  989                  z2 = x + 3.0L;
 968  990                  z = z1 * z2;
 969  991                  zh = CHOPPED((z1));
 970  992                  yh = CHOPPED((z));
 971  993                  z1 = x - (zh - 2.0L);
 972  994                  yl = z1 * (z2 + zh) - (yh - zh * (zh + one));
 973  995                  z2 = z - 2.0L;
 974  996                  x5 = x + 5.0L;
 975  997                  z *= z2;
 976  998                  xh = CHOPPED(z);
 977  999                  zh += 3.0;
 978 1000                  xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L));
 979      -                                                /* xh+xl=(x+1)*...*(x+4) */
 980      -                /* wh+wl=(x+5)*yy */
     1001 +
     1002 +                /*
     1003 +                 * xh+xl=(x+1)*...*(x+4)
     1004 +                 * wh+wl=(x+5)*yy
     1005 +                 */
 981 1006                  wh = CHOPPED((x5 * (yy.h + yy.l)));
 982 1007                  wl = (z1 * yy.h + x5 * yy.l) - (wh - zh * yy.h);
 983 1008                  rr.h = wh * xh;
 984 1009                  rr.l = z * wl + xl * wh;
 985 1010                  break;
 986      -        case 7:         /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
     1011 +        case 7:                 /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
 987 1012                  z1 = x + 3.0L;
 988 1013                  z2 = x + 4.0L;
 989 1014                  z = z2 * z1;
 990 1015                  zh = CHOPPED((z1));
 991 1016                  yh = CHOPPED((z));      /* yh+yl = (x+3)(x+4) */
 992 1017                  yl = (x - (zh - 3.0L)) * (z2 + zh) - (yh - (zh * (zh + one)));
 993 1018                  z1 = x + 6.0L;
 994      -                z2 = z - 2.0L;  /* z2 = (x+2)*(x+5) */
     1019 +                z2 = z - 2.0L;          /* z2 = (x+2)*(x+5) */
 995 1020                  z *= z2;
 996 1021                  xh = CHOPPED((z));
 997 1022                  xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L));
 998      -                                                /* xh+xl=(x+2)*...*(x+5) */
 999      -                /* wh+wl=(x+1)(x+6)*yy */
1000      -                z2 -= 4.0L;     /* z2 = (x+1)(x+6) */
     1023 +
     1024 +                /*
     1025 +                 * xh+xl=(x+2)*...*(x+5)
     1026 +                 * wh+wl=(x+1)(x+6)*yy
     1027 +                 */
     1028 +                z2 -= 4.0L;             /* z2 = (x+1)(x+6) */
1001 1029                  wh = CHOPPED((z2 * (yy.h + yy.l)));
1002 1030                  wl = (z2 * yy.l + yl * yy.h) - (wh - (yh - 6.0L) * yy.h);
1003 1031                  rr.h = wh * xh;
1004 1032                  rr.l = z * wl + xl * wh;
1005 1033          }
     1034 +
1006 1035          return (rr);
1007 1036  }
1008 1037  
1009 1038  long double
1010      -tgammal(long double x) {
     1039 +tgammal(long double x)
     1040 +{
1011 1041          struct LDouble ss, ww;
1012 1042          long double t, t1, t2, t3, t4, t5, w, y, z, z1, z2, z3, z5;
1013 1043          int i, j, m, ix, hx, xk;
1014 1044          unsigned lx;
1015 1045  
1016 1046          hx = H0_WORD(x);
1017 1047          lx = H3_WORD(x);
1018 1048          ix = hx & 0x7fffffff;
1019 1049          y = x;
1020      -        if (ix < 0x3f8e0000) {  /* x < 2**-113 */
     1050 +
     1051 +        if (ix < 0x3f8e0000)            /* x < 2**-113 */
1021 1052                  return (one / x);
1022      -        }
     1053 +
1023 1054          if (ix >= 0x7fff0000)
1024      -                return (x * ((hx < 0)? zero : x));      /* Inf or NaN */
1025      -        if (x > overflow)       /* overflow threshold */
     1055 +                return (x * ((hx < 0) ? zero : x));     /* Inf or NaN */
     1056 +
     1057 +        if (x > overflow)                               /* overflow threshold */
1026 1058                  return (x * 1.0e4932L);
1027      -        if (hx >= 0x40020000) { /* x >= 8 */
     1059 +
     1060 +        if (hx >= 0x40020000) {                         /* x >= 8 */
1028 1061                  ww = large_gam(x, &m);
1029 1062                  w = ww.h + ww.l;
1030 1063                  return (scalbnl(w, m));
1031 1064          }
1032 1065  
1033      -        if (hx > 0) {           /* 0 < x < 8 */
1034      -                i = (int) x;
1035      -                ww = gam_n(i, x - (long double) i);
     1066 +        if (hx > 0) {                   /* 0 < x < 8 */
     1067 +                i = (int)x;
     1068 +                ww = gam_n(i, x - (long double)i);
1036 1069                  return (ww.h + ww.l);
1037 1070          }
1038      -        /* INDENT OFF */
1039      -        /* negative x */
     1071 +
     1072 +        /*
     1073 +         * negative x
     1074 +         */
     1075 +
1040 1076          /*
1041 1077           * compute xk =
1042 1078           *      -2 ... x is an even int (-inf is considered an even #)
1043 1079           *      -1 ... x is an odd int
1044 1080           *      +0 ... x is not an int but chopped to an even int
1045 1081           *      +1 ... x is not an int but chopped to an odd int
1046 1082           */
1047      -        /* INDENT ON */
1048 1083          xk = 0;
1049 1084  #if defined(__x86)
1050      -        if (ix >= 0x403e0000) { /* x >= 2**63 } */
     1085 +        if (ix >= 0x403e0000) {         /* x >= 2**63 } */
1051 1086                  if (ix >= 0x403f0000)
1052 1087                          xk = -2;
1053 1088                  else
1054 1089                          xk = -2 + (lx & 1);
1055 1090  #else
1056      -        if (ix >= 0x406f0000) { /* x >= 2**112 */
     1091 +        if (ix >= 0x406f0000) {         /* x >= 2**112 */
1057 1092                  if (ix >= 0x40700000)
1058 1093                          xk = -2;
1059 1094                  else
1060 1095                          xk = -2 + (lx & 1);
1061 1096  #endif
1062 1097          } else if (ix >= 0x3fff0000) {
1063 1098                  w = -x;
1064 1099                  t1 = floorl(w);
1065 1100                  t2 = t1 * half;
1066 1101                  t3 = floorl(t2);
     1102 +
1067 1103                  if (t1 == w) {
1068 1104                          if (t2 == t3)
1069 1105                                  xk = -2;
1070 1106                          else
1071 1107                                  xk = -1;
1072 1108                  } else {
1073 1109                          if (t2 == t3)
1074 1110                                  xk = 0;
1075 1111                          else
1076 1112                                  xk = 1;
1077 1113                  }
1078 1114          }
1079 1115  
1080 1116          if (xk < 0) {
1081 1117                  /* return NaN. Ideally gamma(-n)= (-1)**(n+1) * inf */
1082      -                return (x - x) / (x - x);
     1118 +                return ((x - x) / (x - x));
1083 1119          }
1084 1120  
1085 1121          /*
1086 1122           * negative underflow thresold -(1774+9ulp)
1087 1123           */
1088 1124          if (x < -1774.0000000000000000000000000000017749370L) {
1089 1125                  z = tiny / x;
     1126 +
1090 1127                  if (xk == 1)
1091 1128                          z = -z;
     1129 +
1092 1130                  return (z * tiny);
1093 1131          }
1094 1132  
1095      -        /* INDENT OFF */
     1133 +
1096 1134          /*
1097 1135           * now compute gamma(x) by  -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x
1098 1136           */
     1137 +
1099 1138          /*
1100 1139           * First compute ss = -sin(pi*y)/pi so that
1101 1140           * gamma(x) = 1/(ss*gamma(1+y))
1102 1141           */
1103      -        /* INDENT ON */
1104 1142          y = -x;
1105      -        j = (int) y;
1106      -        z = y - (long double) j;
1107      -        if (z > 0.3183098861837906715377675L)
     1143 +        j = (int)y;
     1144 +        z = y - (long double)j;
     1145 +
     1146 +        if (z > 0.3183098861837906715377675L) {
1108 1147                  if (z > 0.6816901138162093284622325L)
1109 1148                          ss = kpsin(one - z);
1110 1149                  else
1111 1150                          ss = kpcos(0.5L - z);
1112      -        else
     1151 +        } else {
1113 1152                  ss = kpsin(z);
     1153 +        }
     1154 +
1114 1155          if (xk == 0) {
1115 1156                  ss.h = -ss.h;
1116 1157                  ss.l = -ss.l;
1117 1158          }
1118 1159  
1119 1160          /* Then compute ww = gamma(1+y), note that result scale to 2**m */
1120 1161          m = 0;
     1162 +
1121 1163          if (j < 7) {
1122 1164                  ww = gam_n(j + 1, z);
1123 1165          } else {
1124 1166                  w = y + one;
     1167 +
1125 1168                  if ((lx & 1) == 0) {    /* y+1 exact (note that y<184) */
1126 1169                          ww = large_gam(w, &m);
1127 1170                  } else {
1128 1171                          t = w - one;
     1172 +
1129 1173                          if (t == y) {   /* y+one exact */
1130 1174                                  ww = large_gam(w, &m);
1131 1175                          } else {        /* use y*gamma(y) */
1132 1176                                  if (j == 7)
1133 1177                                          ww = gam_n(j, z);
1134 1178                                  else
1135 1179                                          ww = large_gam(y, &m);
     1180 +
1136 1181                                  t4 = ww.h + ww.l;
1137 1182                                  t1 = CHOPPED((y));
1138 1183                                  t2 = CHOPPED((t4));
1139      -                                                /* t4 will not be too large */
     1184 +                                /* t4 will not be too large */
1140 1185                                  ww.l = y * (ww.l - (t2 - ww.h)) + (y - t1) * t2;
1141 1186                                  ww.h = t1 * t2;
1142 1187                          }
1143 1188                  }
1144 1189          }
1145 1190  
1146 1191          /* compute 1/(ss*ww) */
1147 1192          t3 = ss.h + ss.l;
1148 1193          t4 = ww.h + ww.l;
1149 1194          t1 = CHOPPED((t3));
↓ open down ↓ 15 lines elided ↑ open up ↑
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX