1*24691Smckusick #ifndef lint
2*24691Smckusick static char sccsid[] = "@(#)lgamma.c 4.4 (Berkeley) 09/11/85";
3*24691Smckusick #endif not lint
49931Ssam
59931Ssam /*
6*24691Smckusick C program for floating point log Gamma function
79931Ssam
8*24691Smckusick lgamma(x) computes the log of the absolute
9*24691Smckusick value of the Gamma function.
10*24691Smckusick The sign of the Gamma function is returned in the
119931Ssam external quantity signgam.
129931Ssam
139931Ssam The coefficients for expansion around zero
149931Ssam are #5243 from Hart & Cheney; for expansion
159931Ssam around infinity they are #5404.
169931Ssam
1723563Smiriam Calls log, floor and sin.
189931Ssam */
199931Ssam
209931Ssam #include <math.h>
2122554Smiriam #ifdef VAX
22*24691Smckusick #include <errno.h>
2322554Smiriam #endif
249931Ssam int signgam = 0;
2522554Smiriam static double goobie = 0.9189385332046727417803297; /* log(2*pi)/2 */
269931Ssam static double pi = 3.1415926535897932384626434;
279931Ssam
289931Ssam #define M 6
299931Ssam #define N 8
309931Ssam static double p1[] = {
319931Ssam 0.83333333333333101837e-1,
329931Ssam -.277777777735865004e-2,
339931Ssam 0.793650576493454e-3,
349931Ssam -.5951896861197e-3,
359931Ssam 0.83645878922e-3,
369931Ssam -.1633436431e-2,
379931Ssam };
389931Ssam static double p2[] = {
399931Ssam -.42353689509744089647e5,
409931Ssam -.20886861789269887364e5,
419931Ssam -.87627102978521489560e4,
429931Ssam -.20085274013072791214e4,
439931Ssam -.43933044406002567613e3,
449931Ssam -.50108693752970953015e2,
459931Ssam -.67449507245925289918e1,
469931Ssam 0.0,
479931Ssam };
489931Ssam static double q2[] = {
499931Ssam -.42353689509744090010e5,
509931Ssam -.29803853309256649932e4,
519931Ssam 0.99403074150827709015e4,
529931Ssam -.15286072737795220248e4,
539931Ssam -.49902852662143904834e3,
549931Ssam 0.18949823415702801641e3,
559931Ssam -.23081551524580124562e2,
569931Ssam 0.10000000000000000000e1,
579931Ssam };
589931Ssam
599931Ssam double
lgamma(arg)60*24691Smckusick lgamma(arg)
619931Ssam double arg;
629931Ssam {
639931Ssam double log(), pos(), neg(), asym();
649931Ssam
659931Ssam signgam = 1.;
669931Ssam if(arg <= 0.) return(neg(arg));
679931Ssam if(arg > 8.) return(asym(arg));
689931Ssam return(log(pos(arg)));
699931Ssam }
709931Ssam
719931Ssam static double
asym(arg)729931Ssam asym(arg)
739931Ssam double arg;
749931Ssam {
759931Ssam double log();
769931Ssam double n, argsq;
779931Ssam int i;
789931Ssam
799931Ssam argsq = 1./(arg*arg);
809931Ssam for(n=0,i=M-1; i>=0; i--){
819931Ssam n = n*argsq + p1[i];
829931Ssam }
839931Ssam return((arg-.5)*log(arg) - arg + goobie + n/arg);
849931Ssam }
859931Ssam
869931Ssam static double
neg(arg)879931Ssam neg(arg)
889931Ssam double arg;
899931Ssam {
9023563Smiriam double t;
9123563Smiriam double log(), sin(), floor(), pos();
929931Ssam
939931Ssam arg = -arg;
9422554Smiriam /*
9522554Smiriam * to see if arg were a true integer, the old code used the
9622554Smiriam * mathematically correct observation:
9722554Smiriam * sin(n*pi) = 0 <=> n is an integer.
9822554Smiriam * but in finite precision arithmetic, sin(n*PI) will NEVER
9922554Smiriam * be zero simply because n*PI is a rational number. hence
10022554Smiriam * it failed to work with our newer, more accurate sin()
10122554Smiriam * which uses true pi to do the argument reduction...
10222554Smiriam * temp = sin(pi*arg);
10322554Smiriam */
10423563Smiriam t = floor(arg);
10523563Smiriam if (arg - t > 0.5e0)
106*24691Smckusick t += 1.e0; /* t := integer nearest arg */
10722554Smiriam #ifdef VAX
108*24691Smckusick if (arg == t) {
109*24691Smckusick extern double infnan();
110*24691Smckusick return(infnan(ERANGE)); /* +INF */
1119931Ssam }
11223563Smiriam #endif
11323563Smiriam signgam = (int) (t - 2*floor(t/2)); /* signgam = 1 if t was odd, */
11423563Smiriam /* 0 if t was even */
11523563Smiriam signgam = signgam - 1 + signgam; /* signgam = 1 if t was odd, */
11623563Smiriam /* -1 if t was even */
11723563Smiriam t = arg - t; /* -0.5 <= t <= 0.5 */
11823563Smiriam if (t < 0.e0) {
11923563Smiriam t = -t;
12023563Smiriam signgam = -signgam;
12123563Smiriam }
12223563Smiriam return(-log(arg*pos(arg)*sin(pi*t)/pi));
1239931Ssam }
1249931Ssam
1259931Ssam static double
pos(arg)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