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