← Back to branch summary

~ubuntu-branches/ubuntu/precise/dbacl/precise

« back to all changes in this revision

Viewing changes to src/igam.c

  • Committer: Bazaar Package Importer
  • Author(s): Clint Adams
  • Date: 2005-05-07 12:59:53 UTC
  • Revision ID: james.westby@ubuntu.com-20050507125953-xzy2bwkb2qamglwm
Tags: upstream-1.9
Import upstream version 1.9

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/igam.c
 
1
/*                                                      igam.c
 
2
 *
 
3
 *      Incomplete gamma integral
 
4
 *
 
5
 *
 
6
 *
 
7
 * SYNOPSIS:
 
8
 *
 
9
 * double a, x, y, igam();
 
10
 *
 
11
 * y = igam( a, x );
 
12
 *
 
13
 * DESCRIPTION:
 
14
 *
 
15
 * The function is defined by
 
16
 *
 
17
 *                           x
 
18
 *                            -
 
19
 *                   1       | |  -t  a-1
 
20
 *  igam(a,x)  =   -----     |   e   t   dt.
 
21
 *                  -      | |
 
22
 *                 | (a)    -
 
23
 *                           0
 
24
 *
 
25
 *
 
26
 * In this implementation both arguments must be positive.
 
27
 * The integral is evaluated by either a power series or
 
28
 * continued fraction expansion, depending on the relative
 
29
 * values of a and x.
 
30
 *
 
31
 * ACCURACY:
 
32
 *
 
33
 *                      Relative error:
 
34
 * arithmetic   domain     # trials      peak         rms
 
35
 *    IEEE      0,30       200000       3.6e-14     2.9e-15
 
36
 *    IEEE      0,100      300000       9.9e-14     1.5e-14
 
37
 */
 
38
/*                                                     igamc()
 
39
 *
 
40
 *      Complemented incomplete gamma integral
 
41
 *
 
42
 *
 
43
 *
 
44
 * SYNOPSIS:
 
45
 *
 
46
 * double a, x, y, igamc();
 
47
 *
 
48
 * y = igamc( a, x );
 
49
 *
 
50
 * DESCRIPTION:
 
51
 *
 
52
 * The function is defined by
 
53
 *
 
54
 *
 
55
 *  igamc(a,x)   =   1 - igam(a,x)
 
56
 *
 
57
 *                            inf.
 
58
 *                              -
 
59
 *                     1       | |  -t  a-1
 
60
 *               =   -----     |   e   t   dt.
 
61
 *                    -      | |
 
62
 *                   | (a)    -
 
63
 *                             x
 
64
 *
 
65
 *
 
66
 * In this implementation both arguments must be positive.
 
67
 * The integral is evaluated by either a power series or
 
68
 * continued fraction expansion, depending on the relative
 
69
 * values of a and x.
 
70
 *
 
71
 * ACCURACY:
 
72
 *
 
73
 * Tested at random a, x.
 
74
 *                a         x                      Relative error:
 
75
 * arithmetic   domain   domain     # trials      peak         rms
 
76
 *    IEEE     0.5,100   0,100      200000       1.9e-14     1.7e-15
 
77
 *    IEEE     0.01,0.5  0,100      200000       1.4e-13     1.6e-15
 
78
 */
 
79
 
 
80
/*
 
81
Cephes Math Library Release 2.8:  June, 2000
 
82
Copyright 1985, 1987, 2000 by Stephen L. Moshier
 
83
*/
 
84
 
 
85
#include "mconf.h"
 
86
#ifdef ANSIPROT
 
87
extern double lgam ( double );
 
88
extern double exp ( double );
 
89
extern double log ( double );
 
90
extern double fabs ( double );
 
91
extern double igam ( double, double );
 
92
extern double igamc ( double, double );
 
93
#else
 
94
double lgam(), exp(), log(), fabs(), igam(), igamc();
 
95
#endif
 
96
 
 
97
extern double MACHEP, MAXLOG;
 
98
static double big = 4.503599627370496e15;
 
99
static double biginv =  2.22044604925031308085e-16;
 
100
 
 
101
double igamc( a, x )
 
102
double a, x;
 
103
{
 
104
double ans, ax, c, yc, r, t, y, z;
 
105
double pk, pkm1, pkm2, qk, qkm1, qkm2;
 
106
 
 
107
if( (x <= 0) || ( a <= 0) )
 
108
        return( 1.0 );
 
109
 
 
110
if( (x < 1.0) || (x < a) )
 
111
        return( 1.0 - igam(a,x) );
 
112
 
 
113
ax = a * log(x) - x - lgam(a);
 
114
if( ax < -MAXLOG )
 
115
        {
 
116
        mtherr( "igamc", UNDERFLOW );
 
117
        return( 0.0 );
 
118
        }
 
119
ax = exp(ax);
 
120
 
 
121
/* continued fraction */
 
122
y = 1.0 - a;
 
123
z = x + y + 1.0;
 
124
c = 0.0;
 
125
pkm2 = 1.0;
 
126
qkm2 = x;
 
127
pkm1 = x + 1.0;
 
128
qkm1 = z * x;
 
129
ans = pkm1/qkm1;
 
130
 
 
131
do
 
132
        {
 
133
        c += 1.0;
 
134
        y += 1.0;
 
135
        z += 2.0;
 
136
        yc = y * c;
 
137
        pk = pkm1 * z  -  pkm2 * yc;
 
138
        qk = qkm1 * z  -  qkm2 * yc;
 
139
        if( qk != 0 )
 
140
                {
 
141
                r = pk/qk;
 
142
                t = fabs( (ans - r)/r );
 
143
                ans = r;
 
144
                }
 
145
        else
 
146
                t = 1.0;
 
147
        pkm2 = pkm1;
 
148
        pkm1 = pk;
 
149
        qkm2 = qkm1;
 
150
        qkm1 = qk;
 
151
        if( fabs(pk) > big )
 
152
                {
 
153
                pkm2 *= biginv;
 
154
                pkm1 *= biginv;
 
155
                qkm2 *= biginv;
 
156
                qkm1 *= biginv;
 
157
                }
 
158
        }
 
159
while( t > MACHEP );
 
160
 
 
161
return( ans * ax );
 
162
}
 
163
 
 
164
 
 
165
 
 
166
/* left tail of incomplete gamma function:
 
167
 *
 
168
 *          inf.      k
 
169
 *   a  -x   -       x
 
170
 *  x  e     >   ----------
 
171
 *           -     -
 
172
 *          k=0   | (a+k+1)
 
173
 *
 
174
 */
 
175
 
 
176
double igam( a, x )
 
177
double a, x;
 
178
{
 
179
double ans, ax, c, r;
 
180
 
 
181
if( (x <= 0) || ( a <= 0) )
 
182
        return( 0.0 );
 
183
 
 
184
if( (x > 1.0) && (x > a ) )
 
185
        return( 1.0 - igamc(a,x) );
 
186
 
 
187
/* Compute  x**a * exp(-x) / gamma(a)  */
 
188
ax = a * log(x) - x - lgam(a);
 
189
if( ax < -MAXLOG )
 
190
        {
 
191
        mtherr( "igam", UNDERFLOW );
 
192
        return( 0.0 );
 
193
        }
 
194
ax = exp(ax);
 
195
 
 
196
/* power series */
 
197
r = a;
 
198
c = 1.0;
 
199
ans = 1.0;
 
200
 
 
201
do
 
202
        {
 
203
        r += 1.0;
 
204
        c *= x/r;
 
205
        ans += c;
 
206
        }
 
207
while( c/ans > MACHEP );
 
208
 
 
209
return( ans * ax/a );
 
210
}