In mathematics, '''Lentz's algorithm''' is an algorithm to evaluate continued fractions,<ref name="numerical-recipes-c++">{{cite book|title=Numerical Recipes in C++| pages=177–179|isbn= 0-521-75033-4}}</ref>{{full citation needed|date=September 2025}} and was originally devised to compute tables of spherical Bessel functions.<ref name=":0">{{Cite report|last=Lentz|first=W. J.|date=September 1973|title=A Method of Computing Spherical Bessel Functions of Complex Argument with Tables|url=https://apps.dtic.mil/sti/pdfs/AD0767223.pdf|publisher=Atmospheric Sciences Laboratory, US Army Electronics Command|location=White Sands Missile Range, New Mexico|type=Research and Development Technical Report ECOM-5509 }}</ref><ref>{{cite journal |last1=Lentz |first1=William |title=Continued fraction calculation of spherical Bessel functions |journal=Computers in Physics |date=July 1990 |volume=4 |issue=4 |pages=403–407 |DOI=10.1063/1.168382}}</ref>
The version often employed now is the simplification due to Thompson and Barnett.<ref name="Thompson-and-Barnett" />
== History == The idea was introduced in 1973 by William J. Lentz<ref name=":0" /> and was simplified by him in 1982.<ref>{{Cite book |last=Lentz |first=W. J. |title=A Simplification of Lentz's Algorithm |date=August 1982 |publisher=Defense Technical Information Center |oclc=227549426}}</ref> Lentz suggested that calculating ratios of spherical Bessel functions of complex arguments over a wide range of values can be difficult. He developed a new continued fraction technique for calculating the ratios of spherical Bessel functions of consecutive order. This method was an improvement compared to other methods because it starts from the beginning of the continued fraction rather than the tail, had a built-in check for convergence, and is numerically stable. The original algorithm uses algebra to bypass a zero or near-zero independently in either the numerator or denominator.<ref name="Lentz 668–671">{{Cite journal |last=Lentz |first=William J. |date=1976-03-01 |title=Generating Bessel functions in Mie scattering calculations using continued fractions |journal=Applied Optics |volume=15 |issue=3 |pages=668–671 |bibcode=1976ApOpt..15..668L |doi=10.1364/ao.15.000668 |issn=0003-6935 |pmid=20165036}}</ref> Simpler Improvements to overcome unwanted zero terms include an altered recurrence relation<ref>{{Cite journal |last1=Jaaskelainen |first1=T. |last2=Ruuskanen |first2=J. |date=1981-10-01 |title=Note on Lentz's algorithm |journal=Applied Optics |volume=20 |issue=19 |pages=3289–3290 |bibcode=1981ApOpt..20.3289J |doi=10.1364/ao.20.003289 |issn=0003-6935 |pmid=20333144}}</ref> suggested by Jaaskelainen and Ruuskanen in 1981 or a simple shift of the denominator by a very small number as suggested by Thompson and Barnett in 1986.,<ref name="Thompson-and-Barnett">{{Cite journal |last1=Thompson |first1=I.J. |last2=Barnett |first2=A.R. |date=1986 |title=Coulomb and Bessel functions of complex arguments and order |journal=Journal of Computational Physics |volume=64 |issue=2 |pages=490–509 |bibcode=1986JCoPh..64..490T |doi=10.1016/0021-9991(86)90046-x |issn=0021-9991}}</ref> or Lentz's simplification.
== Initial work == This theory was initially motivated by Lentz's need for accurate calculation of ratios of spherical Bessel functions of consecutive order and complex argument necessary for Mie scattering. He created a new continued fraction algorithm that starts from the beginning of the continued fraction and not at the tail-end. This eliminates guessing how many terms of the continued fraction are needed for convergence. In addition, continued fraction representations for both ratios of Bessel functions and spherical Bessel functions of consecutive order themselves can be computed with Lentz's algorithm.<ref name="Lentz 668–671"/> The algorithm suggested that it is possible to terminate the evaluation of continued fractions when <math>|f_j-f_{j-1} |</math> is relatively small.<ref>{{Cite book|last1=Masmoudi|first1=Atef|last2=Bouhlel|first2=Med Salim|last3=Puech|first3=William|title=2012 6th International Conference on Sciences of Electronics, Technologies of Information and Telecommunications (SETIT) |chapter=Image encryption using chaotic standard map and engle continued fractions map |date=March 2012|chapter-url=http://dx.doi.org/10.1109/setit.2012.6481959|pages=474–480 |publisher=IEEE|doi=10.1109/setit.2012.6481959|isbn=978-1-4673-1658-3 |s2cid=15380706 }}</ref>
== Algorithm == Lentz's algorithm is based on the Wallis-Euler relations. John Wallis independently rediscovered the recursion relations about 500 years after the Indian mathematician Bhas-Cara II.<ref>{{cite book |last1=Brezinski |first1=Claude |title=History of Continued Fractions and Pade Approximants |date=1991 |publisher=Springer-Verlag |location=Berlin Heidelberg New York |isbn=3-540-15286-5 |pages=32–33}}</ref>
For continued fractions, the definitive standard notation is found under "Elementary Analytical Methods", page 19 and throughout the text for each function.<ref>{{cite book |last1=Abramowitz |first1=Milton |title=Handbook of Mathematical Functions |date=1972 |publisher=Dover |location=New York |isbn=978-0486612720|edition=9th}}</ref> :<math>f_0 = b_0</math> :<math>f_1 = b_0 + \frac{a_1}{b_1}</math> :<math>f_2 = b_0 + \frac{a_1}{b_1 + \frac{a_2}{b_2}}</math> :<math>f_3 = b_0 + \frac{a_1}{b_1 + \frac{a_2}{b_2 + \frac{a_3}{b_3}}}</math> etc., or using the big-K notation, if :<math>f_n = b_0 + \underset{j = 1}\overset{n}\operatorname{K}\frac{a_j}{b_j}</math> is the <math>n</math>th convergent to <math>f</math> then :<math>f_n = \frac{A_n}{B_n}</math> where <math>A_n</math> and <math>B_n</math> are given by the Wallis-Euler recurrence relations :<math> \begin{align} A_{-1} & = 1 & B_{-1} & = 0\\ A_0 & = b_0 & B_0 & = 1\\ A_n & = b_n A_{n-1} + a_n A_{n-2} & B_n & = b_n B_{n-1} + a_n B_{n-2} \end{align} </math>
Lentz's method defines :<math>C_n = \frac{A_n}{A_{n - 1}}</math> :<math>D_n = \frac{B_{n - 1}}{B_n}</math> so that the <math>n</math>th convergent is <math>f_n = C_n D_n f_{n - 1}</math> with <math display="inline">f_0 = \frac{A_0}{B_0} = b_0</math> and uses the recurrence relations :<math> \begin{align} C_0 & = \frac{A_0}{A_{-1}} = b_0 & D_0 & = \frac{B_{-1}}{B_0} = 0\\ C_n & = b_n + \frac{a_n}{C_{n-1}} & D_n & = \frac1{b_n + a_n D_{n-1}} \end{align} </math>
When the product <math>C_n D_n</math> reaches unity with increasing <math>n</math> to the accuracy of the computer, then <math>f_n</math> has converged to <math>f</math>.<ref name="Numerical-Recipes">{{Cite book |last1=Press |first1=W.H. |title=Numerical Recipes: The Art of Scientific Computing |last2=Teukolsky |first2=S.A. |last3=Vetterling |first3=W.T. |last4=Flannery |first4=B. P. |publisher=Cambridge University Press |year=2007 |edition=3rd |pages=207–208}}</ref>
Lentz's algorithm has the advantage of side-stepping an inconvenience of the Wallis-Euler relations, namely that the numerators <math>A_n</math> and denominators <math>B_n</math> are prone to grow or diminish very rapidly with increasing <math>n</math>. In direct numerical application of the Wallis-Euler relations, this means that <math>A_{n-1}</math>, <math>A_{n-2}</math>, <math>B_{n-1}</math>, <math>B_{n-2}</math> must be periodically checked and rescaled to avoid floating-point overflow or underflow.<ref name="Numerical-Recipes" />
== Thompson and Barnett modification == In Lentz's original algorithm, it can happen that <math>C_n = 0</math>, resulting in division by zero at the next step. A simpler method than that proposed by Lentz remedied the problem can be simply by setting <math>C_n = \varepsilon</math> for some sufficiently small <math>\varepsilon</math>. This gives <math display="inline">C_{n + 1} = b_{n + 1} + \frac{a_{n + 1}}{\varepsilon} = \frac{a_{n + 1}}{\varepsilon}</math> to within floating-point precision, and the product <math>C_n C_{n + 1} = a_{n + 1}</math> irrespective of the precise value of ε. Accordingly, the value of <math>f_0 = C_0 = b_0</math> is also set to <math>\varepsilon</math> in the case of <math>b_0 = 0</math>.
Similarly, if the denominator in <math display="inline">D_n = \frac{1}{b_n + a_n D_{n - 1}}</math> is zero, then setting <math display="inline">D_n = \frac{1}{\varepsilon}</math> for small enough <math>\varepsilon</math> gives <math display="inline">D_n D_{n + 1} = \frac{1}{a_{n + 1}}</math> irrespective of the value of <math>\varepsilon</math>.<ref name="Thompson-and-Barnett" /><ref name="Numerical-Recipes" />
== Applications == Lentz's algorithm was used widely in the late twentieth century. It was suggested that it doesn't have any rigorous analysis of error propagation. However, a few empirical tests suggest that it's at least as good as the other methods.<ref>{{Cite book |last1=Press |first1=W.H. |title=Numerical Recipes in Fortran, The Art of Scientific Computing|last2=Teukolsky |first2=S.A. |last3=Vetterling |first3=W.T. |last4=Flannery |first4=B. P. |publisher=Cambridge University Press |year=1992 |edition=2nd |page=165}}</ref> As an example, it was applied to evaluate exponential integral functions. This application was then called modified Lentz algorithm.<ref>{{Cite journal|last1=Press|first1=William H.|last2=Teukolsky|first2=Saul A.|date=1988|title=Evaluating Continued Fractions and Computing Exponential Integrals|journal=Computers in Physics|volume=2|issue=5|pages=88|doi=10.1063/1.4822777|bibcode=1988ComPh...2...88P |issn=0894-1866|doi-access=free}}</ref> It's also stated that the Lentz algorithm is not applicable for every calculation, and convergence can be quite rapid for some continued fractions and vice versa for others.<ref>{{Cite journal|last1=Wand|first1=Matt P.|last2=Ormerod|first2=John T.|date=2012-09-18|title=Continued fraction enhancement of Bayesian computing|url=http://dx.doi.org/10.1002/sta4.4|journal=Stat|volume=1|issue=1|pages=31–41|doi=10.1002/sta4.4|pmid=22533111 |s2cid=119636237 |issn=2049-1573|url-access=subscription}}</ref> If the terms of the continued fraction are constant or increasing, it is stable, but if the consecutive terms are decreasing, the algorithm may lose accuracy. This is because the forward recursion from which the algorithm is derived is unstable for decreasing terms.
==References== {{reflist}}
Category:Special hypergeometric functions