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