1*22554Smiriam /* @(#)lgamma.c 4.2 06/06/85 */ 29931Ssam 39931Ssam /* 49931Ssam C program for floating point log gamma function 59931Ssam 69931Ssam gamma(x) computes the log of the absolute 79931Ssam value of the gamma function. 89931Ssam The sign of the gamma function is returned in the 99931Ssam external quantity signgam. 109931Ssam 119931Ssam The coefficients for expansion around zero 129931Ssam are #5243 from Hart & Cheney; for expansion 139931Ssam around infinity they are #5404. 149931Ssam 15*22554Smiriam Calls log, drem and sin. 169931Ssam */ 179931Ssam 189931Ssam #include <errno.h> 199931Ssam #include <math.h> 20*22554Smiriam #ifdef VAX 21*22554Smiriam static long NaN_[] = {0x8000, 0x0}; 22*22554Smiriam #define NaN (*(double *) NaN_) 23*22554Smiriam #endif 249931Ssam 25*22554Smiriam 269931Ssam int errno; 279931Ssam int signgam = 0; 28*22554Smiriam static double goobie = 0.9189385332046727417803297; /* log(2*pi)/2 */ 299931Ssam static double pi = 3.1415926535897932384626434; 309931Ssam 319931Ssam #define M 6 329931Ssam #define N 8 339931Ssam static double p1[] = { 349931Ssam 0.83333333333333101837e-1, 359931Ssam -.277777777735865004e-2, 369931Ssam 0.793650576493454e-3, 379931Ssam -.5951896861197e-3, 389931Ssam 0.83645878922e-3, 399931Ssam -.1633436431e-2, 409931Ssam }; 419931Ssam static double p2[] = { 429931Ssam -.42353689509744089647e5, 439931Ssam -.20886861789269887364e5, 449931Ssam -.87627102978521489560e4, 459931Ssam -.20085274013072791214e4, 469931Ssam -.43933044406002567613e3, 479931Ssam -.50108693752970953015e2, 489931Ssam -.67449507245925289918e1, 499931Ssam 0.0, 509931Ssam }; 519931Ssam static double q2[] = { 529931Ssam -.42353689509744090010e5, 539931Ssam -.29803853309256649932e4, 549931Ssam 0.99403074150827709015e4, 559931Ssam -.15286072737795220248e4, 569931Ssam -.49902852662143904834e3, 579931Ssam 0.18949823415702801641e3, 589931Ssam -.23081551524580124562e2, 599931Ssam 0.10000000000000000000e1, 609931Ssam }; 619931Ssam 629931Ssam double 639931Ssam gamma(arg) 649931Ssam double arg; 659931Ssam { 669931Ssam double log(), pos(), neg(), asym(); 679931Ssam 689931Ssam signgam = 1.; 699931Ssam if(arg <= 0.) return(neg(arg)); 709931Ssam if(arg > 8.) return(asym(arg)); 719931Ssam return(log(pos(arg))); 729931Ssam } 739931Ssam 749931Ssam static double 759931Ssam asym(arg) 769931Ssam double arg; 779931Ssam { 789931Ssam double log(); 799931Ssam double n, argsq; 809931Ssam int i; 819931Ssam 829931Ssam argsq = 1./(arg*arg); 839931Ssam for(n=0,i=M-1; i>=0; i--){ 849931Ssam n = n*argsq + p1[i]; 859931Ssam } 869931Ssam return((arg-.5)*log(arg) - arg + goobie + n/arg); 879931Ssam } 889931Ssam 899931Ssam static double 909931Ssam neg(arg) 919931Ssam double arg; 929931Ssam { 939931Ssam double temp; 94*22554Smiriam double log(), sin(), drem(), pos(); 959931Ssam 969931Ssam arg = -arg; 97*22554Smiriam /* 98*22554Smiriam * to see if arg were a true integer, the old code used the 99*22554Smiriam * mathematically correct observation: 100*22554Smiriam * sin(n*pi) = 0 <=> n is an integer. 101*22554Smiriam * but in finite precision arithmetic, sin(n*PI) will NEVER 102*22554Smiriam * be zero simply because n*PI is a rational number. hence 103*22554Smiriam * it failed to work with our newer, more accurate sin() 104*22554Smiriam * which uses true pi to do the argument reduction... 105*22554Smiriam * temp = sin(pi*arg); 106*22554Smiriam */ 107*22554Smiriam temp = drem(arg, 1.e0); 1089931Ssam if(temp == 0.) { 1099931Ssam errno = EDOM; 110*22554Smiriam #ifdef VAX 111*22554Smiriam return (NaN); 112*22554Smiriam #else 1139931Ssam return(HUGE); 114*22554Smiriam #endif 1159931Ssam } 116*22554Smiriam 117*22554Smiriam temp = drem(arg, 2.e0); 118*22554Smiriam if (temp < 0.) 119*22554Smiriam temp = -temp; 120*22554Smiriam else 121*22554Smiriam signgam = -1; 122*22554Smiriam return(-log(arg*pos(arg)*sin(pi*temp)/pi)); 1239931Ssam } 1249931Ssam 1259931Ssam static double 1269931Ssam pos(arg) 1279931Ssam double arg; 1289931Ssam { 1299931Ssam double n, d, s; 1309931Ssam register i; 1319931Ssam 1329931Ssam if(arg < 2.) return(pos(arg+1.)/arg); 1339931Ssam if(arg > 3.) return((arg-1.)*pos(arg-1.)); 1349931Ssam 1359931Ssam s = arg - 2.; 1369931Ssam for(n=0,d=0,i=N-1; i>=0; i--){ 1379931Ssam n = n*s + p2[i]; 1389931Ssam d = d*s + q2[i]; 1399931Ssam } 1409931Ssam return(n/d); 1419931Ssam } 142