The Museum of HP Calculators

HP Articles Forum

[Return to the Index ]
[ Previous | Next ]


Fast Incomplete Gamma/Exponential Integral/Error Functions With Lentz's Method

Posted by Steven Thomas Smith on 23 Sept 2009, 1:28 p.m.

The incomplete gamma function gamma(a,x) = Int0x e-tta-1dt, error function erf(x) = (2/sqrt pi)Int0x e-t2dt and exponential integrals En(x) = xn-1Intxinf e-tt-ndt are all closely related. The exponential integral E1(x) is also important because it provides the exponential integral Ei(x) = Int-infx e-tt-1dt = -Re(E1(-x)) (ignoring the imaginary -i*pi from integrating halfway around the pole at the origin) and the logarithmic integral li(x) = Ei(log x). The error and incomplete gamma functions P(a,x) = gamma(a,x)/Gamma(a) [lower tail] and Q(a,x) = 1 - P(a,x) [upper tail, "regularized" by dividing by Gamma(t)], are commonly encountered in the cumulative probability distribution of the normal (Gaussian) density and the chi-squared density functions, respectively. En(x) is essentially Gamma(-n+1,x) with the singularities at the nonpositive integers removed.

The continued fraction expansions of all these functions can be derived by integrating by parts the integral f(a,b) = (1/Gamma(a))Int0inf e-tta-1(1+xt)-bdt (see Wall, Analytic Theory of Continued Fractions, p. 354ff), which results in the following similar continued fractions:

                                      -x a  1/x  (1-a)/x 1/x  (2-a)/x 2/x
Gamma(a,x) = Gamma(a)-gamma(a,x) =   e  x   ---  ------- ---- ------- --- ...
                                            1 +   1 +    1 +   1 +    1 +

                            -x^2   1  1/2  1  3/2  2
erfc(x) = 1 - erf(x) =     e      --- --- --- --- --- ...
                       /sqrt(pi)  x + x + x + x + x +

            -x 1/x n/x 1/x (n+1)/x 2/x
E_n(x) =   e   --- --- --- ------- --- ...
               1 + 1 + 1 +  1 +    1 +

The simple, similar form of these continued fractions provides a fast, efficient algorithm for computing these functions for x larger than 2 or so. This is especially important for the exponential integrals, whose series representation alternates for positive x and therefore causes cancellation errors for moderately sized x.

Though backward recursion through the continued fraction is the most efficient, it is not possible to pre-compute the number of terms required to achieve the desired precision with two-parameter functions such as gamma(a,x) and En(x). Therefore, forward recursion is used, the best option being Lentz's Method (see Numerical Recipes, 2d ed., p. 207ff). Briefly, it is easy to show that the numerator An and denominator Bn of the continued fraction

     a1  a2
b0 + --- --- ...
     b1+ b2+

satisfy the same second-order linear recurrence relation:

An = bnAn-1 + anAn-2; A0 = b0, A-1 = 1
Bn = bnBn-1 + anBn-2; B0 = 1, B-1 = 0

Rather than computing An, Bn, and the nth order approximate An/Bn directly, the ratios Cn = An/An-1, Dn = Bn/Bn-1 are used, which satisfy the simpler first-order nonlinear recursion:

Cn = bn + an/Cn-1; C0 = HUGE_VAL
Dn = bn + an/Dn-1; D0 = 1

The nth approximate is then given by repeated multiplication (An/Bn) = (Cn/Dn)(An-1/Bn-1). Convergence is achieved when Cn/Dn = 1 to machine precision, or Cn = Dn to machine precision. The value for C0 = HUGE_VAL need simply be large enough so that C1 = b1 to machine precision, which for these functions at moderately sized xand larger means that HUGE_VAL = E99 on the HP41 is sufficient. In general, the modified Lentz's algorithm requires checking that neither Cn nor Dn becomes too small and adjusting them to control for underflow if they do, but these continued fractions are sufficiently well-behaved and this robustness checking is not required here.

The following HP41 algorithms compute all these functions using either a series (small-ish x), a backward recursion continued fraction if the number of necessary terms is known, or a forward recursion continued fraction using Lentz's Method, all in status registers XYZT/LMNO. Some of the series code is borrowed or modified from other sources, such as Jean-Marc Baillard's excellent collection of functions. It requires a GAMMA function (such as this one), and a digamma function for the higher order exponential integrals (such as Jean-Marc Baillard's PSI). Actual numerical software would use Chebyshev and/or rational approximations over many ranges of arguments -- but these practical approximation methods are much less interesting to understand and code!

The incomplete gamma function GAMMINC computes both upper and lower tails without underflow errors for small or large arguments, and for both ordinary and regularized functions.

a ENTER x XEQ GAMMINC produces the stacked results

T: Gamma(a,x)
Z: gamma(a,x)
Y: Q(a,x)
X: P(a,x)

Due to the efficient continued fraction algorithm performed in the status registers, GAMMINC runs about twice as fast as IGF.

The error function computes both upper and lower tails without underflow errors for small or large arguments. x XEQ ERF produces the stacked results

Z: erfcx(x) = exp(x^2)*erfc(x) for larger x
Y: erfc(x)
X: erf(x)

The exponential integrals E1(x) and En(x) are handled separately so that the special case Ei(x) = -Re(E1(-x)) with negative arguments can be addressed. The continued fraction for E1(x) is not convergent for negative x; therefore, an asymptotic expansion is used for large negative x, and both Ei and li may then be computed with the simple subroutines:

01 LBL "EI"
02 CHS
03 XEQ "E1"
04 CHS
05 END

01 LBL "LI"
02 LN
03 XEQ "EI"
04 END

The HP41 code for Gamma(a,x), gamma(a,x), Q(a,x), and P(a,x) [a ENTER x XEQ GAMMINC]:

01 LBL "GAMMINC"
02 STO 00
03 X<>Y
04 STO 01
05 XEQ "GAMMA"
06 SF 05
07 CF 06
08 X<> 01
09 X<0?
10 SF 06
11 STO N
12 RCL 00
13 STO M
14 X<=Y?
15 CF 05
16 5
17 FS?C 06
18 GTO 00
19 X>Y?
20 CF 05
21 LBL 00
22 1/X
23 X>Y?
24 CF 05
25 LBL 00
26 RDN
27 X!=0?
28 GTO 00
29 RCL 01
30 .
31 E
32 .
33 RTN
34 LBL 00
35 RCL N
36 RCL M
37 LN
38 *
39 RCL M
40 -
41 E^X
42 FS? 05
43 GTO 02
44 RCL N
45 /
46 STO O
47 LBL 01
48 RCL M
49 RCL N
50 E
51 +
52 STO N
53 /
54 RCL O
55 *
56 STO O
57 +
58 X!=Y?
59 GTO 01
60 GTO 07
61 LBL 02
62 E
63 X<> N
64 STO L
65 X<>Y
66 E99
67 X<> M
68 RCL N
69 RCL Z
70 RCL Z
71 ST/ Y
72 RDN
73 STO O
74 CLX
75 RCL M
76 RCL N
77 LBL 03
78 X=Y?
79 GTO 06
80 R^
81 R^
82 RCL X
83 LASTX
84 -
85 R^
86 ST/ Y
87 CLX
88 RCL Y
89 SF 06
90 LBL 04
91 X<> M
92 ST/ M
93 X<>Y
94 X<> N
95 ST/ N
96 CLX
97 E
98 ST+ M
99 ST+ N
100 R^
101 R^
102 RCL M
103 ST* O
104 RCL N
105 ST/ O
106 X=Y?
107 GTO 06
108 FC?C 06
109 GTO 05
110 R^
111 R^
112 RCL X
113 RCL Z
114 ST/ Y
115 CLX
116 RCL Y
117 GTO 04
118 LBL 05
119 ISG Z
120 ""     ;F0
121 GTO 03
122 LBL 06
123 R^
124 STO M
125 LASTX
126 STO N
127 RCL O
128 LBL 07
129 RCL 01
130 X<>Y
131 -
132 LASTX
133 FS?C 05
134 X<>Y
135 STO O
136 RCL Y
137 RCL 01
138 ST/ Y
139 ST/ O
140 CLX
141 RCL O
142 END

The HP41 code for erf(x) [x XEQ ERF]:

01 LBL "ERF"
02 CF 05
03 CF 06
04 X<=0?
05 SF 06
06 ABS
07 STO M
08 6
09 LN
10 X<=Y?
11 GTO 03
12 RCL M
13 X^2
14 STO N
15 E
16 STO O
17 .
18 LASTX
19 .
20 LBL 01
21 X=Y?
22 GTO 02
23 ISG Z
24 ""
25 CLX
26 RCL M
27 RCL N
28 *
29 RCL Z
30 /
31 RCL O
32 *
33 2
34 ST+ O
35 CLX
36 RCL O
37 /
38 CHS
39 STO M
40 X<>Y
41 ST+ Y
42 GTO 01
43 LBL 02
44 ST+ X
45 PI
46 SQRT
47 /
48 FS?C 06
49 CHS
50 GTO 09
51 LBL 03
52 SF 05
53 CLX
54 4
55 X<=Y?
56 GTO 05
57 24
58 .
59 LBL 04
60 RCL M
61 +
62 RCL Y
63 X<>Y
64 ST+ X
65 /
66 DSE Y
67 GTO 04
68 RCL M
69 +
70 1/X
71 GTO 08
72 LBL 05
73 E99
74 X<> M
75 STO N
76 ENTER
77 1/X
78 STO O
79 CLX
80 SIGN
81 RCL M
82 RCL Y
83 LBL 06
84 X=Y?
85 GTO 07
86 R^
87 R^
88 ENTER
89 ENTER
90 2
91 /
92 ENTER
93 X<> M
94 ST/ M
95 X<>Y
96 X<> N
97 ST/ N
98 R^
99 ST+ M
100 ST+ N
101 R^
102 RCL M
103 RCL N
104 /
105 ST* O
106 E
107 ST+ Z
108 GTO 06
109 LBL 07
110 R^
111 STO M
112 RCL O
113 LBL 08
114 PI
115 SQRT
116 /
117 STO Y
118 RCL M
119 X^2
120 CHS
121 E^X
122 *
123 FC?C 06
124 GTO 09
125 2
126 -
127 CHS
128 LBL 09
129 ENTER
130 ENTER
131 E
132 -
133 CHS
134 FC?C 05
135 X<>Y
136 END

The function E1 also uses the synthetic text line "\05\77\21\56\64\99\99" (bytes F7 05 77 21 56 64 99 99), RCL M, which stores Euler's constant 0.5772156649 in the X register. The HP41 code for E1(x) [x XEQ E1]:

01 LBL "E1"
02 STO M
03 5
04 LN
05 X<=Y?
06 GTO 03
07 CLX
08 32
09 CHS
10 X>Y?
11 GTO 08
12 RCL M
13 CHS
14 STO M
15 STO O
16 E
17 STO N
18 .
19 RCL O
20 LBL 01
21 X=Y?
22 GTO 02
23 RCL M
24 RCL N
25 *
26 ISG N
27 ""     ;F0
28 RCL N
29 X^2
30 /
31 ST* O
32 RDN
33 RCL X
34 RCL O
35 +
36 GTO 01
37 LBL 02
38 CHS
39 RCL M
40 ABS
41 LN
42 -
43 " w!Vd "     ;"\05\77\21\56\64\99\99"
44 RCL M
45 -
46 RTN
47 LBL 03
48 CLX
49 3
50 X<=Y?
51 GTO 06
52 27
53 E
54 LBL 04
55 SF 06
56 LBL 05
57 RCL M
58 *
59 RCL Y
60 X<>Y
61 /
62 E
63 +
64 FS?C 06
65 GTO 05
66 DSE Y
67 GTO 04
68 RCL M
69 E^X
70 LASTX
71 *
72 *
73 1/X
74 RTN
75 LBL 06
76 RCL M
77 CHS
78 E^X
79 RCL M
80 /
81 STO O
82 E99
83 X<> M
84 E
85 ENTER
86 STO N
87 RCL M
88 LBL 07
89 X=Y?
90 GTO 10
91 R^
92 R^
93 ENTER
94 ENTER
95 RCL T
96 /
97 STO L
98 RDN
99 SF 06
100 LBL 08
101 LASTX
102 ENTER
103 X<> M
104 ST/ M
105 X<>Y
106 X<> N
107 ST/ N
108 CLX
109 E
110 ST+ M
111 ST+ N
112 R^
113 R^
114 RCL M
115 LASTX
116 X<> N
117 /
118 ST* O
119 LASTX
120 X<> N
121 STO L
122 CLX
123 E
124 X=Y?
125 GTO 10
126 FC?C 06
127 GTO 09
128 R^
129 R^
130 GTO 08
131 LBL 09
132 ISG Z
133 ""     ;F0
134 GTO 07
135 LBL 10
136 RCL O
137 RTN
138 LBL 11
139 13
140 STO N
141 RCL M
142 CHS
143 STO M
144 1/X
145 STO O
146 E
147 GTO 13
148 LBL 12
149 RCL N
150 *
151 RCL O
152 *
153 E
154 +
155 LBL 13
156 DSE N
157 GTO 12
158 LBL 13
159 RCL M
160 E^X
161 *
162 RCL M
163 /
164 CHS
165 END

The HP41 code for En(x) [n ENTER x XEQ EN]:

01 LBL "EN"
02 X<>Y
03 X!=0?
04 GTO 00
05 X<>Y
06 E^X
07 LASTX
08 *
09 1/X
10 RTN
11 LBL 00
12 E
13 X!=Y?
14 GTO 00
15 RCL Z
16 XEQ "E1"
17 RTN
18 LBL 00
19 RDN
20 X<>Y
21 5
22 LN
23 X<=Y?
24 GTO 04
25 RDN
26 STO 00
27 X<>Y
28 STO 01
29 XEQ "PSI"
30 RCL 00
31 STO M
32 LN
33 -
34 RCL M
35 CHS
36 RCL 01
37 E
38 STO O
39 -
40 CHS
41 STO N
42 CHS
43 Y^X
44 LASTX
45 FACT
46 /
47 *
48 STO 01
49 .
50 RCL N
51 1/X
52 RCL Y
53 LBL 01
54 X=Y?
55 GTO 03
56 CLX
57 SIGN
58 ST+ Z
59 ST+ N
60 CLX
61 RCL N
62 X!=0?
63 GTO 02
64 CLX
65 RCL M
66 CHS
67 RCL Z
68 /
69 ST* O
70 CLX
71 GTO 01
72 LBL 02
73 CLX
74 RCL M
75 CHS
76 RCL Z
77 /
78 RCL N
79 /
80 ST* O
81 CLX
82 RCL Y
83 RCL O
84 +
85 X<>Y
86 RCL N
87 ST* O
88 RDN
89 GTO 01
90 LBL 03
91 CHS
92 RCL 01
93 +
94 RTN
95 LBL 04
96 RDN
97 CHS
98 E^X
99 LASTX
100 CHS
101 /
102 STO O
103 CLX
104 LASTX
105 X<>Y
106 STO L
107 CLX
108 E
109 E99
110 STO M
111 RCL Y
112 STO N
113 LBL 05
114 X=Y?
115 GTO 08
116 R^
117 R^
118 X<> L
119 ENTER
120 ENTER
121 RCL T
122 ST/ Y
123 RDN
124 ENTER
125 SF 06
126 LBL 06
127 X<> M
128 ST/ M
129 X<>Y
130 X<> N
131 ST/ N
132 CLX
133 E
134 ST+ M
135 ST+ N
136 R^
137 R^
138 RCL M
139 ST* O
140 RCL N
141 ST/ O
142 X=Y?
143 GTO 08
144 FC?C 06
145 GTO 07
146 R^
147 R^
148 X<> L
149 ENTER
150 ENTER
151 RCL T
152 ST/ Y
153 RDN
154 ENTER
155 GTO 06
156 LBL 07
157 ISG Z
158 ""     ;F0
159 ISG L
160 ""     ;F0
161 GTO 05
162 LBL 08
163 RCL O
164 END

Edited: 23 Sept 2009, 3:03 p.m.

Password:

[ Return to the Message Index ]

Go back to the main exhibit hall