1 #ifndef lint 2 static char sccsid[] = "@(#)lgamma.c 4.4 (Berkeley) 09/11/85"; 3 #endif not lint 4 5 /* 6 C program for floating point log Gamma function 7 8 lgamma(x) computes the log of the absolute 9 value of the Gamma function. 10 The sign of the Gamma function is returned in the 11 external quantity signgam. 12 13 The coefficients for expansion around zero 14 are #5243 from Hart & Cheney; for expansion 15 around infinity they are #5404. 16 17 Calls log, floor and sin. 18 */ 19 20 #include <math.h> 21 #ifdef VAX 22 #include <errno.h> 23 #endif 24 int signgam = 0; 25 static double goobie = 0.9189385332046727417803297; /* log(2*pi)/2 */ 26 static double pi = 3.1415926535897932384626434; 27 28 #define M 6 29 #define N 8 30 static double p1[] = { 31 0.83333333333333101837e-1, 32 -.277777777735865004e-2, 33 0.793650576493454e-3, 34 -.5951896861197e-3, 35 0.83645878922e-3, 36 -.1633436431e-2, 37 }; 38 static double p2[] = { 39 -.42353689509744089647e5, 40 -.20886861789269887364e5, 41 -.87627102978521489560e4, 42 -.20085274013072791214e4, 43 -.43933044406002567613e3, 44 -.50108693752970953015e2, 45 -.67449507245925289918e1, 46 0.0, 47 }; 48 static double q2[] = { 49 -.42353689509744090010e5, 50 -.29803853309256649932e4, 51 0.99403074150827709015e4, 52 -.15286072737795220248e4, 53 -.49902852662143904834e3, 54 0.18949823415702801641e3, 55 -.23081551524580124562e2, 56 0.10000000000000000000e1, 57 }; 58 59 double 60 lgamma(arg) 61 double arg; 62 { 63 double log(), pos(), neg(), asym(); 64 65 signgam = 1.; 66 if(arg <= 0.) return(neg(arg)); 67 if(arg > 8.) return(asym(arg)); 68 return(log(pos(arg))); 69 } 70 71 static double 72 asym(arg) 73 double arg; 74 { 75 double log(); 76 double n, argsq; 77 int i; 78 79 argsq = 1./(arg*arg); 80 for(n=0,i=M-1; i>=0; i--){ 81 n = n*argsq + p1[i]; 82 } 83 return((arg-.5)*log(arg) - arg + goobie + n/arg); 84 } 85 86 static double 87 neg(arg) 88 double arg; 89 { 90 double t; 91 double log(), sin(), floor(), pos(); 92 93 arg = -arg; 94 /* 95 * to see if arg were a true integer, the old code used the 96 * mathematically correct observation: 97 * sin(n*pi) = 0 <=> n is an integer. 98 * but in finite precision arithmetic, sin(n*PI) will NEVER 99 * be zero simply because n*PI is a rational number. hence 100 * it failed to work with our newer, more accurate sin() 101 * which uses true pi to do the argument reduction... 102 * temp = sin(pi*arg); 103 */ 104 t = floor(arg); 105 if (arg - t > 0.5e0) 106 t += 1.e0; /* t := integer nearest arg */ 107 #ifdef VAX 108 if (arg == t) { 109 extern double infnan(); 110 return(infnan(ERANGE)); /* +INF */ 111 } 112 #endif 113 signgam = (int) (t - 2*floor(t/2)); /* signgam = 1 if t was odd, */ 114 /* 0 if t was even */ 115 signgam = signgam - 1 + signgam; /* signgam = 1 if t was odd, */ 116 /* -1 if t was even */ 117 t = arg - t; /* -0.5 <= t <= 0.5 */ 118 if (t < 0.e0) { 119 t = -t; 120 signgam = -signgam; 121 } 122 return(-log(arg*pos(arg)*sin(pi*t)/pi)); 123 } 124 125 static double 126 pos(arg) 127 double arg; 128 { 129 double n, d, s; 130 register i; 131 132 if(arg < 2.) return(pos(arg+1.)/arg); 133 if(arg > 3.) return((arg-1.)*pos(arg-1.)); 134 135 s = arg - 2.; 136 for(n=0,d=0,i=N-1; i>=0; i--){ 137 n = n*s + p2[i]; 138 d = d*s + q2[i]; 139 } 140 return(n/d); 141 } 142