← Back to branch summary

~ubuntu-branches/ubuntu/karmic/grace/karmic

« back to all changes in this revision

Viewing changes to cephes/igam.c

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-19 14:19:58 UTC
  • Revision ID: james.westby@ubuntu.com-20020319141958-5gxna6vo1ek3zjml
Tags: upstream-5.1.7
Import upstream version 5.1.7

Show diffs side-by-side

added added

removed removed

Lines of Context:
cephes/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.0:  April, 1987
 
82
Copyright 1985, 1987 by Stephen L. Moshier
 
83
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 
84
*/
 
85
 
 
86
#include "mconf.h"
 
87
#include "cephes.h"
 
88
 
 
89
extern double MACHEP, MAXLOG;
 
90
static double big = 4.503599627370496e15;
 
91
static double biginv =  2.22044604925031308085e-16;
 
92
 
 
93
double igamc( a, x )
 
94
double a, x;
 
95
{
 
96
double ans, ax, c, yc, r, t, y, z;
 
97
double pk, pkm1, pkm2, qk, qkm1, qkm2;
 
98
 
 
99
if( (x <= 0) || ( a <= 0) )
 
100
        return( 1.0 );
 
101
 
 
102
if( (x < 1.0) || (x < a) )
 
103
        return( 1.0 - igam(a,x) );
 
104
 
 
105
ax = a * log(x) - x - lgam(a);
 
106
if( ax < -MAXLOG )
 
107
        {
 
108
        mtherr( "igamc", UNDERFLOW );
 
109
        return( 0.0 );
 
110
        }
 
111
ax = exp(ax);
 
112
 
 
113
/* continued fraction */
 
114
y = 1.0 - a;
 
115
z = x + y + 1.0;
 
116
c = 0.0;
 
117
pkm2 = 1.0;
 
118
qkm2 = x;
 
119
pkm1 = x + 1.0;
 
120
qkm1 = z * x;
 
121
ans = pkm1/qkm1;
 
122
 
 
123
do
 
124
        {
 
125
        c += 1.0;
 
126
        y += 1.0;
 
127
        z += 2.0;
 
128
        yc = y * c;
 
129
        pk = pkm1 * z  -  pkm2 * yc;
 
130
        qk = qkm1 * z  -  qkm2 * yc;
 
131
        if( qk != 0 )
 
132
                {
 
133
                r = pk/qk;
 
134
                t = fabs( (ans - r)/r );
 
135
                ans = r;
 
136
                }
 
137
        else
 
138
                t = 1.0;
 
139
        pkm2 = pkm1;
 
140
        pkm1 = pk;
 
141
        qkm2 = qkm1;
 
142
        qkm1 = qk;
 
143
        if( fabs(pk) > big )
 
144
                {
 
145
                pkm2 *= biginv;
 
146
                pkm1 *= biginv;
 
147
                qkm2 *= biginv;
 
148
                qkm1 *= biginv;
 
149
                }
 
150
        }
 
151
while( t > MACHEP );
 
152
 
 
153
return( ans * ax );
 
154
}
 
155
 
 
156
 
 
157
 
 
158
/* left tail of incomplete gamma function:
 
159
 *
 
160
 *          inf.      k
 
161
 *   a  -x   -       x
 
162
 *  x  e     >   ----------
 
163
 *           -     -
 
164
 *          k=0   | (a+k+1)
 
165
 *
 
166
 */
 
167
 
 
168
double igam( a, x )
 
169
double a, x;
 
170
{
 
171
double ans, ax, c, r;
 
172
 
 
173
if( (x <= 0) || ( a <= 0) )
 
174
        return( 0.0 );
 
175
 
 
176
if( (x > 1.0) && (x > a ) )
 
177
        return( 1.0 - igamc(a,x) );
 
178
 
 
179
/* Compute  x**a * exp(-x) / gamma(a)  */
 
180
ax = a * log(x) - x - lgam(a);
 
181
if( ax < -MAXLOG )
 
182
        {
 
183
        mtherr( "igam", UNDERFLOW );
 
184
        return( 0.0 );
 
185
        }
 
186
ax = exp(ax);
 
187
 
 
188
/* power series */
 
189
r = a;
 
190
c = 1.0;
 
191
ans = 1.0;
 
192
 
 
193
do
 
194
        {
 
195
        r += 1.0;
 
196
        c *= x/r;
 
197
        ans += c;
 
198
        }
 
199
while( c/ans > MACHEP );
 
200
 
 
201
return( ans * ax/a );
 
202
}