src/math/tgamma.c

   1 /*
   2 "A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964)
   3 "Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001)
   4 "An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004)
   5
   6 approximation method:
   7
   8                         (x - 0.5)         S(x)
   9 Gamma(x) = (x + g - 0.5)         *  ----------------
  10                                     exp(x + g - 0.5)
  11
  12 with
  13                  a1      a2      a3            aN
  14 S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ]
  15                x + 1   x + 2   x + 3         x + N
  16
  17 with a0, a1, a2, a3,.. aN constants which depend on g.
  18
  19 for x < 0 the following reflection formula is used:
  20
  21 Gamma(x)*Gamma(-x) = -pi/(x sin(pi x))
  22
  23 most ideas and constants are from boost and python
  24 */
  25 #include "libm.h"
  26
  27 static const double pi = 3.141592653589793238462643383279502884;
  28
  29 /* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */
  30 static double sinpi(double x)
  31 {
  32         int n;
  33
  34         /* argument reduction: x = |x| mod 2 */
  35         /* spurious inexact when x is odd int */
  36         x = x * 0.5;
  37         x = 2 * (x - floor(x));
  38
  39         /* reduce x into [-.25,.25] */
  40         n = 4 * x;
  41         n = (n+1)/2;
  42         x -= n * 0.5;
  43
  44         x *= pi;
  45         switch (n) {
  46         default: /* case 4 */
  47         case 0:
  48                 return __sin(x, 0, 0);
  49         case 1:
  50                 return __cos(x, 0);
  51         case 2:
  52                 return __sin(-x, 0, 0);
  53         case 3:
  54                 return -__cos(x, 0);
  55         }
  56 }
  57
  58 #define N 12
  59 //static const double g = 6.024680040776729583740234375;
  60 static const double gmhalf = 5.524680040776729583740234375;
  61 static const double Snum[N+1] = {
  62         23531376880.410759688572007674451636754734846804940,
  63         42919803642.649098768957899047001988850926355848959,
  64         35711959237.355668049440185451547166705960488635843,
  65         17921034426.037209699919755754458931112671403265390,
  66         6039542586.3520280050642916443072979210699388420708,
  67         1439720407.3117216736632230727949123939715485786772,
  68         248874557.86205415651146038641322942321632125127801,
  69         31426415.585400194380614231628318205362874684987640,
  70         2876370.6289353724412254090516208496135991145378768,
  71         186056.26539522349504029498971604569928220784236328,
  72         8071.6720023658162106380029022722506138218516325024,
  73         210.82427775157934587250973392071336271166969580291,
  74         2.5066282746310002701649081771338373386264310793408,
  75 };
  76 static const double Sden[N+1] = {
  77         0, 39916800, 120543840, 150917976, 105258076, 45995730, 13339535,
  78         2637558, 357423, 32670, 1925, 66, 1,
  79 };
  80 /* n! for small integer n */
  81 static const double fact[] = {
  82         1, 1, 2, 6, 24, 120, 720, 5040.0, 40320.0, 362880.0, 3628800.0, 39916800.0,
  83         479001600.0, 6227020800.0, 87178291200.0, 1307674368000.0, 20922789888000.0,
  84         355687428096000.0, 6402373705728000.0, 121645100408832000.0,
  85         2432902008176640000.0, 51090942171709440000.0, 1124000727777607680000.0,
  86 };
  87
  88 /* S(x) rational function for positive x */
  89 static double S(double x)
  90 {
  91         double_t num = 0, den = 0;
  92         int i;
  93
  94         /* to avoid overflow handle large x differently */
  95         if (x < 8)
  96                 for (i = N; i >= 0; i--) {
  97                         num = num * x + Snum[i];
  98                         den = den * x + Sden[i];
  99                 }
 100         else
 101                 for (i = 0; i <= N; i++) {
 102                         num = num / x + Snum[i];
 103                         den = den / x + Sden[i];
 104                 }
 105         return num/den;
 106 }
 107
 108 double tgamma(double x)
 109 {
 110         union {double f; uint64_t i;} u = {x};
 111         double absx, y;
 112         double_t dy, z, r;
 113         uint32_t ix = u.i>>32 & 0x7fffffff;
 114         int sign = u.i>>63;
 115
 116         /* special cases */
 117         if (ix >= 0x7ff00000)
 118                 /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */
 119                 return x + INFINITY;
 120         if (ix < (0x3ff-54)<<20)
 121                 /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */
 122                 return 1/x;
 123
 124         /* integer arguments */
 125         /* raise inexact when non-integer */
 126         if (x == floor(x)) {
 127                 if (sign)
 128                         return 0/0.0;
 129                 if (x <= sizeof fact/sizeof *fact)
 130                         return fact[(int)x - 1];
 131         }
 132
 133         /* x >= 172: tgamma(x)=inf with overflow */
 134         /* x =< -184: tgamma(x)=+-0 with underflow */
 135         if (ix >= 0x40670000) { /* |x| >= 184 */
 136                 if (sign) {
 137                         FORCE_EVAL((float)(0x1p-126/x));
 138                         if (floor(x) * 0.5 == floor(x * 0.5))
 139                                 return 0;
 140                         return -0.0;
 141                 }
 142                 x *= 0x1p1023;
 143                 return x;
 144         }
 145
 146         absx = sign ? -x : x;
 147
 148         /* handle the error of x + g - 0.5 */
 149         y = absx + gmhalf;
 150         if (absx > gmhalf) {
 151                 dy = y - absx;
 152                 dy -= gmhalf;
 153         } else {
 154                 dy = y - gmhalf;
 155                 dy -= absx;
 156         }
 157
 158         z = absx - 0.5;
 159         r = S(absx) * exp(-y);
 160         if (x < 0) {
 161                 /* reflection formula for negative x */
 162                 /* sinpi(absx) is not 0, integers are already handled */
 163                 r = -pi / (sinpi(absx) * absx * r);
 164                 dy = -dy;
 165                 z = -z;
 166         }
 167         r += dy * (gmhalf+0.5) * r / y;
 168         z = pow(y, 0.5*z);
 169         y = r * z * z;
 170         return y;
 171 }
 172
 173 #if 0
 174 double __lgamma_r(double x, int *sign)
 175 {
 176         double r, absx;
 177
 178         *sign = 1;
 179
 180         /* special cases */
 181         if (!isfinite(x))
 182                 /* lgamma(nan)=nan, lgamma(+-inf)=inf */
 183                 return x*x;
 184
 185         /* integer arguments */
 186         if (x == floor(x) && x <= 2) {
 187                 /* n <= 0: lgamma(n)=inf with divbyzero */
 188                 /* n == 1,2: lgamma(n)=0 */
 189                 if (x <= 0)
 190                         return 1/0.0;
 191                 return 0;
 192         }
 193
 194         absx = fabs(x);
 195
 196         /* lgamma(x) ~ -log(|x|) for tiny |x| */
 197         if (absx < 0x1p-54) {
 198                 *sign = 1 - 2*!!signbit(x);
 199                 return -log(absx);
 200         }
 201
 202         /* use tgamma for smaller |x| */
 203         if (absx < 128) {
 204                 x = tgamma(x);
 205                 *sign = 1 - 2*!!signbit(x);
 206                 return log(fabs(x));
 207         }
 208
 209         /* second term (log(S)-g) could be more precise here.. */
 210         /* or with stirling: (|x|-0.5)*(log(|x|)-1) + poly(1/|x|) */
 211         r = (absx-0.5)*(log(absx+gmhalf)-1) + (log(S(absx)) - (gmhalf+0.5));
 212         if (x < 0) {
 213                 /* reflection formula for negative x */
 214                 x = sinpi(absx);
 215                 *sign = 2*!!signbit(x) - 1;
 216                 r = log(pi/(fabs(x)*absx)) - r;
 217         }
 218         return r;
 219 }
 220
 221 weak_alias(__lgamma_r, lgamma_r);
 222 #endif