I am writing a program which will compute $\exp(z)$.
Originally I used the Taylor series, which worked fine. However, continued fractions can converge more quickly than power series, so I decided to go that route.
I found this continued fraction in multiple places. It's from The Application of Continued Fractions and Their Generalizations to Problems in Approximation Theory by A. N. Khovanskii (1963), pg 114.
$${e^{z}=1+{\cfrac {2z}{2-z+{\cfrac {z^{2}}{6+{\cfrac {z^{2}}{10+{\cfrac {z^{2}}{14+\ddots }}}}}}}}}$$
It can be represented as
$${e^{z}=1+\cfrac{2z}{2-z +}\cfrac{z^2/6}{1 +}\sum_{m=3}^{\infty}\left({\cfrac{{a_m}^{z^2}}{1}}\right)}$$
(Sorry, first time using Mathjax.)
Annie A.M. Cuyt, Vigdis Petersen, Brigitte Verdonk, Haakon Waadeland, and William B. Jones. 2008. Handbook of Continued Fractions for Special Functions (1 ed.). Springer Publishing Company, Incorporated, pg 194.
I'm computing the continued fraction using the modified Lentz algorithm from Numerical Recipes in C: The Art of Scientific Computing (ISBN 0-521-43105-5), pg 171.
The issue: it only works for a small set of numbers. (From what I can tell, [-30, 30].)
So, my question: is this expected? I'm relatively new to continued fractions, so while I think I grasp them, I'm not entirely sure.
Given a "generator", C++'s boost library can compute continued fractions. Essentially each call to the "generator" returns the next term in the CF.
Here's what I used (where $z$ is the input and $m$ is the current term index):
- Term 0:
- A: $0$
- B: $1$
- Term 1:
- A: $2z$
- B: $2-z$
- Term 2:
- A: $z^2$
- B: $1$
- Term 3..N:
- A: $1 / (4 * (2m - 3) * (2m - 1))) * z^2$
- B: $1$
Given $z = 38.5$, the Lentz algorithm provides the following (each line is $f_j$):
-1.1095890410958904
1.3657233326736593
-1.8636200592602417
2.81700268519061
-4.711808282946167
8.70960008496205
-17.765408200530924
39.92118587023142
-98.65273729066934
267.59129323043885
-795.1270569517195
2582.9903072474267
-9154.51726459279
35324.19248620738
-148091.74385797494
673154.6572372171
-3.310843427214524e+06
1.758455218773504e+07
-1.0065603099002995e+08
6.197709770714103e+08
-4.097272719854685e+09
2.9029669729927063e+10
-2.2004739456295242e+11
1.7812408854416082e+12
-1.5389163725064727e+13
1.40311430366499e+14
-1.483409360122725e+15
8.317237533312957e+15
2.273093139324771e+16
1.964661705307002e+16
1.9877542316211204e+16
1.985862955594415e+16
1.9860080303623144e+16
1.985997463174032e+16
1.9859981941770148e+16
1.9859981460940616e+16
1.9859981491045844e+16
1.9859981489249772e+16
1.9859981489351984e+16
1.9859981489346428e+16
1.9859981489346716e+16
1.98599814893467e+16
So, it converges at $1.98599814893467e+16$ when the actual answer is supposed to be: $\approx 52521552285925160$
This seems to be a computational problem. In particular, it may be a round-off representation error, which typically occurs when we attempt to represent a number using a finite string of digits. In this case, rounding errors even in the decimal digits of the first numbers in the progression shown in the OP finally results in a large error of convergence.
The continued fraction is correct and works well for any $z $ (also note that it is only a variation of other Euler's continued fractions giving $e^z $). However, using standard calculation tools that work with 16 digits of precision, it provides a very good approximation only until about $z=15$. For higher values of $z $, it begins to significantly underestimate the true value of $e^z$. Notably, this error does not change if we increase the length of the continued fraction (i.e. increasing $m $), supporting the hypothesis of its round-off nature.
For instance, calculating at a $16$-digit precision, for $z=15\,$ we have a good approximation (absolute error $0.0004\,\,\,$), but for $z=20\,$ the underestimation error is $\approx 4.9\,$, raising to $\approx 208585\,\,\,$ for $z=25\,$ and to $\approx 9218010394\,\,\,$ for $z=30\,$.
Interestingly, in these first values, the magnitude of the absolute error seems to growth exponentially: see for example this linear fitting by WA for the relation between $z $ and the log of the absolute underestimation error, calculated for the values $z=20\,$, $z=25\,$, and $z=30\,$ (further analysis of the logarithm of the error for higher values of $z $, however, seems to suggest a slightly lower growth rate).