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

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/m9x/tgammaf.c
          +++ new/usr/src/lib/libm/common/m9x/tgammaf.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 __tgammaf = tgammaf
  31   32  
  32   33  /*
  33   34   * True gamma function
  34   35   *
↓ open down ↓ 5 lines elided ↑ open up ↑
  40   41   */
  41   42  
  42   43  #include "libm.h"
  43   44  #include <math.h>
  44   45  #if defined(__SUNPRO_C)
  45   46  #include <sunmath.h>
  46   47  #endif
  47   48  #include <sys/isa_defs.h>
  48   49  
  49   50  #if defined(_BIG_ENDIAN)
  50      -#define HIWORD  0
  51      -#define LOWORD  1
       51 +#define HIWORD          0
       52 +#define LOWORD          1
  52   53  #else
  53      -#define HIWORD  1
  54      -#define LOWORD  0
       54 +#define HIWORD          1
       55 +#define LOWORD          0
  55   56  #endif
  56      -#define __HI(x) ((int *) &x)[HIWORD]
  57      -#define __LO(x) ((unsigned *) &x)[LOWORD]
       57 +#define __HI(x)         ((int *)&x)[HIWORD]
       58 +#define __LO(x)         ((unsigned *)&x)[LOWORD]
  58   59  
  59   60  /* Coefficients for primary intervals GTi() */
  60   61  static const double cr[] = {
  61   62          /* p1 */
  62   63          +7.09087253435088360271451613398019280077561279443e-0001,
  63   64          -5.17229560788652108545141978238701790105241761089e-0001,
  64   65          +5.23403394528150789405825222323770647162337764327e-0001,
  65   66          -4.54586308717075010784041566069480411732634814899e-0001,
  66   67          +4.20596490915239085459964590559256913498190955233e-0001,
  67   68          -3.57307589712377520978332185838241458642142185789e-0001,
  68      -
  69   69          /* p2 */
  70   70          +4.28486983980295198166056119223984284434264344578e-0001,
  71   71          -1.30704539487709138528680121627899735386650103914e-0001,
  72   72          +1.60856285038051955072861219352655851542955430871e-0001,
  73   73          -9.22285161346010583774458802067371182158937943507e-0002,
  74   74          +7.19240511767225260740890292605070595560626179357e-0002,
  75   75          -4.88158265593355093703112238534484636193260459574e-0002,
  76      -
  77   76          /* p3 */
  78   77          +3.82409531118807759081121479786092134814808872880e-0001,
  79   78          +2.65309888180188647956400403013495759365167853426e-0002,
  80   79          +8.06815109775079171923561169415370309376296739835e-0002,
  81   80          -1.54821591666137613928840890835174351674007764799e-0002,
  82   81          +1.76308239242717268530498313416899188157165183405e-0002,
  83      -
  84   82          /* GZi and TZi */
  85   83          +0.9382046279096824494097535615803269576988,    /* GZ1 */
  86   84          +0.8856031944108887002788159005825887332080,    /* GZ2 */
  87   85          +0.9367814114636523216188468970808378497426,    /* GZ3 */
  88      -        -0.3517214357852935791015625,   /* TZ1 */
  89      -        +0.280530631542205810546875,    /* TZ3 */
       86 +        -0.3517214357852935791015625,                   /* TZ1 */
       87 +        +0.280530631542205810546875,                    /* TZ3 */
  90   88  };
  91   89  
  92      -#define P10     cr[0]
  93      -#define P11     cr[1]
  94      -#define P12     cr[2]
  95      -#define P13     cr[3]
  96      -#define P14     cr[4]
  97      -#define P15     cr[5]
  98      -#define P20     cr[6]
  99      -#define P21     cr[7]
 100      -#define P22     cr[8]
 101      -#define P23     cr[9]
 102      -#define P24     cr[10]
 103      -#define P25     cr[11]
 104      -#define P30     cr[12]
 105      -#define P31     cr[13]
 106      -#define P32     cr[14]
 107      -#define P33     cr[15]
 108      -#define P34     cr[16]
 109      -#define GZ1     cr[17]
 110      -#define GZ2     cr[18]
 111      -#define GZ3     cr[19]
 112      -#define TZ1     cr[20]
 113      -#define TZ3     cr[21]
       90 +#define P10             cr[0]
       91 +#define P11             cr[1]
       92 +#define P12             cr[2]
       93 +#define P13             cr[3]
       94 +#define P14             cr[4]
       95 +#define P15             cr[5]
       96 +#define P20             cr[6]
       97 +#define P21             cr[7]
       98 +#define P22             cr[8]
       99 +#define P23             cr[9]
      100 +#define P24             cr[10]
      101 +#define P25             cr[11]
      102 +#define P30             cr[12]
      103 +#define P31             cr[13]
      104 +#define P32             cr[14]
      105 +#define P33             cr[15]
      106 +#define P34             cr[16]
      107 +#define GZ1             cr[17]
      108 +#define GZ2             cr[18]
      109 +#define GZ3             cr[19]
      110 +#define TZ1             cr[20]
      111 +#define TZ3             cr[21]
 114  112  
 115  113  /* compute gamma(y) for y in GT1 = [1.0000, 1.2845] */
 116  114  static double
 117      -GT1(double y) {
      115 +GT1(double y)
      116 +{
 118  117          double z, r;
 119  118  
 120  119          z = y * y;
 121  120          r = TZ1 * y + z * ((P10 + y * P11 + z * P12) + (z * y) * (P13 + y *
 122      -                P14 + z * P15));
      121 +            P14 + z * P15));
 123  122          return (GZ1 + r);
 124  123  }
 125  124  
 126  125  /* compute gamma(y) for y in GT2 = [1.2844, 1.6374] */
 127  126  static double
 128      -GT2(double y) {
      127 +GT2(double y)
      128 +{
 129  129          double z;
 130  130  
 131  131          z = y * y;
 132  132          return (GZ2 + z * ((P20 + y * P21 + z * P22) + (z * y) * (P23 + y *
 133      -                P24 + z * P25)));
      133 +            P24 + z * P25)));
 134  134  }
 135  135  
 136  136  /* compute gamma(y) for y in GT3 = [1.6373, 2.0000] */
 137  137  static double
 138      -GT3(double y) {
 139      -double z, r;
      138 +GT3(double y)
      139 +{
      140 +        double z, r;
 140  141  
 141  142          z = y * y;
 142  143          r = TZ3 * y + z * ((P30 + y * P31 + z * P32) + (z * y) * (P33 + y *
 143      -                P34));
      144 +            P34));
 144  145          return (GZ3 + r);
 145  146  }
 146  147  
 147      -/* INDENT OFF */
 148  148  static const double c[] = {
 149      -+1.0,
 150      -+2.0,
 151      -+0.5,
 152      -+1.0e-300,
 153      -+6.666717231848518054693623697539230e-0001,                     /* A1=T3[0] */
 154      -+8.33333330959694065245736888749042811909994573178e-0002,       /* GP[0] */
 155      --2.77765545601667179767706600890361535225507762168e-0003,       /* GP[1] */
 156      -+7.77830853479775281781085278324621033523037489883e-0004,       /* GP[2] */
 157      -+4.18938533204672741744150788368695779923320328369e-0001,       /* hln2pi   */
 158      -+2.16608493924982901946e-02,                                    /* ln2_32 */
 159      -+4.61662413084468283841e+01,                                    /* invln2_32 */
 160      -+5.00004103388988968841156421415669985414073453720e-0001,       /* Et1 */
 161      -+1.66667656752800761782778277828110208108687545908e-0001,       /* Et2 */
      149 +        +1.0,
      150 +        +2.0,
      151 +        +0.5,
      152 +        +1.0e-300,
      153 +        +6.666717231848518054693623697539230e-0001,     /* A1=T3[0] */
      154 +        +8.33333330959694065245736888749042811909994573178e-0002, /* GP[0] */
      155 +        -2.77765545601667179767706600890361535225507762168e-0003, /* GP[1] */
      156 +        +7.77830853479775281781085278324621033523037489883e-0004, /* GP[2] */
      157 +        +4.18938533204672741744150788368695779923320328369e-0001, /* hln2pi   */
      158 +        +2.16608493924982901946e-02, /* ln2_32 */
      159 +        +4.61662413084468283841e+01, /* invln2_32 */
      160 +        +5.00004103388988968841156421415669985414073453720e-0001, /* Et1 */
      161 +        +1.66667656752800761782778277828110208108687545908e-0001, /* Et2 */
 162  162  };
 163  163  
 164  164  #define one             c[0]
 165  165  #define two             c[1]
 166  166  #define half            c[2]
 167  167  #define tiny            c[3]
 168  168  #define A1              c[4]
 169  169  #define GP0             c[5]
 170  170  #define GP1             c[6]
 171  171  #define GP2             c[7]
 172  172  #define hln2pi          c[8]
 173  173  #define ln2_32          c[9]
 174  174  #define invln2_32       c[10]
 175  175  #define Et1             c[11]
 176  176  #define Et2             c[12]
 177  177  
 178  178  /* S[j] = 2**(j/32.) for the final computation of exp(w) */
 179  179  static const double S[] = {
 180      -+1.00000000000000000000e+00,    /* 3FF0000000000000 */
 181      -+1.02189714865411662714e+00,    /* 3FF059B0D3158574 */
 182      -+1.04427378242741375480e+00,    /* 3FF0B5586CF9890F */
 183      -+1.06714040067682369717e+00,    /* 3FF11301D0125B51 */
 184      -+1.09050773266525768967e+00,    /* 3FF172B83C7D517B */
 185      -+1.11438674259589243221e+00,    /* 3FF1D4873168B9AA */
 186      -+1.13878863475669156458e+00,    /* 3FF2387A6E756238 */
 187      -+1.16372485877757747552e+00,    /* 3FF29E9DF51FDEE1 */
 188      -+1.18920711500272102690e+00,    /* 3FF306FE0A31B715 */
 189      -+1.21524735998046895524e+00,    /* 3FF371A7373AA9CB */
 190      -+1.24185781207348400201e+00,    /* 3FF3DEA64C123422 */
 191      -+1.26905095719173321989e+00,    /* 3FF44E086061892D */
 192      -+1.29683955465100964055e+00,    /* 3FF4BFDAD5362A27 */
 193      -+1.32523664315974132322e+00,    /* 3FF5342B569D4F82 */
 194      -+1.35425554693689265129e+00,    /* 3FF5AB07DD485429 */
 195      -+1.38390988196383202258e+00,    /* 3FF6247EB03A5585 */
 196      -+1.41421356237309514547e+00,    /* 3FF6A09E667F3BCD */
 197      -+1.44518080697704665027e+00,    /* 3FF71F75E8EC5F74 */
 198      -+1.47682614593949934623e+00,    /* 3FF7A11473EB0187 */
 199      -+1.50916442759342284141e+00,    /* 3FF82589994CCE13 */
 200      -+1.54221082540794074411e+00,    /* 3FF8ACE5422AA0DB */
 201      -+1.57598084510788649659e+00,    /* 3FF93737B0CDC5E5 */
 202      -+1.61049033194925428347e+00,    /* 3FF9C49182A3F090 */
 203      -+1.64575547815396494578e+00,    /* 3FFA5503B23E255D */
 204      -+1.68179283050742900407e+00,    /* 3FFAE89F995AD3AD */
 205      -+1.71861929812247793414e+00,    /* 3FFB7F76F2FB5E47 */
 206      -+1.75625216037329945351e+00,    /* 3FFC199BDD85529C */
 207      -+1.79470907500310716820e+00,    /* 3FFCB720DCEF9069 */
 208      -+1.83400808640934243066e+00,    /* 3FFD5818DCFBA487 */
 209      -+1.87416763411029996256e+00,    /* 3FFDFC97337B9B5F */
 210      -+1.91520656139714740007e+00,    /* 3FFEA4AFA2A490DA */
 211      -+1.95714412417540017941e+00,    /* 3FFF50765B6E4540 */
      180 +        +1.00000000000000000000e+00,    /* 3FF0000000000000 */
      181 +        +1.02189714865411662714e+00,    /* 3FF059B0D3158574 */
      182 +        +1.04427378242741375480e+00,    /* 3FF0B5586CF9890F */
      183 +        +1.06714040067682369717e+00,    /* 3FF11301D0125B51 */
      184 +        +1.09050773266525768967e+00,    /* 3FF172B83C7D517B */
      185 +        +1.11438674259589243221e+00,    /* 3FF1D4873168B9AA */
      186 +        +1.13878863475669156458e+00,    /* 3FF2387A6E756238 */
      187 +        +1.16372485877757747552e+00,    /* 3FF29E9DF51FDEE1 */
      188 +        +1.18920711500272102690e+00,    /* 3FF306FE0A31B715 */
      189 +        +1.21524735998046895524e+00,    /* 3FF371A7373AA9CB */
      190 +        +1.24185781207348400201e+00,    /* 3FF3DEA64C123422 */
      191 +        +1.26905095719173321989e+00,    /* 3FF44E086061892D */
      192 +        +1.29683955465100964055e+00,    /* 3FF4BFDAD5362A27 */
      193 +        +1.32523664315974132322e+00,    /* 3FF5342B569D4F82 */
      194 +        +1.35425554693689265129e+00,    /* 3FF5AB07DD485429 */
      195 +        +1.38390988196383202258e+00,    /* 3FF6247EB03A5585 */
      196 +        +1.41421356237309514547e+00,    /* 3FF6A09E667F3BCD */
      197 +        +1.44518080697704665027e+00,    /* 3FF71F75E8EC5F74 */
      198 +        +1.47682614593949934623e+00,    /* 3FF7A11473EB0187 */
      199 +        +1.50916442759342284141e+00,    /* 3FF82589994CCE13 */
      200 +        +1.54221082540794074411e+00,    /* 3FF8ACE5422AA0DB */
      201 +        +1.57598084510788649659e+00,    /* 3FF93737B0CDC5E5 */
      202 +        +1.61049033194925428347e+00,    /* 3FF9C49182A3F090 */
      203 +        +1.64575547815396494578e+00,    /* 3FFA5503B23E255D */
      204 +        +1.68179283050742900407e+00,    /* 3FFAE89F995AD3AD */
      205 +        +1.71861929812247793414e+00,    /* 3FFB7F76F2FB5E47 */
      206 +        +1.75625216037329945351e+00,    /* 3FFC199BDD85529C */
      207 +        +1.79470907500310716820e+00,    /* 3FFCB720DCEF9069 */
      208 +        +1.83400808640934243066e+00,    /* 3FFD5818DCFBA487 */
      209 +        +1.87416763411029996256e+00,    /* 3FFDFC97337B9B5F */
      210 +        +1.91520656139714740007e+00,    /* 3FFEA4AFA2A490DA */
      211 +        +1.95714412417540017941e+00,    /* 3FFF50765B6E4540 */
 212  212  };
 213      -/* INDENT ON */
 214  213  
 215      -/* INDENT OFF */
      214 +
 216  215  /*
 217  216   * return tgammaf(x) in double for 8<x<=35.040096283... using Stirling's formula
 218  217   *     log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x))
 219  218   */
      219 +
 220  220  /*
 221  221   * compute ss = log(x)-1
 222  222   *
 223  223   *  log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y,  1<=y<2,
 224  224   *  j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
 225  225   *       T1(n-3) = n*log(2)-1,  n=3,4,5
 226  226   *       T2(j) = log(z[j]),
 227  227   *       T3(s) = 2s + A1*s^3
 228  228   *  Note
 229  229   *  (1) Remez error for T3(s) is bounded by 2**(-35.8)
 230  230   *      (see mpremez/work/Log/tgamma_log_2_outr1)
 231  231   */
 232      -
 233      -static const double T1[] = { /* T1[j]=(j+3)*log(2)-1 */
 234      -+1.079441541679835928251696364375e+00,
 235      -+1.772588722239781237668928485833e+00,
 236      -+2.465735902799726547086160607291e+00,
      232 +static const double T1[] = {            /* T1[j]=(j+3)*log(2)-1 */
      233 +        +1.079441541679835928251696364375e+00,
      234 +        +1.772588722239781237668928485833e+00,
      235 +        +2.465735902799726547086160607291e+00,
 237  236  };
 238  237  
 239      -static const double T2[] = {   /* T2[j]=log(1+j/64+1/128) */
 240      -+7.782140442054948947462900061137e-03,
 241      -+2.316705928153437822879916096229e-02,
 242      -+3.831886430213659919375532512380e-02,
 243      -+5.324451451881228286587019378653e-02,
 244      -+6.795066190850774939456527777263e-02,
 245      -+8.244366921107459126816006866831e-02,
 246      -+9.672962645855111229557105648746e-02,
 247      -+1.108143663402901141948061693232e-01,
 248      -+1.247034785009572358634065153809e-01,
 249      -+1.384023228591191356853258736016e-01,
 250      -+1.519160420258419750718034248969e-01,
 251      -+1.652495728953071628756114492772e-01,
 252      -+1.784076574728182971194002415109e-01,
 253      -+1.913948529996294546092988075613e-01,
 254      -+2.042155414286908915038203861962e-01,
 255      -+2.168739383006143596190895257443e-01,
 256      -+2.293741010648458299914807250461e-01,
 257      -+2.417199368871451681443075159135e-01,
 258      -+2.539152099809634441373232979066e-01,
 259      -+2.659635484971379413391259265375e-01,
 260      -+2.778684510034563061863500329234e-01,
 261      -+2.896332925830426768788930555257e-01,
 262      -+3.012613305781617810128755382338e-01,
 263      -+3.127557100038968883862465596883e-01,
 264      -+3.241194686542119760906707604350e-01,
 265      -+3.353555419211378302571795798142e-01,
 266      -+3.464667673462085809184621884258e-01,
 267      -+3.574558889218037742260094901409e-01,
 268      -+3.683255611587076530482301540504e-01,
 269      -+3.790783529349694583908533456310e-01,
 270      -+3.897167511400252133704636040035e-01,
 271      -+4.002431641270127069293251019951e-01,
 272      -+4.106599249852683859343062031758e-01,
 273      -+4.209692946441296361288671615068e-01,
 274      -+4.311734648183713408591724789556e-01,
 275      -+4.412745608048752294894964416613e-01,
 276      -+4.512746441394585851446923830790e-01,
 277      -+4.611757151221701663679999255979e-01,
 278      -+4.709797152187910125468978560564e-01,
 279      -+4.806885293457519076766184554480e-01,
 280      -+4.903039880451938381503461596457e-01,
 281      -+4.998278695564493298213314152470e-01,
 282      -+5.092619017898079468040749192283e-01,
 283      -+5.186077642080456321529769963648e-01,
 284      -+5.278670896208423851138922177783e-01,
 285      -+5.370414658968836545667292441538e-01,
 286      -+5.461324375981356503823972092312e-01,
 287      -+5.551415075405015927154803595159e-01,
 288      -+5.640701382848029660713842900902e-01,
 289      -+5.729197535617855090927567266263e-01,
 290      -+5.816917396346224825206107537254e-01,
 291      -+5.903874466021763746419167081236e-01,
 292      -+5.990081896460833993816000244617e-01,
 293      -+6.075552502245417955010851527911e-01,
 294      -+6.160298772155140196475659281967e-01,
 295      -+6.244332880118935010425387440547e-01,
 296      -+6.327666695710378295457864685036e-01,
 297      -+6.410311794209312910556013344054e-01,
 298      -+6.492279466251098188908399699053e-01,
 299      -+6.573580727083600301418900232459e-01,
 300      -+6.654226325450904489500926100067e-01,
 301      -+6.734226752121667202979603888010e-01,
 302      -+6.813592248079030689480715595681e-01,
 303      -+6.892332812388089803249143378146e-01,
      238 +static const double T2[] = {            /* T2[j]=log(1+j/64+1/128) */
      239 +        +7.782140442054948947462900061137e-03,
      240 +        +2.316705928153437822879916096229e-02,
      241 +        +3.831886430213659919375532512380e-02,
      242 +        +5.324451451881228286587019378653e-02,
      243 +        +6.795066190850774939456527777263e-02,
      244 +        +8.244366921107459126816006866831e-02,
      245 +        +9.672962645855111229557105648746e-02,
      246 +        +1.108143663402901141948061693232e-01,
      247 +        +1.247034785009572358634065153809e-01,
      248 +        +1.384023228591191356853258736016e-01,
      249 +        +1.519160420258419750718034248969e-01,
      250 +        +1.652495728953071628756114492772e-01,
      251 +        +1.784076574728182971194002415109e-01,
      252 +        +1.913948529996294546092988075613e-01,
      253 +        +2.042155414286908915038203861962e-01,
      254 +        +2.168739383006143596190895257443e-01,
      255 +        +2.293741010648458299914807250461e-01,
      256 +        +2.417199368871451681443075159135e-01,
      257 +        +2.539152099809634441373232979066e-01,
      258 +        +2.659635484971379413391259265375e-01,
      259 +        +2.778684510034563061863500329234e-01,
      260 +        +2.896332925830426768788930555257e-01,
      261 +        +3.012613305781617810128755382338e-01,
      262 +        +3.127557100038968883862465596883e-01,
      263 +        +3.241194686542119760906707604350e-01,
      264 +        +3.353555419211378302571795798142e-01,
      265 +        +3.464667673462085809184621884258e-01,
      266 +        +3.574558889218037742260094901409e-01,
      267 +        +3.683255611587076530482301540504e-01,
      268 +        +3.790783529349694583908533456310e-01,
      269 +        +3.897167511400252133704636040035e-01,
      270 +        +4.002431641270127069293251019951e-01,
      271 +        +4.106599249852683859343062031758e-01,
      272 +        +4.209692946441296361288671615068e-01,
      273 +        +4.311734648183713408591724789556e-01,
      274 +        +4.412745608048752294894964416613e-01,
      275 +        +4.512746441394585851446923830790e-01,
      276 +        +4.611757151221701663679999255979e-01,
      277 +        +4.709797152187910125468978560564e-01,
      278 +        +4.806885293457519076766184554480e-01,
      279 +        +4.903039880451938381503461596457e-01,
      280 +        +4.998278695564493298213314152470e-01,
      281 +        +5.092619017898079468040749192283e-01,
      282 +        +5.186077642080456321529769963648e-01,
      283 +        +5.278670896208423851138922177783e-01,
      284 +        +5.370414658968836545667292441538e-01,
      285 +        +5.461324375981356503823972092312e-01,
      286 +        +5.551415075405015927154803595159e-01,
      287 +        +5.640701382848029660713842900902e-01,
      288 +        +5.729197535617855090927567266263e-01,
      289 +        +5.816917396346224825206107537254e-01,
      290 +        +5.903874466021763746419167081236e-01,
      291 +        +5.990081896460833993816000244617e-01,
      292 +        +6.075552502245417955010851527911e-01,
      293 +        +6.160298772155140196475659281967e-01,
      294 +        +6.244332880118935010425387440547e-01,
      295 +        +6.327666695710378295457864685036e-01,
      296 +        +6.410311794209312910556013344054e-01,
      297 +        +6.492279466251098188908399699053e-01,
      298 +        +6.573580727083600301418900232459e-01,
      299 +        +6.654226325450904489500926100067e-01,
      300 +        +6.734226752121667202979603888010e-01,
      301 +        +6.813592248079030689480715595681e-01,
      302 +        +6.892332812388089803249143378146e-01,
 304  303  };
 305      -/* INDENT ON */
 306  304  
 307  305  static double
 308      -large_gam(double x) {
      306 +large_gam(double x)
      307 +{
 309  308          double ss, zz, z, t1, t2, w, y, u;
 310  309          unsigned lx;
 311  310          int k, ix, j, m;
 312  311  
 313  312          ix = __HI(x);
 314  313          lx = __LO(x);
 315  314          m = (ix >> 20) - 0x3ff;                 /* exponent of x, range:3-5 */
 316  315          ix = (ix & 0x000fffff) | 0x3ff00000;    /* y = scale x to [1,2] */
 317  316          __HI(y) = ix;
 318  317          __LO(y) = lx;
 319  318          __HI(z) = (ix & 0xffffc000) | 0x2000;   /* z[j]=1+j/64+1/128 */
 320  319          __LO(z) = 0;
 321  320          j = (ix >> 14) & 0x3f;
 322  321          t1 = y + z;
 323  322          t2 = y - z;
 324  323          u = t2 / t1;
 325  324          ss = T1[m - 3] + T2[j] + u * (two + A1 * (u * u));
 326      -                                                        /* ss = log(x)-1 */
      325 +        /* ss = log(x)-1 */
      326 +
 327  327          /*
 328  328           * compute ww = (x-.5)*(log(x)-1) + .5*(log(2pi)-1) + 1/x*(P(1/x^2)))
 329  329           * where ss = log(x) - 1
 330  330           */
 331  331          z = one / x;
 332  332          zz = z * z;
 333  333          w = ((x - half) * ss + hln2pi) + z * (GP0 + zz * GP1 + (zz * zz) * GP2);
 334      -        k = (int) (w * invln2_32 + half);
      334 +        k = (int)(w * invln2_32 + half);
 335  335  
 336  336          /* compute the exponential of w */
 337  337          j = k & 0x1f;
 338  338          m = k >> 5;
 339      -        z = w - (double) k *ln2_32;
      339 +        z = w - (double)k * ln2_32;
 340  340          zz = S[j] * (one + z + (z * z) * (Et1 + z * Et2));
 341  341          __HI(zz) += m << 20;
 342  342          return (zz);
 343  343  }
 344      -/* INDENT OFF */
      344 +
      345 +
 345  346  /*
 346  347   * kpsin(x)= sin(pi*x)/pi
 347  348   *                 3        5        7        9
 348  349   *      = x+ks[0]*x +ks[1]*x +ks[2]*x +ks[3]*x
 349  350   */
 350  351  static const double ks[] = {
 351      --1.64493404985645811354476665052005342839447790544e+0000,
 352      -+8.11740794458351064092797249069438269367389272270e-0001,
 353      --1.90703144603551216933075809162889536878854055202e-0001,
 354      -+2.55742333994264563281155312271481108635575331201e-0002,
      352 +        -1.64493404985645811354476665052005342839447790544e+0000,
      353 +        +8.11740794458351064092797249069438269367389272270e-0001,
      354 +        -1.90703144603551216933075809162889536878854055202e-0001,
      355 +        +2.55742333994264563281155312271481108635575331201e-0002,
 355  356  };
 356      -/* INDENT ON */
 357  357  
 358  358  static double
 359      -kpsin(double x) {
      359 +kpsin(double x)
      360 +{
 360  361          double z;
 361  362  
 362  363          z = x * x;
 363  364          return (x + (x * z) * ((ks[0] + z * ks[1]) + (z * z) * (ks[2] + z *
 364      -                ks[3])));
      365 +            ks[3])));
 365  366  }
 366  367  
 367      -/* INDENT OFF */
      368 +
 368  369  /*
 369  370   * kpcos(x)= cos(pi*x)/pi
 370  371   *                     2        4        6
 371  372   *      = kc[0]+kc[1]*x +kc[2]*x +kc[3]*x
 372  373   */
 373  374  static const double kc[] = {
 374      -+3.18309886183790671537767526745028724068919291480e-0001,
 375      --1.57079581447762568199467875065854538626594937791e+0000,
 376      -+1.29183528092558692844073004029568674027807393862e+0000,
 377      --4.20232949771307685981015914425195471602739075537e-0001,
      375 +        +3.18309886183790671537767526745028724068919291480e-0001,
      376 +        -1.57079581447762568199467875065854538626594937791e+0000,
      377 +        +1.29183528092558692844073004029568674027807393862e+0000,
      378 +        -4.20232949771307685981015914425195471602739075537e-0001,
 378  379  };
 379      -/* INDENT ON */
 380  380  
 381  381  static double
 382      -kpcos(double x) {
      382 +kpcos(double x)
      383 +{
 383  384          double z;
 384  385  
 385  386          z = x * x;
 386  387          return (kc[0] + z * (kc[1] + z * kc[2] + (z * z) * kc[3]));
 387  388  }
 388  389  
 389      -/* INDENT OFF */
 390      -static const double
 391      -t0z1 = 0.134861805732790769689793935774652917006,
 392      -t0z2 = 0.461632144968362341262659542325721328468,
 393      -t0z3 = 0.819773101100500601787868704921606996312;
 394      -        /* 1.134861805732790769689793935774652917006 */
 395      -/* INDENT ON */
      390 +static const double t0z1 = 0.134861805732790769689793935774652917006,
      391 +        t0z2 = 0.461632144968362341262659542325721328468,
      392 +        t0z3 = 0.819773101100500601787868704921606996312;
      393 +
      394 +/*
      395 + * 1.134861805732790769689793935774652917006
      396 + */
 396  397  
 397  398  /*
 398  399   * gamma(x+i) for 0 <= x < 1
 399  400   */
 400  401  static double
 401      -gam_n(int i, double x) {
      402 +gam_n(int i, double x)
      403 +{
 402  404          double rr = 0.0L, yy;
 403  405          double z1, z2;
 404  406  
 405  407          /* compute yy = gamma(x+1) */
 406  408          if (x > 0.2845) {
 407  409                  if (x > 0.6374)
 408  410                          yy = GT3(x - t0z3);
 409  411                  else
 410  412                          yy = GT2(x - t0z2);
 411      -        } else
      413 +        } else {
 412  414                  yy = GT1(x - t0z1);
      415 +        }
 413  416  
 414  417          /* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */
 415  418          switch (i) {
 416      -        case 0:         /* yy/x */
      419 +        case 0:                         /* yy/x */
 417  420                  rr = yy / x;
 418  421                  break;
 419      -        case 1:         /* yy */
      422 +        case 1:                         /* yy */
 420  423                  rr = yy;
 421  424                  break;
 422      -        case 2:         /* (x+1)*yy */
      425 +        case 2:                         /* (x+1)*yy */
 423  426                  rr = (x + one) * yy;
 424  427                  break;
 425      -        case 3:         /* (x+2)*(x+1)*yy */
      428 +        case 3:                         /* (x+2)*(x+1)*yy */
 426  429                  rr = (x + one) * (x + two) * yy;
 427  430                  break;
 428  431  
 429      -        case 4:         /* (x+1)*(x+3)*(x+2)*yy */
      432 +        case 4:                         /* (x+1)*(x+3)*(x+2)*yy */
 430  433                  rr = (x + one) * (x + two) * ((x + 3.0) * yy);
 431  434                  break;
 432      -        case 5:         /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */
      435 +        case 5:                         /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */
 433  436                  z1 = (x + two) * (x + 3.0) * yy;
 434  437                  z2 = (x + one) * (x + 4.0);
 435  438                  rr = z1 * z2;
 436  439                  break;
 437      -        case 6:         /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
      440 +        case 6:                         /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
 438  441                  z1 = (x + two) * (x + 3.0);
 439  442                  z2 = (x + 5.0) * yy;
 440  443                  rr = z1 * (z1 - two) * z2;
 441  444                  break;
 442      -        case 7:         /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
      445 +        case 7:                 /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
 443  446                  z1 = (x + two) * (x + 3.0);
 444  447                  z2 = (x + 5.0) * (x + 6.0) * yy;
 445  448                  rr = z1 * (z1 - two) * z2;
 446  449                  break;
 447  450          }
      451 +
 448  452          return (rr);
 449  453  }
 450  454  
 451  455  float
 452      -tgammaf(float xf) {
      456 +tgammaf(float xf)
      457 +{
 453  458          float zf;
 454  459          double ss, ww;
 455  460          double x, y, z;
 456  461          int i, j, k, ix, hx, xk;
 457  462  
 458      -        hx = *(int *) &xf;
      463 +        hx = *(int *)&xf;
 459  464          ix = hx & 0x7fffffff;
 460  465  
 461      -        x = (double) xf;
      466 +        x = (double)xf;
      467 +
 462  468          if (ix < 0x33800000)
 463      -                return (1.0F / xf);     /* |x| < 2**-24 */
      469 +                return (1.0F / xf);                     /* |x| < 2**-24 */
 464  470  
 465  471          if (ix >= 0x7f800000)
 466      -                return (xf * ((hx < 0)? 0.0F : xf)); /* +-Inf or NaN */
      472 +                return (xf * ((hx < 0) ? 0.0F : xf));   /* +-Inf or NaN */
 467  473  
 468      -        if (hx > 0x420C290F)    /* x > 35.040096283... overflow */
 469      -                return (float)(x / tiny);
      474 +        if (hx > 0x420C290F)    /* x > 35.040096283... overflow */
      475 +                return ((float)(x / tiny));
 470  476  
 471      -        if (hx >= 0x41000000)   /* x >= 8 */
 472      -                return ((float) large_gam(x));
      477 +        if (hx >= 0x41000000)                           /* x >= 8 */
      478 +                return ((float)large_gam(x));
 473  479  
 474      -        if (hx > 0) {           /* 0 < x < 8 */
 475      -                i = (int) xf;
 476      -                return ((float) gam_n(i, x - (double) i));
      480 +        if (hx > 0) {                                   /* 0 < x < 8 */
      481 +                i = (int)xf;
      482 +                return ((float)gam_n(i, x - (double)i));
 477  483          }
 478  484  
 479      -        /* negative x */
 480      -        /* INDENT OFF */
      485 +        /*
      486 +         * negative x
      487 +         */
      488 +
 481  489          /*
 482  490           * compute xk =
 483  491           *      -2 ... x is an even int (-inf is considered even)
 484  492           *      -1 ... x is an odd int
 485  493           *      +0 ... x is not an int but chopped to an even int
 486  494           *      +1 ... x is not an int but chopped to an odd int
 487  495           */
 488      -        /* INDENT ON */
 489  496          xk = 0;
      497 +
 490  498          if (ix >= 0x4b000000) {
 491  499                  if (ix > 0x4b000000)
 492  500                          xk = -2;
 493  501                  else
 494  502                          xk = -2 + (ix & 1);
 495  503          } else if (ix >= 0x3f800000) {
 496  504                  k = (ix >> 23) - 0x7f;
 497  505                  j = ix >> (23 - k);
      506 +
 498  507                  if ((j << (23 - k)) == ix)
 499  508                          xk = -2 + (j & 1);
 500  509                  else
 501  510                          xk = j & 1;
 502  511          }
      512 +
 503  513          if (xk < 0) {
 504  514                  /* 0/0 invalid NaN, ideally gamma(-n)= (-1)**(n+1) * inf */
 505  515                  zf = xf - xf;
 506  516                  return (zf / zf);
 507  517          }
 508  518  
 509  519          /* negative underflow thresold */
 510      -        if (ix > 0x4224000B) {  /* x < -(41+11ulp) */
      520 +        if (ix > 0x4224000B) {          /* x < -(41+11ulp) */
 511  521                  if (xk == 0)
 512  522                          z = -tiny;
 513  523                  else
 514  524                          z = tiny;
      525 +
 515  526                  return ((float)z);
 516  527          }
 517  528  
 518      -        /* INDENT OFF */
 519      -        /* now compute gamma(x) by  -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x */
      529 +        /*
      530 +         * now compute gamma(x) by  -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x
      531 +         */
      532 +
 520  533          /*
 521  534           * First compute ss = -sin(pi*y)/pi , so that
 522  535           * gamma(x) = 1/(ss*gamma(1+y))
 523  536           */
 524      -        /* INDENT ON */
 525  537          y = -x;
 526      -        j = (int) y;
 527      -        z = y - (double) j;
 528      -        if (z > 0.3183098861837906715377675)
      538 +        j = (int)y;
      539 +        z = y - (double)j;
      540 +
      541 +        if (z > 0.3183098861837906715377675) {
 529  542                  if (z > 0.6816901138162093284622325)
 530  543                          ss = kpsin(one - z);
 531  544                  else
 532  545                          ss = kpcos(0.5 - z);
 533      -        else
      546 +        } else {
 534  547                  ss = kpsin(z);
      548 +        }
      549 +
 535  550          if (xk == 0)
 536  551                  ss = -ss;
 537  552  
 538  553          /* Then compute ww = gamma(1+y)  */
 539  554          if (j < 7)
 540  555                  ww = gam_n(j + 1, z);
 541  556          else
 542  557                  ww = large_gam(y + one);
 543  558  
 544  559          /* return 1/(ss*ww) */
 545      -        return ((float) (one / (ww * ss)));
      560 +        return ((float)(one / (ww * ss)));
 546  561  }
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX