1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
/* rgamma.c
*
* Reciprocal Gamma function
*
*
*
* SYNOPSIS:
*
* double x, y, rgamma();
*
* y = rgamma( x );
*
*
*
* DESCRIPTION:
*
* Returns one divided by the Gamma function of the argument.
*
* The function is approximated by a Chebyshev expansion in
* the interval [0,1]. Range reduction is by recurrence
* for arguments between -34.034 and +34.84425627277176174.
* 0 is returned for positive arguments outside this
* range. For arguments less than -34.034 the cosecant
* reflection formula is applied; lograrithms are employed
* to avoid unnecessary overflow.
*
* The reciprocal Gamma function has no singularities,
* but overflow and underflow may occur for large arguments.
* These conditions return either NPY_INFINITY or 0 with
* appropriate sign.
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -30,+30 30000 1.1e-15 2.0e-16
* For arguments less than -34.034 the peak error is on the
* order of 5e-15 (DEC), excepting overflow or underflow.
*/
/*
* Cephes Math Library Release 2.0: April, 1987
* Copyright 1985, 1987 by Stephen L. Moshier
* Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include "mconf.h"
/* Chebyshev coefficients for reciprocal Gamma function
* in interval 0 to 1. Function is 1/(x Gamma(x)) - 1
*/
static double R[] = {
3.13173458231230000000E-17,
-6.70718606477908000000E-16,
2.20039078172259550000E-15,
2.47691630348254132600E-13,
-6.60074100411295197440E-12,
5.13850186324226978840E-11,
1.08965386454418662084E-9,
-3.33964630686836942556E-8,
2.68975996440595483619E-7,
2.96001177518801696639E-6,
-8.04814124978471142852E-5,
4.16609138709688864714E-4,
5.06579864028608725080E-3,
-6.41925436109158228810E-2,
-4.98558728684003594785E-3,
1.27546015610523951063E-1
};
static char name[] = "rgamma";
extern double MAXLOG;
double rgamma(x)
double x;
{
double w, y, z;
int sign;
if (x > 34.84425627277176174) {
return exp(-lgam(x));
}
if (x < -34.034) {
w = -x;
z = sin(NPY_PI * w);
if (z == 0.0)
return (0.0);
if (z < 0.0) {
sign = 1;
z = -z;
}
else
sign = -1;
y = log(w * z) - log(NPY_PI) + lgam(w);
if (y < -MAXLOG) {
mtherr(name, UNDERFLOW);
return (sign * 0.0);
}
if (y > MAXLOG) {
mtherr(name, OVERFLOW);
return (sign * NPY_INFINITY);
}
return (sign * exp(y));
}
z = 1.0;
w = x;
while (w > 1.0) { /* Downward recurrence */
w -= 1.0;
z *= w;
}
while (w < 0.0) { /* Upward recurrence */
z /= w;
w += 1.0;
}
if (w == 0.0) /* Nonpositive integer */
return (0.0);
if (w == 1.0) /* Other integer */
return (1.0 / z);
y = w * (1.0 + chbevl(4.0 * w - 2.0, R, 16)) / z;
return (y);
}
|
