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

Question

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

Answer 1

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.

Added:

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)]
http://books.google.com/books?id=DQtpJaEs4NIC&dq=exponential+function+continued+fraction+convergence

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.

Answer 2

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.

Answer 3

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?

Answer 4

Fastest way is to use Newton method of calculation of square root $x_{n+1}=0.5 (x_n + a/x_n)$ from $\mathsf e$, it will give you $$\exp(^2.1),$$ where $^2$ denotes a binary number.

Then you calculate square root of square root, giving you $$\exp(^2.01)$$

And so on. This will allow you to create tables of logarithms and exponents. This is the fastest way to do by hand calculation and it was used historically.

Pade approximations are good, but they do not converge as fast as Newton method, which effectively doubles your "accurate number count" by single iteration.

Then you multiply the required amount of "roots" you have, effectively stacking any binary number you need to exponentiate, and you also is able to reuse your roots to build an "acceleration table", where you keep exponents for often used numbers like $\exp(^2.101)$

Note that you can shift between different exponentiation digits by "rooting" of your number, so if you take square root you obtain $$\sqrt{\exp(^2.101)}=\exp(^2.0101)$$

It is a very strong cost-efficient method, but it requires from calculator to keep square roots tables and "often used" exponents.

Especially if you are allowed to have large enough exponentiation table with separate roots of high precision, you can achieve super strong results with this. Modern computers don't do it because they are using some Taylor series or variation of Taylor series, because they are "saving" memory for each call. Also they can be using generalized continued fractions like this https://en.wikipedia.org/wiki/Exponential_function#Continued_fractions_for_ex, but in my expirience GCFs more often a compromise between "table" approach and "series". Something in between. Not super speed.

PS. You can further accelerate it by good searching algo which will determine patterns in your binary number and re-use table ot exponents for those patterns. Typical patterns are $^211, ^2101, ^2111, ^21111, ^211011$ and so on. Also you can further accelerate it by using stronger approximations for starting iteration for your Newton root-finder, something based on GCFs like that $$\sqrt{z} = \sqrt{x^2+y} = x+\cfrac{2x \cdot y} {2(2z - y)-y-\cfrac{y^2} {2(2z - y)-\cfrac{y^2} {2(2z - y)-\ddots}}}$$

This will be super efficient computation-wise if you prepare some tables of squares and re-use those every time for your initial GCF guess.

PPS. Note that binary search algo could be enhanced to not only look for additive patterns, but also for negating patterns. You can think of binary sequence like that $^2.1011101=^21-.0100010$ minimizing required roots amount from 5 to 2.

Answer 5

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.
     */
    public int advance(Rational beta) {
        if (verify && beta.signum()0) throw new IllegalArgumentException(Et_t+" exceeds bound "+beta);
        h00=h00.mul(Et_t);
        i=i.add(BigInteger.ONE);

        // Produce (if we can: without a converging bound we can't).
        int rv=0;
        if (beta.compareTo(Rational.ONE)0; i++) numdigits-=advance(beta);
        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;
            for (; numdigits>0; i++) numdigits-=advance(termRatio().abs());
            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.

Answer 6

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.