5.2 Evaluation of Continued Fractions 171 where the denominator in (5.2.6)approaches zero,so that Di and Afi are very large.The next Afj+will typically cancel this large change,but with loss of accuracy in the numerical running sum of the fi's.It is awkward to program around this,so Steed's method can be recommended only for cases where you know in advance that no denominator can vanish.We will use it for a special purpose in the routine bessik (86.7). The best general method for evaluating continued fractions seems to be the modified Lentz's method [61.The need for rescaling intermediate results is avoided by using both the ratios Cj=Aj/Aj-1;Dj=Bj-1/Bj (5.2.8) and calculating f;by f方=fi-1C5D5 (5.2.9) 满 From equation(5.2.5),one easily shows that the ratios satisfy the recurrence relations D=1/(b+aD5-1),C3=b+a/Ci-1 (5.2.10 令 In this algorithm there is the danger that the denominator in the expression for Dj, or the quantity C;itself,might approach zero.Either of these conditions invalidates (5.2.10).However,Thompson and Barnett [5]show how to modify Lentz's algorithm to fix this:Just shift the offending term by a small amount,e.g.,10-30.If you work through a cycle of the algorithm with this prescription,you will see that fj+ 9」 is accurately calculated. In detail,the modified Lentz's algorithm is this: Set fo bo;if bo =0 set fo tiny ·Set Co=fo: to dir 。Set Do=0. ·Forj=1,2, OF SCIENTIFIC COMPUTING (ISBN Set Dj=bj ajDj-1. If Di=0,set Di tiny. Set Cj=bj aj/Cj-1. If Cj=0 set Cj=tiny. Set Di=1/Di v@cambridge.org 1988-1992 by Numerical Recipes Set△，=CiDj 12-:621-43106-50 Set fi=fi-1△i Ifl△j-l<eps then exit. (outside Here eps is your floating-point precision,say 10-7or 10-15.The parameter tiny North Software. should be less than typical values of eps say 10-30. The above algorithm assumes that you can terminate the evaluation of the continued fraction when f;-fi-1 is sufficiently small.This is usually the case, but by no means guaranteed.Jones [7]gives a list of theorems that can be used to justify this termination criterion for various kinds of continued fractions. There is at present no rigorous analysis oferror propagation in Lentz's algorithm. However,empirical tests suggest that it is at least as good as other methods.5.2 Evaluation of Continued Fractions 171 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). where the denominator in (5.2.6) approaches zero, so that D j and ∆fj are very large. The next ∆fj+1 will typically cancel this large change, but with loss of accuracy in the numerical running sum of the fj’s. It is awkward to program around this, so Steed’s method can be recommended only for cases where you know in advance that no denominator can vanish. We will use it for a special purpose in the routine bessik (§6.7). The best general method for evaluating continued fractions seems to be the modified Lentz’s method [6]. The need for rescaling intermediate results is avoided by using both the ratios Cj = Aj/Aj−1, Dj = Bj−1/Bj (5.2.8) and calculating fj by fj = fj−1CjDj (5.2.9) From equation (5.2.5), one easily shows that the ratios satisfy the recurrence relations Dj = 1/(bj + ajDj−1), Cj = bj + aj/Cj−1 (5.2.10) In this algorithm there is the danger that the denominator in the expression for D j , or the quantity Cj itself, might approach zero. Either of these conditions invalidates (5.2.10). However, Thompson and Barnett [5] show how to modify Lentz’s algorithm to fix this: Just shift the offending term by a small amount, e.g., 10 −30. If you work through a cycle of the algorithm with this prescription, you will see that f j+1 is accurately calculated. In detail, the modified Lentz’s algorithm is this: • Set f0 = b0; if b0 = 0 set f0 = tiny. • Set C0 = f0. • Set D0 = 0. • For j = 1, 2,... Set Dj = bj + ajDj−1. If Dj = 0, set Dj = tiny. Set Cj = bj + aj/Cj−1. If Cj = 0 set Cj = tiny. Set Dj = 1/Dj. Set ∆j = CjDj. Set fj = fj−1∆j . If |∆j − 1| < eps then exit. Here eps is your floating-point precision, say 10−7 or 10−15. The parameter tiny should be less than typical values of eps|bj|, say 10−30. The above algorithm assumes that you can terminate the evaluation of the continued fraction when |fj − fj−1| is sufficiently small. This is usually the case, but by no means guaranteed. Jones [7] gives a list of theorems that can be used to justify this termination criterion for various kinds of continued fractions. There is at present no rigorous analysis of error propagationin Lentz’s algorithm. However, empirical tests suggest that it is at least as good as other methods