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
128
129
130
|
/* rgammaf.c
*
* Reciprocal gamma function
*
*
*
* SYNOPSIS:
*
* float x, y, rgammaf();
*
* y = rgammaf( 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.
* 1/MAXNUMF is returned for positive arguments outside this
* range.
*
* The reciprocal gamma function has no singularities,
* but overflow and underflow may occur for large arguments.
* These conditions return either MAXNUMF or 1/MAXNUMF with
* appropriate sign.
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -34,+34 100000 8.9e-7 1.1e-7
*/
/*
Cephes Math Library Release 2.2: June, 1992
Copyright 1985, 1987, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include <math.h>
/* Chebyshev coefficients for reciprocal gamma function
* in interval 0 to 1. Function is 1/(x gamma(x)) - 1
*/
static float R[] = {
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[] = "rgammaf";
extern float PIF, MAXLOGF, MAXNUMF;
float chbevlf(float, float *, int);
float expf(float), logf(float), sinf(float), lgamf(float);
float rgammaf(float xx)
{
float x, w, y, z;
int sign;
x = xx;
if( x > 34.84425627277176174)
{
mtherr( name, UNDERFLOW );
return(1.0/MAXNUMF);
}
if( x < -34.034 )
{
w = -x;
z = sinf( PIF*w );
if( z == 0.0 )
return(0.0);
if( z < 0.0 )
{
sign = 1;
z = -z;
}
else
sign = -1;
y = logf( w * z / PIF ) + lgamf(w);
if( y < -MAXLOGF )
{
mtherr( name, UNDERFLOW );
return( sign * 1.0 / MAXNUMF );
}
if( y > MAXLOGF )
{
mtherr( name, OVERFLOW );
return( sign * MAXNUMF );
}
return( sign * expf(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 + chbevlf( 4.0*w-2.0, R, 10 ) ) / z;
return(y);
}
|
