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