1*23563Smiriam /* @(#)lgamma.c 4.3 06/19/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*23563Smiriam Calls log, floor and sin. 169931Ssam */ 179931Ssam 189931Ssam #include <errno.h> 199931Ssam #include <math.h> 2022554Smiriam #ifdef VAX 2122554Smiriam static long NaN_[] = {0x8000, 0x0}; 2222554Smiriam #define NaN (*(double *) NaN_) 2322554Smiriam #endif 249931Ssam 2522554Smiriam 269931Ssam int errno; 279931Ssam int signgam = 0; 2822554Smiriam 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 { 93*23563Smiriam double t; 94*23563Smiriam double log(), sin(), floor(), pos(); 959931Ssam 969931Ssam arg = -arg; 9722554Smiriam /* 9822554Smiriam * to see if arg were a true integer, the old code used the 9922554Smiriam * mathematically correct observation: 10022554Smiriam * sin(n*pi) = 0 <=> n is an integer. 10122554Smiriam * but in finite precision arithmetic, sin(n*PI) will NEVER 10222554Smiriam * be zero simply because n*PI is a rational number. hence 10322554Smiriam * it failed to work with our newer, more accurate sin() 10422554Smiriam * which uses true pi to do the argument reduction... 10522554Smiriam * temp = sin(pi*arg); 10622554Smiriam */ 107*23563Smiriam t = floor(arg); 108*23563Smiriam if (arg - t > 0.5e0) 109*23563Smiriam t += 1.e0; /* t := integer nearest arg */ 110*23563Smiriam #if !(IEEE|NATIONAL) 111*23563Smiriam if(arg == t) { 1129931Ssam errno = EDOM; 11322554Smiriam #ifdef VAX 11422554Smiriam return (NaN); 11522554Smiriam #else 1169931Ssam return(HUGE); 11722554Smiriam #endif 1189931Ssam } 119*23563Smiriam #endif 12022554Smiriam 121*23563Smiriam signgam = (int) (t - 2*floor(t/2)); /* signgam = 1 if t was odd, */ 122*23563Smiriam /* 0 if t was even */ 123*23563Smiriam signgam = signgam - 1 + signgam; /* signgam = 1 if t was odd, */ 124*23563Smiriam /* -1 if t was even */ 125*23563Smiriam t = arg - t; /* -0.5 <= t <= 0.5 */ 126*23563Smiriam if (t < 0.e0) { 127*23563Smiriam t = -t; 128*23563Smiriam signgam = -signgam; 129*23563Smiriam } 130*23563Smiriam return(-log(arg*pos(arg)*sin(pi*t)/pi)); 1319931Ssam } 1329931Ssam 1339931Ssam static double 1349931Ssam pos(arg) 1359931Ssam double arg; 1369931Ssam { 1379931Ssam double n, d, s; 1389931Ssam register i; 1399931Ssam 1409931Ssam if(arg < 2.) return(pos(arg+1.)/arg); 1419931Ssam if(arg > 3.) return((arg-1.)*pos(arg-1.)); 1429931Ssam 1439931Ssam s = arg - 2.; 1449931Ssam for(n=0,d=0,i=N-1; i>=0; i--){ 1459931Ssam n = n*s + p2[i]; 1469931Ssam d = d*s + q2[i]; 1479931Ssam } 1489931Ssam return(n/d); 1499931Ssam } 150