Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

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