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