lib/gamln.c

   1 #include "xmath.h"
   2
   3 /*
   4 log gamma function from Numerical Recipes
   5 */
   6
   7 static double cof[6] = {
   8    76.18009173,
   9   -86.50532033,
  10    24.01409822,
  11    -1.231739516,
  12     0.120858003e-2,
  13    -0.536382e-5
  14 };
  15
  16 double gamma(xx)
  17         double xx;
  18 {
  19   double x, tmp, ser;
  20   int j;
  21
  22   if (xx < 1.0) return(gamma(1.0 + xx) - log(xx));
  23   else {
  24     x = xx - 1.0;
  25     tmp = x + 5.5;
  26     tmp -= (x + 0.5) * log(tmp);
  27     ser = 1.0;
  28     for (j = 0; j < 6; j++) {
  29       x += 1.0;
  30       ser += cof[j] / x;
  31     }
  32     return(-tmp + log(2.50662827465 * ser));
  33   }
  34 }
  35