Lentz's algorithm

In mathematics, Lentz's algorithm is an algorithm to evaluate continued fractions, and was originally devised to compute tables of spherical Bessel functions.^{[1] [2] [full citation needed ][3]}

The version often employed now is the simplification due to Thompson and Barnett.^[4]

The idea was introduced in 1973 by William J. Lentz^[1] and was simplified by him in 1982.^[5] 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.^[6] Simpler Improvements to overcome unwanted zero terms include an altered recurrence relation^[7] 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.^[4] , or Lentz's simplification.

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.^[6] The algorithm suggested that it is possible to terminate the evaluation of continued fractions when ${\displaystyle |f_{j}-f_{j-1}|}${\displaystyle |f_{j}-f_{j-1}|} is relatively small.^[8]

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.^[9]

For continued fractions, the definitive standard notation is found under "Elementary Analytical Methods", page 19 and throughout the text for each function.^[10]

${\displaystyle f_{0}=b_{0}}${\displaystyle f_{0}=b_{0}}

${\displaystyle f_{1}=b_{0}+{\frac {a_{1}}{b_{1}}}}${\displaystyle f_{1}=b_{0}+{\frac {a_{1}}{b_{1}}}}

${\displaystyle f_{2}=b_{0}+{\frac {a_{1}}{b_{1}+{\frac {a_{2}}{b_{2}}}}}}${\displaystyle f_{2}=b_{0}+{\frac {a_{1}}{b_{1}+{\frac {a_{2}}{b_{2}}}}}}

${\displaystyle f_{3}=b_{0}+{\frac {a_{1}}{b_{1}+{\frac {a_{2}}{b_{2}+{\frac {a_{3}}{b_{3}}}}}}}}${\displaystyle f_{3}=b_{0}+{\frac {a_{1}}{b_{1}+{\frac {a_{2}}{b_{2}+{\frac {a_{3}}{b_{3}}}}}}}}

etc., or using the big-K notation, if

${\displaystyle f_{n}=b_{0}+{\underset {j=1}{\overset {n}{\operatorname {K} }}}{\frac {a_{j}}{b_{j}}}}${\displaystyle f_{n}=b_{0}+{\underset {j=1}{\overset {n}{\operatorname {K} }}}{\frac {a_{j}}{b_{j}}}}

is the ${\displaystyle n}${\displaystyle n}th convergent to ${\displaystyle f}${\displaystyle f} then

${\displaystyle f_{n}={\frac {A_{n}}{B_{n}}}}${\displaystyle f_{n}={\frac {A_{n}}{B_{n}}}}

where ${\displaystyle A_{n}}${\displaystyle A_{n}} and ${\displaystyle B_{n}}${\displaystyle B_{n}} are given by the Wallis-Euler recurrence relations

${\displaystyle {\begin{aligned}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{aligned}}}${\displaystyle {\begin{aligned}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{aligned}}}

Lentz's method defines

${\displaystyle C_{n}={\frac {A_{n}}{A_{n-1}}}}${\displaystyle C_{n}={\frac {A_{n}}{A_{n-1}}}}

${\displaystyle D_{n}={\frac {B_{n-1}}{B_{n}}}}${\displaystyle D_{n}={\frac {B_{n-1}}{B_{n}}}}

so that the ${\displaystyle n}${\displaystyle n}th convergent is ${\displaystyle f_{n}=C_{n}D_{n}f_{n-1}}${\displaystyle f_{n}=C_{n}D_{n}f_{n-1}} with ${\textstyle f_{0}={\frac {A_{0}}{B_{0}}}=b_{0}}${\textstyle f_{0}={\frac {A_{0}}{B_{0}}}=b_{0}} and uses the recurrence relations

${\displaystyle {\begin{aligned}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}&={\frac {1}{b_{n}+a_{n}D_{n-1}}}\end{aligned}}}${\displaystyle {\begin{aligned}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}&={\frac {1}{b_{n}+a_{n}D_{n-1}}}\end{aligned}}}

When the product ${\displaystyle C_{n}D_{n}}${\displaystyle C_{n}D_{n}} reaches unity with increasing ${\displaystyle n}${\displaystyle n} to the accuracy of the computer, then ${\displaystyle f_{n}}${\displaystyle f_{n}} has converged to ${\displaystyle f}${\displaystyle f}.^[11]

Lentz's algorithm has the advantage of side-stepping an inconvenience of the Wallis-Euler relations, namely that the numerators ${\displaystyle A_{n}}${\displaystyle A_{n}} and denominators ${\displaystyle B_{n}}${\displaystyle B_{n}} are prone to grow or diminish very rapidly with increasing ${\displaystyle n}${\displaystyle n}. In direct numerical application of the Wallis-Euler relations, this means that ${\displaystyle A_{n-1}}${\displaystyle A_{n-1}}, ${\displaystyle A_{n-2}}${\displaystyle A_{n-2}}, ${\displaystyle B_{n-1}}${\displaystyle B_{n-1}}, ${\displaystyle B_{n-2}}${\displaystyle B_{n-2}} must be periodically checked and rescaled to avoid floating-point overflow or underflow.^[11]

In Lentz's original algorithm, it can happen that ${\displaystyle C_{n}=0}${\displaystyle C_{n}=0}, 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 ${\displaystyle C_{n}=\varepsilon }${\displaystyle C_{n}=\varepsilon } for some sufficiently small ${\displaystyle \varepsilon }${\displaystyle \varepsilon }. This gives ${\textstyle C_{n+1}=b_{n+1}+{\frac {a_{n+1}}{\varepsilon }}={\frac {a_{n+1}}{\varepsilon }}}${\textstyle C_{n+1}=b_{n+1}+{\frac {a_{n+1}}{\varepsilon }}={\frac {a_{n+1}}{\varepsilon }}} to within floating-point precision, and the product ${\displaystyle C_{n}C_{n+1}=a_{n+1}}${\displaystyle C_{n}C_{n+1}=a_{n+1}} irrespective of the precise value of ε. Accordingly, the value of ${\displaystyle f_{0}=C_{0}=b_{0}}${\displaystyle f_{0}=C_{0}=b_{0}} is also set to ${\displaystyle \varepsilon }${\displaystyle \varepsilon } in the case of ${\displaystyle b_{0}=0}${\displaystyle b_{0}=0}.

Similarly, if the denominator in ${\textstyle D_{n}={\frac {1}{b_{n}+a_{n}D_{n-1}}}}${\textstyle D_{n}={\frac {1}{b_{n}+a_{n}D_{n-1}}}} is zero, then setting ${\textstyle D_{n}={\frac {1}{\varepsilon }}}${\textstyle D_{n}={\frac {1}{\varepsilon }}} for small enough ${\displaystyle \varepsilon }${\displaystyle \varepsilon } gives ${\textstyle D_{n}D_{n+1}={\frac {1}{a_{n+1}}}}${\textstyle D_{n}D_{n+1}={\frac {1}{a_{n+1}}}} irrespective of the value of ${\displaystyle \varepsilon }${\displaystyle \varepsilon }.^{[4] [11]}

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.^[12] As an example, it was applied to evaluate exponential integral functions. This application was then called modified Lentz algorithm.^[13] 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.^[14] 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.

• Special hypergeometric functions

• All articles with incomplete citations
• Articles with incomplete citations from September 2025