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 } 117