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$