Evaluating continued fractions

Question

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$

Answer 1

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).

Answer 2

English is not my mother tongue, but I feel I can make me understood:

A better way to calculate $e^z$ is using a a Matrix $M_k$:

$$ \begin{pmatrix} (2k+1)+z & (2k+1) \\ (2k+1) & (2k+1)-z \\ \end{pmatrix} $$

or in shorter notation $( (2k+1)+z, (2k+1); (2k+1), (2k+1)-z )$

If $z$ is rational (in practice, a real number used in computer calculations is always rational, maybe as approximation to an irrational number), for $e^{p/q}$ you can transform $M_k$ to

$( (2k+1)q+p, (2k+1)q, (2k+1)q, (2k+1)q-p )$

The Product of the $M_k$ series from $1...n$ gives a Matrix $(a,b; c,d)$, and if all $4$ numbers in this matrix are positive (which may take a somewhat longer series for a big $z=p/q$), $a/c$ and $b/d$ are approximations of $e^{p/q}$, the correct value always between them.

To limit the growth of numbers you can

(1) treat the Matrices as fractions, i.e. cancel down by factors common to all $4$ terms

(2) extract a simple continued fraction, i.e transform $(a,b; c,d)$ to $(c, d; a-n*c, b-n*d)$ with $n$ as the next term in the simple continued fraction. If you then use the simple continued fraction to produce a Matrix.

If you use a something (program, library, etc.) which provides integers of virtually unlimited length (in Java, Integer - with capital I! - can be as large as $2^{31}$ binary digits), you will get good results for almost all realistic scenarios. In the size of integers is limited (as standard types are), you will get soon at the limits.