sys/src/ape/lib/ap/math/gamma.c

   1 #include <math.h>
   2 #include <errno.h>
   3
   4 /*
   5         C program for floating point log gamma function
   6
   7         gamma(x) computes the log of the absolute
   8         value of the gamma function.
   9         The sign of the gamma function is returned in the
  10         external quantity signgam.
  11
  12         The coefficients for expansion around zero
  13         are #5243 from Hart & Cheney; for expansion
  14         around infinity they are #5404.
  15
  16         Calls log and sin.
  17 */
  18
  19 int     signgam;
  20 static double goobie    = 0.9189385332046727417803297;
  21 static double pi        = 3.1415926535897932384626434;
  22
  23 #define M 6
  24 #define N 8
  25 static double p1[] = {
  26         0.83333333333333101837e-1,
  27         -.277777777735865004e-2,
  28         0.793650576493454e-3,
  29         -.5951896861197e-3,
  30         0.83645878922e-3,
  31         -.1633436431e-2,
  32 };
  33 static double p2[] = {
  34         -.42353689509744089647e5,
  35         -.20886861789269887364e5,
  36         -.87627102978521489560e4,
  37         -.20085274013072791214e4,
  38         -.43933044406002567613e3,
  39         -.50108693752970953015e2,
  40         -.67449507245925289918e1,
  41         0.0,
  42 };
  43 static double q2[] = {
  44         -.42353689509744090010e5,
  45         -.29803853309256649932e4,
  46         0.99403074150827709015e4,
  47         -.15286072737795220248e4,
  48         -.49902852662143904834e3,
  49         0.18949823415702801641e3,
  50         -.23081551524580124562e2,
  51         0.10000000000000000000e1,
  52 };
  53
  54 static double
  55 asym(double arg)
  56 {
  57         double n, argsq;
  58         int i;
  59
  60         argsq = 1 / (arg*arg);
  61         n = 0;
  62         for(i=M-1; i>=0; i--)
  63                 n = n*argsq + p1[i];
  64         return (arg-.5)*log(arg) - arg + goobie + n/arg;
  65 }
  66
  67 static double
  68 pos(double arg)
  69 {
  70         double n, d, s;
  71         int i;
  72
  73         if(arg < 2)
  74                 return pos(arg+1)/arg;
  75         if(arg > 3)
  76                 return (arg-1)*pos(arg-1);
  77
  78         s = arg - 2;
  79         n = 0;
  80         d = 0;
  81         for(i=N-1; i>=0; i--){
  82                 n = n*s + p2[i];
  83                 d = d*s + q2[i];
  84         }
  85         return n/d;
  86 }
  87
  88 static double
  89 neg(double arg)
  90 {
  91         double temp;
  92
  93         arg = -arg;
  94         temp = sin(pi*arg);
  95         if(temp == 0) {
  96                 errno = EDOM;
  97                 return HUGE_VAL;
  98         }
  99         if(temp < 0)
 100                 temp = -temp;
 101         else
 102                 signgam = -1;
 103         return -log(arg*pos(arg)*temp/pi);
 104 }
 105
 106 double
 107 gamma(double arg)
 108 {
 109
 110         signgam = 1;
 111         if(arg <= 0)
 112                 return neg(arg);
 113         if(arg > 8)
 114                 return asym(arg);
 115         return log(pos(arg));
 116 }