Quelle  gamma.c

/* gamma.c
 *
 * Gamma function
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, gamma();
 * extern int sgngam;
 *
 * y = gamma( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns gamma function of the argument.  The result is
 * correctly signed, and the sign (+1 or -1) is also
 * returned in a global (extern) variable named sgngam.
 * This variable is also filled in by the logarithmic gamma
 * function lgam().
 *
 * Arguments |x| <= 34 are reduced by recurrence and the function
 * approximated by a rational function of degree 6/7 in the
 * interval (2,3).  Large arguments are handled by Stirling's
 * formula. Large negative arguments are made positive using
 * a reflection formula.
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC      -34, 34      10000       1.3e-16     2.5e-17
 *    IEEE    -170,-33      20000       2.3e-15     3.3e-16
 *    IEEE     -33,  33     20000       9.4e-16     2.2e-16
 *    IEEE      33, 171.6   20000       2.3e-15     3.2e-16
 *
 * Error for arguments outside the test range will be larger
 * owing to error amplification by the exponential function.
 *
 */
/* lgam()
 *
 * Natural logarithm of gamma function
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, lgam();
 * extern int sgngam;
 *
 * y = lgam( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the base e (2.718...) logarithm of the absolute
 * value of the gamma function of the argument.
 * The sign (+1 or -1) of the gamma function is returned in a
 * global (extern) variable named sgngam.
 *
 * For arguments greater than 13, the logarithm of the gamma
 * function is approximated by the logarithmic version of
 * Stirling's formula using a polynomial approximation of
 * degree 4. Arguments between -33 and +33 are reduced by
 * recurrence to the interval [2,3] of a rational approximation.
 * The cosecant reflection formula is employed for arguments
 * less than -33.
 *
 * Arguments greater than MAXLGM return MAXNUM and an error
 * message.  MAXLGM = 2.035093e36 for DEC
 * arithmetic or 2.556348e305 for IEEE arithmetic.
 *
 *
 *
 * ACCURACY:
 *
 *
 * arithmetic      domain        # trials     peak         rms
 *    DEC     0, 3                  7000     5.2e-17     1.3e-17
 *    DEC     2.718, 2.035e36       5000     3.9e-17     9.9e-18
 *    IEEE    0, 3                 28000     5.4e-16     1.1e-16
 *    IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
 * The error criterion was relative when the function magnitude
 * was greater than one but absolute when it was less than one.
 *
 * The following test used the relative error criterion, though
 * at certain points the relative error could be much higher than
 * indicated.
 *    IEEE    -200, -4             10000     4.8e-16     1.3e-16
 *
 */

/* gamma.c */
/* gamma function */

/*
Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
*/

#include "mconf.h"

#ifdef UNK
static double P[] = {1.60119522476751861407E-4, 1.19135147006586384913E-3,
                     1.04213797561761569935E-2, 4.76367800457137231464E-2,
                     2.07448227648435975150E-1, 4.94214826801497100753E-1,
                     9.99999999999999996796E-1};
static double Q[] = {-2.31581873324120129819E-5, 5.39605580493303397842E-4,
                     -4.45641913851797240494E-3, 1.18139785222060435552E-2,
                     3.58236398605498653373E-2,  -2.34591795718243348568E-1,
                     7.14304917030273074085E-2,  1.00000000000000000320E0};
#define MAXGAM 171.624376956302725
static double LOGPI = 1.14472988584940017414;
#endif

#ifdef DEC
static unsigned short P[] = {
    0035047, 0162701, 0146301, 0005234, 0035634, 0023437, 0032065,
    0176530, 0036452, 0137157, 0047330, 0122574, 0037103, 0017310,
    0143041, 0017232, 0037524, 0066516, 0162563, 0164605, 0037775,
    0004671, 0146237, 0014222, 0040200, 0000000, 0000000, 0000000};
static unsigned short Q[] = {
    0134302, 0041724, 0020006, 0116565, 0035415, 0072121, 0044251, 0025634,
    0136222, 0003447, 0035205, 0121114, 0036501, 0107552, 0154335, 0104271,
    0037022, 0135717, 0014776, 0171471, 0137560, 0034324, 0165024, 0037021,
    0037222, 0045046, 0047151, 0161213, 0040200, 0000000, 0000000, 0000000};
#define MAXGAM 34.84425627277176174
static unsigned short LPI[4] = {
    0040222,
    0103202,
    0043475,
    0006750,
};
#define LOGPI *(double *)LPI
#endif

#ifdef IBMPC
static unsigned short P[] = {0x2153, 0x3998, 0xfcb8, 0x3f24, 0xbfab, 0xe686,
                             0x84e3, 0x3f53, 0x14b0, 0xe9db, 0x57cd, 0x3f85,
                             0x23d3, 0x18c4, 0x63d9, 0x3fa8, 0x7d31, 0xdcae,
                             0x8da9, 0x3fca, 0xe312, 0x3993, 0xa137, 0x3fdf,
                             0x0000, 0x0000, 0x0000, 0x3ff0};
static unsigned short Q[] = {
    0xd3af, 0x8400, 0x487a, 0xbef8, 0x2573, 0x2915, 0xae8a, 0x3f41,
    0xb44a, 0xe750, 0x40e4, 0xbf72, 0xb117, 0x5b1b, 0x31ed, 0x3f88,
    0xde67, 0xe33f, 0x5779, 0x3fa2, 0x87c2, 0x9d42, 0x071a, 0xbfce,
    0x3c51, 0xc9cd, 0x4944, 0x3fb2, 0x0000, 0x0000, 0x0000, 0x3ff0};
#define MAXGAM 171.624376956302725
static unsigned short LPI[4] = {
    0xa1bd,
    0x48e7,
    0x50d0,
    0x3ff2,
};
#define LOGPI *(double *)LPI
#endif

#ifdef MIEEE
static unsigned short P[] = {0x3f24, 0xfcb8, 0x3998, 0x2153, 0x3f53, 0x84e3,
                             0xe686, 0xbfab, 0x3f85, 0x57cd, 0xe9db, 0x14b0,
                             0x3fa8, 0x63d9, 0x18c4, 0x23d3, 0x3fca, 0x8da9,
                             0xdcae, 0x7d31, 0x3fdf, 0xa137, 0x3993, 0xe312,
                             0x3ff0, 0x0000, 0x0000, 0x0000};
static unsigned short Q[] = {
    0xbef8, 0x487a, 0x8400, 0xd3af, 0x3f41, 0xae8a, 0x2915, 0x2573,
    0xbf72, 0x40e4, 0xe750, 0xb44a, 0x3f88, 0x31ed, 0x5b1b, 0xb117,
    0x3fa2, 0x5779, 0xe33f, 0xde67, 0xbfce, 0x071a, 0x9d42, 0x87c2,
    0x3fb2, 0x4944, 0xc9cd, 0x3c51, 0x3ff0, 0x0000, 0x0000, 0x0000};
#define MAXGAM 171.624376956302725
static unsigned short LPI[4] = {
    0x3ff2,
    0x50d0,
    0x48e7,
    0xa1bd,
};
#define LOGPI *(double *)LPI
#endif

/* Stirling's formula for the gamma function */
#if UNK
static double STIR[5] = {
    7.87311395793093628397E-4,  -2.29549961613378126380E-4,
    -2.68132617805781232825E-3, 3.47222221605458667310E-3,
    8.33333333333482257126E-2,
};
#define MAXSTIR 143.01608
static double SQTPI = 2.50662827463100050242E0;
#endif
#if DEC
static unsigned short STIR[20] = {
    0035516, 0061622, 0144553, 0112224, 0135160, 0131531, 0037460,
    0165740, 0136057, 0134460, 0037242, 0077270, 0036143, 0107070,
    0156306, 0027751, 0037252, 0125252, 0125252, 0146064,
};
#define MAXSTIR 26.77
static unsigned short SQT[4] = {
    0040440,
    0066230,
    0177661,
    0034055,
};
#define SQTPI *(double *)SQT
#endif
#if IBMPC
static unsigned short STIR[20] = {
    0x7293, 0x592d, 0xcc72, 0x3f49, 0x1d7c, 0x27e6, 0x166b,
    0xbf2e, 0x4fd7, 0x07d4, 0xf726, 0xbf65, 0xc5fd, 0x1b98,
    0x71c7, 0x3f6c, 0x5986, 0x5555, 0x5555, 0x3fb5,
};
#define MAXSTIR 143.01608
static unsigned short SQT[4] = {
    0x2706,
    0x1ff6,
    0x0d93,
    0x4004,
};
#define SQTPI *(double *)SQT
#endif
#if MIEEE
static unsigned short STIR[20] = {
    0x3f49, 0xcc72, 0x592d, 0x7293, 0xbf2e, 0x166b, 0x27e6,
    0x1d7c, 0xbf65, 0xf726, 0x07d4, 0x4fd7, 0x3f6c, 0x71c7,
    0x1b98, 0xc5fd, 0x3fb5, 0x5555, 0x5555, 0x5986,
};
#define MAXSTIR 143.01608
static unsigned short SQT[4] = {
    0x4004,
    0x0d93,
    0x1ff6,
    0x2706,
};
#define SQTPI *(double *)SQT
#endif

int sgngam = 0;
extern int sgngam;
extern double MAXLOG, MAXNUM, PI;
#ifdef ANSIPROT
extern double pow(double, double);
extern double log(double);
extern double exp(double);
extern double sin(double);
extern double polevl(double, void *, int);
extern double p1evl(double, void *, int);
extern double floor(double);
extern double fabs(double);
extern int isnan(double);
extern int isfinite(double);
static double stirf(double);
double lgam(double);
#else
double pow(), log(), exp(), sin(), polevl(), p1evl(), floor(), fabs();
int isnan(), isfinite();
static double stirf();
double lgam();
#endif
#ifdef INFINITIES
extern double INFINITY;
#endif
#ifdef NANS
extern double NAN;
#endif

/* Gamma function computed by Stirling's formula.
 * The polynomial STIR is valid for 33 <= x <= 172.
 */
static double stirf(x) double x;
{
  double y, w, v;

  w = 1.0 / x;
  w = 1.0 + w * polevl(w, STIR, 4);
  y = exp(x);
  if (x > MAXSTIR) { /* Avoid overflow in pow() */
    v = pow(x, 0.5 * x - 0.25);
    y = v * (v / y);
  } else {
    y = pow(x, x - 0.5) / y;
  }
  y = SQTPI * y * w;
  return (y);
}

double gamma(x) double x;
{
  double p, q, z;
  int i;

  sgngam = 1;
#ifdef NANS
  if (isnan(x))
    return (x);
#endif
#ifdef INFINITIES
#ifdef NANS
  if (x == INFINITY)
    return (x);
  if (x == -INFINITY)
    return (NAN);
#else
  if (!isfinite(x))
    return (x);
#endif
#endif
  q = fabs(x);

  if (q > 33.0) {
    if (x < 0.0) {
      p = floor(q);
      if (p == q) {
#ifdef NANS
      gamnan:
        mtherr("gamma", DOMAIN);
        return (NAN);
#else
        goto goverf;
#endif
      }
      i = p;
      if ((i & 1) == 0)
        sgngam = -1;
      z = q - p;
      if (z > 0.5) {
        p += 1.0;
        z = q - p;
      }
      z = q * sin(PI * z);
      if (z == 0.0) {
#ifdef INFINITIES
        return (sgngam * INFINITY);
#else
      goverf:
        mtherr("gamma", OVERFLOW);
        return (sgngam * MAXNUM);
#endif
      }
      z = fabs(z);
      z = PI / (z * stirf(q));
    } else {
      z = stirf(x);
    }
    return (sgngam * z);
  }

  z = 1.0;
  while (x >= 3.0) {
    x -= 1.0;
    z *= x;
  }

  while (x < 0.0) {
    if (x > -1.E-9)
      goto small;
    z /= x;
    x += 1.0;
  }

  while (x < 2.0) {
    if (x < 1.e-9)
      goto small;
    z /= x;
    x += 1.0;
  }

  if (x == 2.0)
    return (z);

  x -= 2.0;
  p = polevl(x, P, 6);
  q = polevl(x, Q, 7);
  return (z * p / q);

small:
  if (x == 0.0) {
#ifdef INFINITIES
#ifdef NANS
    goto gamnan;
#else
    return (INFINITY);
#endif
#else
    mtherr("gamma", SING);
    return (MAXNUM);
#endif
  } else
    return (z / ((1.0 + 0.5772156649015329 * x) * x));
}

/* A[]: Stirling's formula expansion of log gamma
 * B[], C[]: log gamma function between 2 and 3
 */
#ifdef UNK
static double A[] = {8.11614167470508450300E-4, -5.95061904284301438324E-4,
                     7.93650340457716943945E-4, -2.77777777730099687205E-3,
                     8.33333333333331927722E-2};
static double B[] = {-1.37825152569120859100E3, -3.88016315134637840924E4,
                     -3.31612992738871184744E5, -1.16237097492762307383E6,
                     -1.72173700820839662146E6, -8.53555664245765465627E5};
static double C[] = {
    /* 1.00000000000000000000E0, */
    -3.51815701436523470549E2, -1.70642106651881159223E4,
    -2.20528590553854454839E5, -1.13933444367982507207E6,
    -2.53252307177582951285E6, -2.01889141433532773231E6};
/* log( sqrt( 2*pi ) ) */
static double LS2PI = 0.91893853320467274178;
#define MAXLGM 2.556348e305
#endif

#ifdef DEC
static unsigned short A[] = {0035524, 0141201, 0034633, 0031405, 0135433,
                             0176755, 0126007, 0045030, 0035520, 0006371,
                             0003342, 0172730, 0136066, 0005540, 0132605,
                             0026407, 0037252, 0125252, 0125252, 0125132};
static unsigned short B[] = {
    0142654, 0044014, 0077633, 0035410, 0144027, 0110641, 0125335, 0144760,
    0144641, 0165637, 0142204, 0047447, 0145215, 0162027, 0146246, 0155211,
    0145322, 0026110, 0010317, 0110130, 0145120, 0061472, 0120300, 0025363};
static unsigned short C[] = {
    /*0040200,0000000,0000000,0000000*/
    0142257, 0164150, 0163630, 0112622, 0143605, 0050153, 0156116, 0135272,
    0144527, 0056045, 0145642, 0062332, 0145213, 0012063, 0106250, 0001025,
    0145432, 0111254, 0044577, 0115142, 0145366, 0071133, 0050217, 0005122};
/* log( sqrt( 2*pi ) ) */
static unsigned short LS2P[] = {
    040153,
    037616,
    041445,
    0172645,
};
#define LS2PI *(double *)LS2P
#define MAXLGM 2.035093e36
#endif

#ifdef IBMPC
static unsigned short A[] = {0x6661, 0x2733, 0x9850, 0x3f4a, 0xe943,
                             0xb580, 0x7fbd, 0xbf43, 0x5ebb, 0x20dc,
                             0x019f, 0x3f4a, 0xa5a1, 0x16b0, 0xc16c,
                             0xbf66, 0x554b, 0x5555, 0x5555, 0x3fb5};
static unsigned short B[] = {0x6761, 0x8ff3, 0x8901, 0xc095, 0xb93e, 0x355b,
                             0xf234, 0xc0e2, 0x89e5, 0xf890, 0x3d73, 0xc114,
                             0xdb51, 0xf994, 0xbc82, 0xc131, 0xf20b, 0x0219,
                             0x4589, 0xc13a, 0x055e, 0x5418, 0x0c67, 0xc12a};
static unsigned short C[] = {
    /*0x0000,0x0000,0x0000,0x3ff0,*/
    0x12b2, 0x1cf3, 0xfd0d, 0xc075, 0xd757, 0x7b89, 0xaa0d, 0xc0d0,
    0x4c9b, 0xb974, 0xeb84, 0xc10a, 0x0043, 0x7195, 0x6286, 0xc131,
    0xf34c, 0x892f, 0x5255, 0xc143, 0xe14a, 0x6a11, 0xce4b, 0xc13e};
/* log( sqrt( 2*pi ) ) */
static unsigned short LS2P[] = {0xbeb5, 0xc864, 0x67f1, 0x3fed};
#define LS2PI *(double *)LS2P
#define MAXLGM 2.556348e305
#endif

#ifdef MIEEE
static unsigned short A[] = {0x3f4a, 0x9850, 0x2733, 0x6661, 0xbf43,
                             0x7fbd, 0xb580, 0xe943, 0x3f4a, 0x019f,
                             0x20dc, 0x5ebb, 0xbf66, 0xc16c, 0x16b0,
                             0xa5a1, 0x3fb5, 0x5555, 0x5555, 0x554b};
static unsigned short B[] = {0xc095, 0x8901, 0x8ff3, 0x6761, 0xc0e2, 0xf234,
                             0x355b, 0xb93e, 0xc114, 0x3d73, 0xf890, 0x89e5,
                             0xc131, 0xbc82, 0xf994, 0xdb51, 0xc13a, 0x4589,
                             0x0219, 0xf20b, 0xc12a, 0x0c67, 0x5418, 0x055e};
static unsigned short C[] = {0xc075, 0xfd0d, 0x1cf3, 0x12b2, 0xc0d0, 0xaa0d,
                             0x7b89, 0xd757, 0xc10a, 0xeb84, 0xb974, 0x4c9b,
                             0xc131, 0x6286, 0x7195, 0x0043, 0xc143, 0x5255,
                             0x892f, 0xf34c, 0xc13e, 0xce4b, 0x6a11, 0xe14a};
/* log( sqrt( 2*pi ) ) */
static unsigned short LS2P[] = {0x3fed, 0x67f1, 0xc864, 0xbeb5};
#define LS2PI *(double *)LS2P
#define MAXLGM 2.556348e305
#endif

/* Logarithm of gamma function */

double lgam(x) double x;
{
  double p, q, u, w, z;
  int i;

  sgngam = 1;
#ifdef NANS
  if (isnan(x))
    return (x);
#endif

#ifdef INFINITIES
  if (!isfinite(x))
    return (INFINITY);
#endif

  if (x < -34.0) {
    q = -x;
    w = lgam(q); /* note this modifies sgngam! */
    p = floor(q);
    if (p == q) {
    lgsing:
#ifdef INFINITIES
      mtherr("lgam", SING);
      return (INFINITY);
#else
      goto loverf;
#endif
    }
    i = p;
    if ((i & 1) == 0)
      sgngam = -1;
    else
      sgngam = 1;
    z = q - p;
    if (z > 0.5) {
      p += 1.0;
      z = p - q;
    }
    z = q * sin(PI * z);
    if (z == 0.0)
      goto lgsing;
    /* z = log(PI) - log( z ) - w;*/
    z = LOGPI - log(z) - w;
    return (z);
  }

  if (x < 13.0) {
    z = 1.0;
    p = 0.0;
    u = x;
    while (u >= 3.0) {
      p -= 1.0;
      u = x + p;
      z *= u;
    }
    while (u < 2.0) {
      if (u == 0.0)
        goto lgsing;
      z /= u;
      p += 1.0;
      u = x + p;
    }
    if (z < 0.0) {
      sgngam = -1;
      z = -z;
    } else
      sgngam = 1;
    if (u == 2.0)
      return (log(z));
    p -= 2.0;
    x = x + p;
    p = x * polevl(x, B, 5) / p1evl(x, C, 6);
    return (log(z) + p);
  }

  if (x > MAXLGM) {
#ifdef INFINITIES
    return (sgngam * INFINITY);
#else
  loverf:
    mtherr("lgam", OVERFLOW);
    return (sgngam * MAXNUM);
#endif
  }

  q = (x - 0.5) * log(x) - x + LS2PI;
  if (x > 1.0e8)
    return (q);

  p = 1.0 / (x * x);
  if (x >= 1000.0)
    q += ((7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3) * p +
          0.0833333333333333333333) /
         x;
  else
    q += polevl(p, A, 4) / x;
  return (q);
}