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.