Fastest way to calculate $e^x$ up to arbitrary number of decimals?

What are other faster methods of calculating $e^x$ up to any number of decimals other than using the Taylor series formula?

-
I'm told the series converges rapidly and is well-conditioned. Have you tried that? –  Zhen Lin Jan 21 '11 at 17:33
Why not use the standard series expansion $\sum_i x^i/i!$ ? I'm not sure, but it may pay of to first convert to a value in $[0,1[$, e.g. to calculate $e^{10}$, calculate $(e^{1/2})^{20}$ –  Myself Jan 21 '11 at 17:34
I believe Haskell supports it. And many programmings languages have standard libraries which support it. E.g. Java has java.math.BigDecimal. –  Peter Taylor Jan 21 '11 at 18:10
I used the generalized continued fraction method (en.wikipedia.org/wiki/Generalized_continued_fraction) to find e upto large number of decimals, and it is way to fast than this series method... The only problem with is that I dont know how deep I should go to find e^x upto a particular decimal places... I usually have to find it for large value... So if someone can help me with this it would be great!! –  Guanidene Jan 22 '11 at 9:55
To evaluate $e^x$ for large values of $x$, use the idea of range reduction explained in my answer. You'll need to evaluate the constant $\ln 2$ to an accuracy corresponding to what your computation requires, but since it's a constant, you may only need to compute and store it once. I'll add a note on how to get $\ln 2$ to my answer. –  hardmath Jan 22 '11 at 14:39

The theoretically fastest way appears to be to use Newton iteration to reduce the problem to computing the logarithm function and then using an algorithm based on the arithmetic-geometric mean to compute the logarithm. See this wikipedia page. http://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations

In practice the Taylor series should work fine, given a good implementation. The following webpage http://numbers.computation.free.fr/Constants/constants.html has an example impletation of using the taylor series to compute e. They claim it took 0.02 seconds to compute e to a thousand decimals on a PentiumII, 450 MHz processor.

-

If I wanted to implement an arbitrary precision approximation to $e^x$ from "scratch" (building on an arbitrary precision numeric library such as GMP) with faster convergence than Taylor series/polynomial approximations, I would design around continued fraction/rational function approximations.

Just as the precision given by the well-known Taylor series can be extended as required by using additional terms, so too will the continued fraction expansion give any required precision by using additional convergents.

At first it may seem that evaluating the continued fraction approximation entails a significant disadvantage in comparison with the power series. Additional terms can be added to the partial sums of the power series, to obtain left-to-right sequential approximations. There is a way to achieve the same effect with continued fractions, namely evaluating the partial numerators and denominators into the convergents through a recurrence relation.

There are a couple of important ideas worth considering even if another method for evaluating the exponential function were to be used on a bounded interval such as $[0,1]$, as for example power series or even lookup tables(!). One is symmetry of the exponential function, and the other is a method of range reduction.

Although the exponential function is not even (or odd), it satisfies a relation:

$$e^{-x} = 1/e^x$$

that implies it is a simple rational transformation of an even function:

$$e^x = (1 + x*f(x^2))/(1 - x*f(x^2))$$

which may be solved for a power series or continued fraction approximation of f(z). This is the idea of symmetry as it applies to the exponential function. Not only does it reduce the evaluaton of the exponential function on the whole real line to just the positive (or the negative) half, it also "contracts" the number of terms needed for a given accuracy by half (retaining only $x^2$ terms in the expansion).

Although a continued fraction or Taylor series approximation to the exponential function may converge for arbitrary real (or complex) values), the rate of convergence is not uniform. The Taylor series converges faster the closer the argument is to zero. Thus a range reduction is important for expressing the exponential function at an arbitrary real value in terms of a value in some interval like $[0,1]$. For example, if the floating point arithmetic is binary (radix 2), it can be especially convenient to use the familiar law of exponents with multiples of $\ln 2$:

$$e^x = 2^k * e^r \ \text{where} \ x = r + k*\ln 2$$

which allows $e^x$ to be evaluated, up to a change in binary exponent, using an approximation that converges rapidly over $[0,\ln 2]$.

Combining these two ideas (symmetry, range reduction) the speed of convergence can be limited to the interval $[0, \ln 2/2]$. Limiting the interval of evaluation may allow you to determine in advance how many terms of the continued fraction or Taylor series expansion have to be retained to obtain a desired accuracy. When this is known the evaluation can be done more efficiently (backwards recursion for continued fractions or Horner's method for truncated Taylor series/polynomials) than if we were forced to continually introduced further terms until the desired accuracy is reached.

The "faster" continued fraction I had in mind is formula (11.1.2) here:

[Handbook of continued fractions for special functions (Cuyt et al, 2008)]

Their reference is to this modern classic:

The application of continued fractions and their generalizations to problems in approximation theory
A.N. Khovanskii, 1963 (P.Noordhoff)

A neat webMathematic-based site by Andreas Lauschke presents some ideas for accelerating convergence of continued fractions by "contractions". The IP address changes from time to time and cannot be used as a link within StackExchange, but you can Google for it:

[Convergence Acceleration with Canonical Contractions of Continued Fractions: Exponential function -- Andreas Lauschke Consulting]

I have some notes on his formulas (derived by contraction from the one noted above) if that would be helpful. Some related material was contributed to Wolfram Demonstrations.

Computing the constant $\ln 2$

Generations of math students have been introduced to the concept of conditional versus absolute convergence by the example of the alternating harmonic series:

$$\ln 2 = 1 - 1/2 + 1/3 - 1/4 + ...$$

Of course this series, derived from a power series expansion of $\ln x$ around $x = 1$, has such slow convergence that even if we combine consecutive pairs of terms:

$$\ln 2 = 1/2 + 1/12 + 1/30 + ...$$

the resulting absolutely convergent series is useless for obtaining arbitrary precision values of $\ln 2$. For convenience the first seven partial sum of this are:

0.50000000...
0.58333333...
0.61666666...
0.63452381...
0.64563492...
0.65870518...


Since $\ln 2$ is 0.69314718..., we have an error of about a third of a unit in the first decimal place. In other words not much more convergence than one decimal place correct.

It therefore makes a striking contrast with the nice convergence of a continued fraction expansion of the same value:

$$\ln 2 = \cfrac{1}{1 + \cfrac{1}{2 + \cfrac{1}{3 + \cfrac{4}{4 + \cfrac{4}{5 + \cfrac{9}{6 + \cfrac{9}{7 + ...}}}}}}}$$

The first seven convergents are:

0.66666666...
0.70000000...
0.69230769...
0.69333333...
0.69312169...
0.69315245...


The error here is about half a unit in the fifth decimal place.

-
Please explain further the "idea of symmetry" that allows the number of terms needed to be halved. I don't understand it. –  MikeM Oct 1 '13 at 17:22
@MikeM: Briefly the function $e^x$ can be expressed as $(1+g(x))/(1-g(x))$ where $g(x)$ is an odd function of $x$. Try replacing $g(x)$ with $g(-x)$ and you will see it implies the "symmetry" $e^{-x} = 1/e^x$. Now a power series expansion of $g(x)$ will have only odd powers, or $g(x) = x f(x^2)$. Thus a series expansion of $g(x)$ has only half as many nonzero terms as one for $e^x$. If you want a detailed comparison of op counts for equal accuracy, please post that as a Question so it can be addressed at length. –  hardmath Oct 2 '13 at 12:56
@MikeM: Note that this "economy" of using only even powers of $x$ can be realized by generalized continued fraction as well as power series expansions of $e^x$. –  hardmath Oct 2 '13 at 13:28
I have asked a new question here. –  MikeM Oct 2 '13 at 15:12

A very reasonable way to do this is to use

$e^x = 1 + x/1! + x^2/2! + x^3/3!...$

as it converges really fast.

A "trick" I read in some library for calculating $e^1$ is to calculate the series for $x=1/2^k$ where $k$ is the root of the amount of required bits, then you have to square the result k times to get your final number, it avoids that big/small numbers make the result too inaccurate.

Edit: Also this is the exact same topic like Iterative refinement algorithm for computing exp(x) with arbitrary precision

Maybe it should be merged / closed?

-
Well, this is a slow method, I know other methods which are more faster than this... I just want some other method –  Guanidene Jan 21 '11 at 17:51
I think the other question is significantly different as it asks for methods to increase the precision of computed values without starting over "from scratch". Here the emphasis seems to be on the speed with which a given accuracy can be achieved. –  hardmath Jan 21 '11 at 18:20

How is $x$ presented? If you want $e^x$ to (say) 100 significant figures, then you have to provide $x$ to (more or less) $100 + \log x$ significant figures. Are you using a multi-precision library?

-
Take x to be a simple number, say an integer... Now how are you going to proceed? –  Guanidene Jan 22 '11 at 10:35
@Guanidene: If $x$ is a (positive) integer, then I would compute $e$ to $100 + n + \log_{10} n$ significant figures, and raise it to the $n$th power using an exponentiation chain (see en.wikipedia.org/wiki/Addition-chain_exponentiation). –  TonyK Jan 24 '11 at 17:32

I don't know that this is necessarily the fastest way, but it's fun. If you want to know how it works, Google.

package com.akshor.pjt33.math;

import com.akshor.pjt33.math.symbalg.Rational;
import java.math.BigInteger;

public abstract class Spigot
{
/** The relevant components of the prefix matrix. */
Rational h00, h01=Rational.ZERO, h11=Rational.ONE;
BigInteger i=BigInteger.ZERO;
boolean verify=true;
final int base;
private StringBuilder sb=new StringBuilder();

public Spigot(BigInteger multiplier) {
this(multiplier, 10);
}

public Spigot(BigInteger multiplier, int base) {
h00=new Rational(multiplier, BigInteger.ONE);
this.base=base;
}

/**
* Advances this spigot algorithm, producing as many digits as each consumption permits.
* @param beta An upper bound (positive) on the magnitude of all subsequent term ratios.
* @return The number of digits output.
*/
if (verify && beta.signum()0) throw new IllegalArgumentException(Et_t+" exceeds bound "+beta);
h00=h00.mul(Et_t);

// Produce (if we can: without a converging bound we can't).
int rv=0;
return i;
}

protected abstract Rational termRatio();

public String toString() {
return sb.toString();
}

/**
* A spigot implementation which evaluates any hypergeometric sum for which you can
* provide a term bound.
* E.g. PI = 3 F(1/2, 1, 1, 8/5 ; 3/5, 4/3, 5/3 | 2/27) with term limit 2/27.
*/
public static class Hypergeometric extends Spigot
{
private final Rational[] a, b;
protected final Rational z;

public Hypergeometric(Rational[] a, Rational[] b, Rational z) {
this(a,b,z,BigInteger.ONE);
}

public Hypergeometric(Rational[] a, Rational[] b, Rational z, BigInteger mul) {
super(mul);
this.a=a.clone();
this.b=new Rational[b.length+1];
this.b[0]=new Rational(1,1);
System.arraycopy(b,0,this.b,1,b.length);
this.z=z;
}

protected Rational termRatio() {
Rational rv=z;
for (Rational num : a) rv=rv.mul(num.add(i));
for (Rational num : b) rv=rv.div(num.add(i));
return rv;
}
}

public static class Exp extends Hypergeometric
{
public Exp(Rational z) {
super(new Rational[0], new Rational[0], z);
}

// We rely on the fact that the term ratio for Exp always decreases.
public int output(int numdigits) {
int i=0;
return i;
}
}

public static void main(String[] args) {
Exp e=new Exp(Rational.ONE);
int count=e.output(100);
System.out.println(e);
System.out.println("Iterations required: "+(count+1));

Exp e2=new Exp(new Rational(2,1));
count=e2.output(100);
System.out.println(e2);
System.out.println("Iterations required: "+(count+1));

Exp em2=new Exp(new Rational(-2,1));
count = em2.output(100);
System.out.println(em2);
System.out.println("Iterations required: "+(count+1));
}
}
`

Implementing the Rational class is left as an exercise. And yes, this code is an argument for limited operator overloading.

-
The spigot algorithm for $e$ is interesting, but it does not allow one to compute $e^x$ without considerable additional work. –  hardmath Jan 21 '11 at 19:34
Oops, I haven't actually included the implementation of e^x. Most of the additional work is done, though. –  Peter Taylor Jan 21 '11 at 20:08

My webMathematica site hardmath is referring to is indeed down at the moment as I have problems getting webM work with the latest tomcat.

A few comments: - the question of the o/p was about "fastest way". That should be made a bit more precise. Usually people mean fast convergence rates with that. But "fast" in its original sense really means running time. And contractions of continued fractions are a good example of dramatically increasing the convergence rate, but at the expense of getting more and more complicated terms to compute for that. In other words, you need much fewer terms to attain higher precision, but every term is more complex than the simpler ones, of which you need more. Sometimes more cheap steps are faster, and sometimes few expensive terms are faster. It depends on how the computation is done. High-level language? A VM? Machine code? Assember? - hardmath, regarding your third paragraph. Every power series, not just Taylor series, can be converted into an EQUIVALENT continued fraction expansion. The converse is not generally true, especially for delta fractions. The continued fraction expansion is (almost?) always more stable than a power series solution as the terms don't grow as fast. In a power series you are adding many terms of which numerator and denominator are growing at very fast rates, incurring a lot of numerical problems. The c/f expansions are (usually) much more stable, in fact, the parabola theorem states that Stieltjes, C-, and several other c/f expansions are even self-correcting if the backwards recursion is used. And many power series are outright trash if the function to be approximated has poles, which isn't the case here with exp(x), so I won't go into that (just consider that the p/s for tan(x) involves Bernoulli polynomials, converges slowly, converges only in a small interval between poles, whereas the c/f expansion converges rapidly with very simple terms in the whole complex plane EXCEPT at the poles -- which form a null set -- and is very stable with small terms in the convergents). - there's also a few papers on the web that describe how sin, cos, exp, and log are computed by the chip manufacturers. The fastest way to compute these functions is to burn them on a chip and let the semiconductors do it for you, than writing software. The math they use for this is truly mind-boggling. A combination of numerical math, implementation on a chip, and sometimes GPU-based parallelization that is beyond fascinating.

Feel free to contact me through my website, www.lauschkeconsulting.com/contact.html if you need more information, I can provide you with my contraction formulas of the formulas hardmath is referring to (in the Cuyt et al book). The convergence rates are enormous, but the terms also get pretty complicated very fast, which makes numerical computations a bit slower again as more complicated terms have to be evaluated, than if you didn't contract the c/f.

Added on 20110201: My webMathematica site is live again, I solved the tomcat problem.

-