An efficient way to compute mathematical constant e

The standard representatio开发者_Python百科n of constant e as the sum of the infinite series is very inefficient for computation, because of many division operations. So are there any alternative ways to compute the constant efficiently?

---

Since it's not possible to calculate every digit of 'e', you're going to have to pick a stopping point.

double precision: 16 decimal digits

For practical applications, "the 64-bit double precision floating point value that is as close as possible to the true value of 'e' -- approximately 16 decimal digits" is more than adequate.

As KennyTM said, that value has already been pre-calculated for you in the math library. If you want to calculate it yourself, as Hans Passant pointed out, factorial already grows very fast. The first 22 terms in the series is already overkill for calculating to that precision -- adding further terms from the series won't change the result if it's stored in a 64 bit double-precision floating point variable. I think it will take you longer to blink than for your computer to do 22 divides. So I don't see any reason to optimize this further.

thousands, millions, or billions of decimal digits

As Matthieu M. pointed out, this value has already been calculated, and you can download it from Yee's web site.

If you want to calculate it yourself, that many digits won't fit in a standard double-precision floating-point number. You need a "bignum" library. As always, you can either use one of the many free bignum libraries already available, or re-invent the wheel by building your own yet another bignum library with its own special quirks.

The result -- a long file of digits -- is not terribly useful, but programs to calculate it are sometimes used as benchmarks to test the performance and accuracy of "bignum" library software, and as stress tests to check the stability and cooling capacity of new machine hardware.

One page very briefly describes the algorithms Yee uses to calculate mathematical constants.

The Wikipedia "binary splitting" article goes into much more detail. I think the part you are looking for is the number representation: instead of internally storing all numbers as a long series of digits before and after the decimal point (or a binary point), Yee stores each term and each partial sum as a rational number -- as two integers, each of which is a long series of digits. For example, say one of the worker CPUs was assigned the partial sum,

... 1/4! + 1/5! + 1/6! + ... .

Instead of doing the division first for each term, and then adding, and then returning a single million-digit fixed-point result to the manager CPU:

// extended to a million digits
1/24 + 1/120 + 1/720 => 0.0416666 + 0.0083333 + 0.00138888

that CPU can add all the terms in the series together first with rational arithmetic, and return the rational result to the manager CPU: two integers of perhaps a few hundred digits each:

// faster
1/24 + 1/120 + 1/720 => 1/24 + 840/86400 => 106560/2073600

After thousands of terms have been added together in this way, the manager CPU does the one and only division at the very end to get the decimal digits after the decimal point.

Remember to avoid PrematureOptimization, and always ProfileBeforeOptimizing.

---

If you're using double or float, there is an M_E constant in math.h already.

#define M_E         2.71828182845904523536028747135266250   /* e */

There are other representions of e in http://en.wikipedia.org/wiki/Representations_of_e#As_an_infinite_series; all the them will involve division.

---

I'm not aware of any "faster" computation than the Taylor expansion of the series, i.e.:

e = 1/0! + 1/1! + 1/2! + ...

or

1/e = 1/0! - 1/1! + 1/2! - 1/3! + ...

Considering that these were used by A. Yee, who calculated the first 500 billion digits of e, I guess that there's not much optimising to do (or better, it could be optimised, but nobody yet found a way, AFAIK)

EDIT

A very rough implementation

#include <iostream>
#include <iomanip>

using namespace std;

double gete(int nsteps)
{
  // Let's skip the first two terms
  double res = 2.0;
  double fact = 1;

  for (int i=2; i<nsteps; i++)
  {
    fact *= i;
    res += 1/fact;
  }

  return res;
}

int main()
{
  cout << setprecision(50) << gete(10) << endl;
  cout << setprecision(50) << gete(50) << endl;
}

Outputs

2.71828152557319224769116772222332656383514404296875
2.71828182845904553488480814849026501178741455078125

---

This page has a nice rundown of different calculation methods.

This is a tiny C program from Xavier Gourdon to compute 9000 decimal digits of e on your computer. A program of the same kind exists for π and for some other constants defined by mean of hypergeometric series.

[degolfed version from https://codereview.stackexchange.com/a/33019 ]

#include <stdio.h>
int main() {
      int N = 9009, a[9009], x;
      for (int n = N - 1; n > 0; --n) {
          a[n] = 1;
      }
      a[1] = 2;
      while (N > 9) {
          int n = N--;
          while (--n) {
              a[n] = x % n;
              x = 10 * a[n-1] + x/n;
          }
          printf("%d", x);
      }
      return 0;
  }

This program [when code-golfed] has 117 characters. It can be changed to compute more digits (change the value 9009 to more) and to be faster (change the constant 10 to another power of 10 and the printf command). A not so obvious question is to find the algorithm used.

---

I gave this answer at CodeReviews on the question regarding computing e by its definition via Taylor series (so, other methods were not an option). The cross-post here was suggested in the comments. I've removed my remarks relevant to that other topic; Those interested in further explanations migth want to check the original post.

---

The solution in C (should be easy enough to adapt to adapt to C++):

#include <stdio.h>
#include <math.h>

int main ()
{
    long double n = 0, f = 1;
    int i;
    for (i = 28; i >= 1; i--) {
        f *= i;  // f = 28*27*...*i = 28! / (i-1)!
        n += f;  // n = 28 + 28*27 + ... + 28! / (i-1)!
    }  // n = 28! * (1/0! + 1/1! + ... + 1/28!), f = 28!
    n /= f;
    printf("%.64llf\n", n);
    printf("%.64llf\n", expl(1));
    printf("%llg\n", n - expl(1));
    printf("%d\n", n == expl(1));
}

Output:

2.7182818284590452354281681079939403389289509505033493041992187500
2.7182818284590452354281681079939403389289509505033493041992187500
0
1

There are two important points:

1. This code doesn't compute 1, 1*2, 1*2*3,... which is O(n^2), but computes 1*2*3*... in one pass (which is O(n)).

2. It starts from smaller numbers. If we tried to compute

   1/1 + 1/2 + 1/6 + ... + 1/20!

   and tried to add it 1/21!, we'd be adding

   1/21! = 1/51090942171709440000 = 2E-20,

   to 2.something, which has no effect on the result (double holds about 16 significant digits). This effect is called underflow.

   However, when we start with these numbers, i.e., if we compute 1/32!+1/31!+... they all have some impact.

This solution seems in accordance to what C computes with its expl function, on my 64bit machine, compiled with gcc 4.7.2 20120921.

---

You may be able to gain some efficiency. Since each term involves the next factorial, some efficiency may be obtained by remembering the last value of the factorial.

e = 1 + 1/1! + 1/2! + 1/3! ...  

Expanding the equation:

e = 1 + 1/(1 * 1) + 1/(1 * 1 * 2) + 1/(1 * 2 * 3) ...

Instead of computing each factorial, the denominator is multiplied by the next increment. So keeping the denominator as a variable and multiplying it will produce some optimization.

---

If you're ok with an approximation up to seven digits, use

3-sqrt(5/63)
2.7182819

If you want the exact value:

e = (-1)^(1/(j*pi))

where j is the imaginary unit and pi the well-known mathematical constant (Euler's Identity)

---

There are several "spigot" algorithms which compute digits sequentially in an unbounded manner. This is useful because you can simply calculate the "next" digit through a constant number of basic arithmetic operations, without defining beforehand how many digits you wish to produce.

These apply a series of successive transformations such that the next digit comes to the 1's place, so that they are not affected by float rounding errors. The efficiency is high because these transformations can be formulated as matrix multiplications, which reduce to integer addition and multiplication.

In short, the taylor series expansion

e = 1/0! + 1/1! + 1/2! + 1/3! ... + 1/n!

Can be rewritten by factoring out fractional parts of the factorials (note that to make the series regular we've moved 1 to the left side):

(e - 1) = 1 + (1/2)*(1 + (1/3)*(1 + (1/4)...))

We can define a series of functions f1(x) ... fn(x) thus:

f1(x) = 1 + (1/2)x
f2(x) = 1 + (1/3)x
f3(x) = 1 + (1/4)x
...

The value of e is found from the composition of all of these functions:

(e-1) = f1(f2(f3(...fn(x))))

We can observe that the value of x in each function is determined by the next function, and that each of these values is bounded on the range [1,2] - that is, for any of these functions, the value of x will be 1 <= x <= 2

Since this is the case, we can set a lower and upper bound for e by using the values 1 and 2 for x respectively:

lower(e-1) = f1(1) = 1 + (1/2)*1 = 3/2 = 1.5
upper(e-1) = f1(2) = 1 + (1/2)*2 = 2

We can increase precision by composing the functions defined above, and when a digit matches in the lower and upper bound, we know that our computed value of e is precise to that digit:

lower(e-1) = f1(f2(f3(1))) = 1 + (1/2)*(1 + (1/3)*(1 + (1/4)*1)) = 41/24 = 1.708333
upper(e-1) = f1(f2(f3(2))) = 1 + (1/2)*(1 + (1/3)*(1 + (1/4)*2)) = 7/4 = 1.75

Since the 1s and 10ths digits match, we can say that an approximation of (e-1) with precision of 10ths is 1.7. When the first digit matches between the upper and lower bounds, we subtract it off and then multiply by 10 - this way the digit in question is always in the 1's place where floating-point precision is high.

The real optimization comes from the technique in linear algebra of describing a linear function as a transformation matrix. Composing functions maps to matrix multiplication, so all of those nested functions can be reduced to simple integer multiplication and addition. The procedure of subtracting the digit and multiplying by 10 also constitutes a linear transformation, and therefore can also be accomplished by matrix multiplication.

Another explanation of the method: http://www.hulver.com/scoop/story/2004/7/22/153549/352

The paper that describes the algorithm: http://www.cs.ox.ac.uk/people/jeremy.gibbons/publications/spigot.pdf

A quick intro to performing linear transformations via matrix arithmetic: https://people.math.gatech.edu/~cain/notes/cal6.pdf

NB this algorithm makes use of Mobius Transformations which are a type of linear transformation described briefly in the Gibbons paper.

---

From my point of view, the most efficient way to compute e up to a desired precision is to use the following representation:

e := lim (n -> inf): (1 + (1/n))^n

Especially if you choose n = 2^x, you can compute the potency with just x multiplications, since:

a^n = (a^2)^(n/2), if n % 2 = 0

---

The binary splitting method lends itself nicely to a template metaprogram which produces a type which represents a rational corresponding to an approximation of e. 13 iterations seems to be the maximum - any higher will produce a "integral constant overflow" error.

#include <iostream>
#include <iomanip>

template<int NUMER = 0, int DENOM = 1>
struct Rational
{
    enum {NUMERATOR = NUMER};
    enum {DENOMINATOR = DENOM};

    static double value;
};

template<int NUMER, int DENOM>
double Rational<NUMER, DENOM>::value = static_cast<double> (NUMER) / DENOM;

template<int ITERS, class APPROX = Rational<2, 1>, int I = 2>
struct CalcE
{
    typedef Rational<APPROX::NUMERATOR * I + 1, APPROX::DENOMINATOR * I> NewApprox;
    typedef typename CalcE<ITERS, NewApprox, I + 1>::Result Result;
};

template<int ITERS, class APPROX>
struct CalcE<ITERS, APPROX, ITERS>
{
    typedef APPROX Result;
};

int test (int argc, char* argv[])
{
    std::cout << std::setprecision (9);

    // ExpType is the type containing our approximation to e.
    typedef CalcE<13>::Result ExpType;

    // Call result() to produce the double value.
    std::cout << "e ~ " << ExpType::value << std::endl;

    return 0;
}

Another (non-metaprogram) templated variation will, at compile-time, calculate a double approximating e. This one doesn't have the limit on the number of iterations.

#include <iostream>
#include <iomanip>

template<int ITERS, long long NUMERATOR = 2, long long DENOMINATOR = 1, int I = 2>
struct CalcE
{
    static double result ()
    {
        return CalcE<ITERS, NUMERATOR * I + 1, DENOMINATOR * I, I + 1>::result ();
    }
};

template<int ITERS, long long NUMERATOR, long long DENOMINATOR>
struct CalcE<ITERS, NUMERATOR, DENOMINATOR, ITERS>
{
    static double result ()
    {
        return (double)NUMERATOR / DENOMINATOR;
    }
};

int main (int argc, char* argv[])
{
    std::cout << std::setprecision (16);

    std::cout << "e ~ " <<  CalcE<16>::result () << std::endl;

    return 0;
}

In an optimised build the expression CalcE<16>::result () will be replaced by the actual double value.

Both are arguably quite efficient since they calculate e at compile time :-)

---

@nico Re:

..."faster" computation than the Taylor expansion of the series, i.e.:

e = 1/0! + 1/1! + 1/2! + ...

or

1/e = 1/0! - 1/1! + 1/2! - 1/3! + ...

Here are ways to algebraically improve the convergence of Newton’s method:

https://www.researchgate.net/publication/52005980_Improving_the_Convergence_of_Newton's_Series_Approximation_for_e

It appears to be an open question as to whether they can be used in conjunction with binary splitting to computationally speed things up. Nonetheless, here is an example from Damian Conway using Perl that illustrates the improvement in direct computational efficiency for this new approach. It’s in the section titled “