5
$\begingroup$

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$

Improve this question
$\endgroup$
  • $\begingroup$ The Seidel-Stern theorem may be of use here. $\endgroup$ – Simply Beautiful Art Jan 18 '17 at 2:22
4
+50
$\begingroup$

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

Improve this answer
$\endgroup$
  • $\begingroup$ Ah! I think I figured out where the bad rounding was happening. And it makes sense that, given the exponential nature of $\exp$ (duh, lol) rounding errors would cause the result to be skewed by a, well, exponential amount. $\endgroup$ – Eric Lagergren Jan 20 '17 at 19:00
  • $\begingroup$ admittedly this is another question, but is there a way I could determine how many decimal digits of precision I'd need to compute the correct result? $\endgroup$ – Eric Lagergren Jan 20 '17 at 19:01
  • $\begingroup$ This other question is very interesting, but quite difficult to answer. How much reliable is the 500-digit precision software that you tested? $\endgroup$ – Anatoly Jan 20 '17 at 20:08
  • $\begingroup$ my software is reliable. Is there some sort of algorithm that will produce the required precision so that the $\exp$ algorithm doesn't have rounding errors? $\endgroup$ – Eric Lagergren Jan 20 '17 at 20:16
  • $\begingroup$ Just an idea: you could use your software to determine, for relatively small values of $z $, the precision level that gives a predefined level of approximation (or maximal tolerated error). For example, you could predefine that an absolute error $<0.5$ (that is to say, the result is right for all nondecimal digits) is tolerated. Then, launch your software starting with a standard 16-digits precision to calculate the absolute error for increasing values of $z $. As you get the first value of $z $ that violates the tolerance criterion, you can increase the precision to, say, $k $ digits. $\endgroup$ – Anatoly Jan 20 '17 at 20:35
0
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.