Lentz's algorithm: Difference between revisions

Content deleted Content added
new article
 
Bumpf (talk | contribs)
m Algorithm: fix "big K" notation
 
(61 intermediate revisions by 21 users not shown)
Line 1:
In mathematics, '''Lentz's Algorithmalgorithm''' is usedan [[algorithm]] to calculateevaluate [[generalized continued fractionsfraction|continued fraction]]s, and presentwas originally devised to compute tables of spherical [[Bessel function|Bessel functions]]s.<ref name=":0">{{Cite journalreport|last=Lentz|first=W. J.|date=September 1973-09-01|title=A Method of Computing Spherical Bessel Functions of Complex Argument with Tables|url=httphttps://dxapps.doidtic.orgmil/10.21236sti/ad0767223pdfs/AD0767223.pdf|publisher=Atmospheric Sciences Laboratory, US Army Electronics Command|___location=FortWhite BelvoirSands Missile Range, VANew Mexico|type=Research and Development Technical Report ECOM-5509 }}</ref><ref name="numerical-recipes-c++">{{cite book|title=Numerical Recipes in C++| pages=177–179|isbn= 0 521 75033 4}}</ref>
 
The version usually employed now is due to Thompson and Barnett.<ref name="Thompson-and-Barnett" />
=== History ===
The idea was introduced more than thirty years ago by W.J. Lentz. Lentz suggested that calculating ratios of Spherical Bessel functions of complex arguments can be difficult. He developed a new continued fraction technique for calculating them. This method was an improvement compared to other methods because it eliminated errors on certain terms or provided zero as a result<ref>{{Cite book|last=J.|first=Lentz, W.|url=http://worldcat.org/oclc/227549426|title=A Simplification of Lentz's Algorithm.|date=1982-08|publisher=Defense Technical Information Center|oclc=227549426}}</ref>.
 
=== InitialHistory Working ===
ThisThe theoryidea was initiallyintroduced implementedin 1973 by William J. Lentz<ref name=":0" /> and was simplified by him in 1982.<ref>{{Cite book|last=J.|first=Lentz, W.|url=http://worldcat.org/oclc/227549426|title=A Simplification of Lentz's anotherAlgorithm.|date=August research1982|publisher=Defense whenTechnical heInformation calculatedCenter|oclc=227549426}}</ref> Lentz suggested that calculating ratios of spherical Bessel functionfunctions necessaryof forcomplex [[Miearguments scattering]]can be difficult. He demonstrateddeveloped thata thenew algorithmcontinued uses afraction technique involvingfor calculating the evaluationratios continuedof fractionsspherical thatBessel functions of consecutive order. This method was an improvement compared to other methods because it startsstarted from the beginning andof notthe atcontinued thefraction tail.rather Inthan additionthe tail, thathad continueda fractionbuilt-in representationscheck for bothconvergence, ratiosand ofwas Besselnumerically functionsstable. and sphericalThe Besseloriginal functionsalgorithm ofuses consecutivealgebra orderto canbypass bea presentedzero within either the Lentznumerator algorithmor 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|url=http://dx.doi.org/10.1364/ao.15.000668|journal=Applied Optics|volume=15|issue=3|pages=668668–671|doi=10.1364/ao.15.000668|pmid=20165036 |bibcode=1976ApOpt..15..668L |issn=0003-6935|url-access=subscription}}</ref>. The algorithmSimpler suggested that it is possibleImprovements to terminateovercome theunwanted evaluationzero ofterms continuedinclude fractionsan whenaltered recurrence relation<mathref>{{Cite journal|f_jlast1=Jaaskelainen|first1=T.|last2=Ruuskanen|first2=J.|date=1981-f_(j10-1)01|title=Note on Lentz's algorithm|url=http://dx.doi.org/10.1364/ao.20.003289|journal=Applied Optics|volume=20|issue=19|pages=3289–3290|doi=10.1364/ao.20.003289|pmid=20333144 |bibcode=1981ApOpt..20.3289J |issn=0003-6935|url-access=subscription}}</mathref> issuggested relativelyby 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|lastlast1=MasmoudiThompson|firstfirst1=AtefI.J.|last2=BouhlelBarnett|first2=Med Salim|last3=Puech|first3=WilliamA.R.|date=2012-031986|title=ImageCoulomb encryptionand usingBessel chaoticfunctions standardof complex maparguments and engle continued fractions maporder|url=http://dx.doi.org/10.11091016/setit.2012.64819590021-9991(86)90046-x|journal=2012 6th International Conference on SciencesJournal of Electronics, TechnologiesComputational of Information and Telecommunications (SETIT)Physics|publishervolume=IEEE64|issue=2|pages=490–509|doi=10.11091016/setit0021-9991(86)90046-x|bibcode=1986JCoPh.2012.648195964..490T |issn=0021-9991|url-access=subscription}}</ref>.
 
=== ApplicationsInitial work ===
This theory was initially motivated by Lentz's need for accurate calculation of ratios of spherical Bessel function 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>
Lentz's algorithm was used widely in the late 1900s. It was suggested that it doesn't have any rigorous analysis of error propagation. However, a few empirical tests suggest that it's almost as good as the other methods. As an example, it was applied to evaluate exponential integral functions. This application was then called modified Lentz algorithm<ref>{{Cite journal|last=Press|first=William H.|last2=Teukolsky|first2=Saul A.|date=1988|title=Evaluating Continued Fractions and Computing Exponential Integrals|url=http://dx.doi.org/10.1063/1.4822777|journal=Computers in Physics|volume=2|issue=5|pages=88|doi=10.1063/1.4822777|issn=0894-1866}}</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|last=Wand|first=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|issn=2049-1573}}</ref>
----
 
== Algorithm ==
<references />
Lentz's algorithm is based on the [[Wallis-Euler relations]]. If
 
:<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 [[Generalized continued fraction#Notation|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>{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} & = \frac{1}{{b}_{n} + {a}_{n} {D}_{n-1}}
\end{align}
</math>
 
When the product <math>{C}_{n} {D}_{n}</math> approaches unity with increasing <math>n</math>, it is hoped that <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. The problem can be remedied simply by setting <math>{C}_{n} = \varepsilon</math> for some sufficiently small <math>\varepsilon</math>. This gives <math>{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>{D}_{n} = \frac{1}{{b}_{n} + {a}_{n} {D}_{n - 1}}</math> is zero, then setting <math>{D}_{n} = \frac{1}{\varepsilon}</math> for small enough <math>\varepsilon</math> gives <math>{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 1900stwentieth century. It was suggested that it doesn't have any rigorous analysis of error propagation. However, a few empirical tests suggest that it's almostat 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|lastlast1=Press|firstfirst1=William H.|last2=Teukolsky|first2=Saul A.|date=1988|title=Evaluating Continued Fractions and Computing Exponential Integrals|url=http://dx.doi.org/10.1063/1.4822777|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|lastlast1=Wand|firstfirst1=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>
 
==References==
{{reflist}}
 
[[Category:Special hypergeometric functions]]