Mar 11, 2009

The constant $e$ and its computation

http://numbers.computation.free.fr/Constants/E/e.html

The constant e and its computation

e = 2.71828182845904523536028747135266249775724709369995...

(Click here for a Postscript version of this page.)

 

eip+1=0 , all analysis lies here

- Felix Klein (1849-1925).

Gentlemen, we have not the slightest idea what this equation means, but we may be sure that it means something very important (to his students about the formula i-i = Öep.

- Benjamin Peirce (1809-1880).

1  Introduction

1.1  A function equal to it's own derivative

Let a be a real positive number, the exponential function ax (in base a) is a differentiable function, that is the following limit exists

(ax)¢=
lim
h®0 

æ
è
 ax+h-ax

h
ö
ø
=
lim
h®0 

æ
è
 ah-1

h
ö
ø
ax=Cax.

A very important case is given when the derivative of the exponential function is equal to itself, which implies

C=1=
lim
h®0 

æ
è
 ah-1

h
ö
ø
,
and this may also be written as
a=
lim
h®0 
( 1+h) 1/h.

This limit is well defined and it was denoted by the letter e by the Swiss mathematician Leonhard Euler (1707-1783), first around the end of year 1727 in a manuscript entitled Meditatio in Experimenta explosione tormentorum nuper instituta (Meditation upon experiments made recently on the firing of Canon, [21]), then in a letter to a friend Goldbach (1690-1764) in 1731 and later in 1748 in his work [5]. The origin of this choice is maybe due to the first letter of the word exponential or just because it was the first free vowel (a was already used by Euler !). It should be pointed out that Leibniz (1646-1716) was the first, in 1690, to recognize e as a constant and he used the notation b (more about the early origins of e may be found in [4]).

Hence the constant e is defined by the monotone increasing sequence :

Definition 1

e=
lim
n®¥ 
en=
lim
n®¥ 

æ
è
1+  1

n
ö
ø
n

 
.
(1)

But the convergence is extremely slow and not convenient to make an accurate determination of the constant :


ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
e1=2.(000...),
e10=2.(593...),
e100=2.7(048...),
e1000=2.71(692...),
e10000=2.718(145...),
...

A quicker convergence is easily obtained by a modified version of this definition (see [15], [3], [11]) :

e=
lim
n®¥ 

æ
è
 2n+1

2n-1
ö
ø
n

 
.
(2)

It's interesting to note that the constant e holds a key position in the simplest first order differential equation

y¢=y,
while p occurs in the similar second order differential equation.

Links between the definition (1) and the law of compound interest in financial calculations are given in [13] which is a nice book completely dedicated to the constant e and its mathematical background. An on-line essay is also available at [6].

1.2  A famous sequence

Using the Newton's binomial theorem

(1+x)n=1+nx+  n(n-1)

2
x2+...+xn,
for x=1/n and after some careful manipulations gives the famous series expansion
e=1+  1

1
+  1

1.2
+  1

1.2.3
+  1

1.2.3.4
+...=
lim
n®¥ 
sn=
lim
n®¥ 

æ
è

å
k=0 
n  1

k!
ö
ø
,
(3)
given by Euler in 1748 [5] and with x=-1/n it also produces the fast converging and alternating series
e-1=1-  1

1
+  1

1.2
-  1

1.2.3
+  1

1.2.3.4
-...
(4)

Relation (3) is very efficient to compute e because the factorial of a number grows very quickly and it's immediate to show that

sn < e < sn+  1

n.n!
.
Euler used this relation to find 23 digits of e [5].

The constant e is also known as the natural base of logarithm, that is the number whose natural logarithm is 1 :

log(e)= ó
õ
e

1 

 dx

x
=1,
and it is equal to the value of the exponential function exp(z) at 1, where :
exp(z)=1+  z

1!
+  z2

2!
+  z3

3!
+...+  zn

n!
+...
(5)
so
e=exp(1).

1.3  Nature of the constant e

Theorem 2 (Euler 1737). e and e2 are irrational numbers.

Proof. A very easy proof of the irrationality of e can be found in Irrationality proofs.  

And some years later, Johann Heinrich Lambert (1728-1777) established the stronger result :

Theorem 3 (Lambert 1768). Let p/q being a nonzero rational number then ep/q is an irrational number.

Proof. A proof is also detailed in Irrationality proofs.  

But it took more than a century to prove the transcendence of e. Before this major milestone, the stronger result was due to Joseph Liouville (1809-1882) who proved in 1844 that e could not be the root of any quadratic equation with integral coefficients.

Theorem 4 (Hermite 1873, [9]). e is a transcendental number.

Proof. An elegant proof due to Hilbert may be find in [1].  

Here is what Charles Hermite (1822-1901) claimed after this successful proof of transcendence : ''I shall risk nothing on an attempt to prove the transcendence of p. If others undertake this enterprise, no one will be happier than I in their success. But believe me, it will not fail to cost them some effort''.  Ferdinand von Lindemann (1852-1939) succeeded in this hard task only a few years later in 1882.

2  e and continued fraction

2.1  Euler's results

In 1737 and 1748 [5], Euler published three important developments related to the constant e. It is more convenient to use for the continued fractions

x=a0+  1

a1+  1

a2+...

the following more compact notation :
x=[a0;a1,a2,a3,...],
where the ( ak) are the partial quotients [8].

The first of those developments is

e=[2;1,2,1,1,4,1,1,6,1,1,8,1,1,10,1,1,12,1,1,14,...]
leading to the fractions
e= æ
è
2,3,  8

3
,  11

4
,  19

7
,  87

32
,  106

39
,  193

71
,  1264

465
,  1457

536
,  2721

1001
,  23225

8544
,  25946

9545
,¼ ö
ø
.

The second one is

Öe=[1;1,1,1,5,1,1,9,1,1,13,1,1,17,1,1,21,1,1,25,...],
the last one has, as Euler noticed, a very regular and elementary pattern

 e-1

2
=[0;1,6,10,14,18,22,26,30,34,38,42,46,50,...]
(6)
and it also produces the following set of rational approximations for e
e= æ
è
0,3,  19

7
,  193

71
,  2721

1001
,  49171

18089
,  1084483

398959
,  28245729

10391023
,  848456353

312129649
,¼ ö
ø

where the last of those fractions is already a very good approximation

 848456353

312129649
=2.71828182845904523(475...).

Notice that, unlike p, the developments involving e and its powers are dramatically regular.

2.2  Generalization

There is an interesting result about the continued fractions of expressions of the form :


 e±1/p-1

2
.

Theorem 5 Let p be a positive integer and p ³ 1:


ì
ï
ï
í
ï
ï
î

 1

2
( e1/p-1) = [0;2p-1,6p,10p,14p,18p,22p,...]

 1

2
( 1-e-1/p) = [0;2p+1,6p,10p,14p,18p,22p,...]

Proof. Left as exercise. Hint : Use the continued and generalized fraction development of the tanh(x) function at x=1/p and conclude with the relation (ex-1)/2=( 1/tanh(x/2)-1) -1.  

Relation (6) is the special case p=1 of this theorem. Such an easy pattern may be convenient to use for values of the form p=2m, and for example with m=5,p=32 it will produce the fast converging development :


 e1/32-1

2
=[0;63,192,320,448,576,704,832,960,1088,...]

and e will then be obtained by computing m=5 consecutive squares.

2.3  An algorithm using continued fraction

Let pk/qk be the convergent of order k of a simple continued fraction of a number x and let x=[a0;a1,a2,...,ak,...], we know that the convergents may be computed at any order thanks to the following recurrence relations


ì
í
î
pk=akpk-1+pk-2
qk=akqk-1+qk-2

with the initial conditions

ì
í
î
p-1=1,p0=a0
q-1=0,q0=1

It's a well known properties on continued fraction that



lim
k®¥ 

æ
è
 pk

qk
ö
ø
=x
with the bound

ê
ê
 pk

qk
-x ê
ê
<  1

qk2
.

If we use the continued fraction (6) for (e-1)/2, then

a0=0,a1=1,
and ak=4(k-1)+2 for k > 1. It's therefore very easy to build an efficient algorithm to compute a fast converging sequence to the constant e. With k=1500, e may be computed to 104 decimal places, with k=12000, 105 digits of e will be reached.

3  Using Padé approximant

If we consider the series expansion

s(t)= ¥
å
k=0 
sktk,sk Î R,
we are looking for a rational approximation of this sequence with numerator of degree m and denominator of degree n such as :
s(t)- æ
ç
ç
è
m
å
k=0 
mktk

n
å
k=0 
nktk
ö
÷
÷
ø
=O(tm+n+1),       t®0.

When such an approximation exists, it's called the Padé approximant of degree (m,n) of the series s. There are many properties associated with those approximants and they are often used to compute series with bad convergence properties. Here we will just make use of the Padé approximants of the et function which have an explicit representation :


m
å
k=0 

 (m+n-k)!m!

(m+n)!k!(m-k)!
tk

n
å
k=0 

 (m+n-k)!n!

(m+n)!k!(n-k)!
(-t)k
.
(7)

To illustrate relation (7), here are the first approximants with n=m and of increasing degree for the function et :


é
ë
 2+t

2-t
,  12+6t+t2

12-6t+t2
,  120+60t+12t2+t3

120-60t+12t2-t3
,  1680+840t+180t2+20t3+t4

1680-840t+180t2-20t3+t4
,... ù
û
.

From this we may deduce that for small values of t, e=(et)1/twill be well approximated by this set of functions


é
ë
æ
è
 2+t

2-t
ö
ø
1/t

 
, æ
è
 12+6t+t2

12-6t+t2
ö
ø
1/t

 
, æ
è
 120+60t+12t2+t3

120-60t+12t2-t3
ö
ø
1/t

 
,... ù
û
,
(8)
the error being respectively O(t2),O(t4),O(t6),... while the famous formula (1+t)1/t converges with an error O(t). Notice that relation (2) is equivalent to the first approximation of the list (8) with t=1/n.

To illustrate the efficiency of those new expressions, the last formula of the list (8) for different values of t=1/n produces


ê
ê
ê
ê
ê
ê
ê
ê
ê
n=20,e » 2.718(309...),
n=21,e » 2.71828(225...),
n=22,e » 2.7182818(350...),
n=25,e » 2.7182818284590(703...),
n=28,e » 2.718281828459045235(456...).

4  Infinite products

There are several formulae to compute e as infinite products, one of them was given in 1873 by Catalan :

e=  2

1

æ
è
 4

3
ö
ø
1/2

 

æ
è
 6.8

5.7
ö
ø
1/4

 

æ
è
 10.12.14.16

9.11.13.15
ö
ø
1/8

 
¼,
and another similar relation is [1] :
e=2 æ
è
 2

1
ö
ø
1/2

 

æ
è
 2.4

3.3
ö
ø
1/4

 

æ
è
 4.6.6.8

5.5.7.7
ö
ø
1/8

 
¼

Like Wallis formula for p, the convergence is extremely slow.

Another product may be easily deduced from formula (3); let sn be the partial sum of this sequence

sn=1+  1

1!
+...+  1

n!

then
sn+1=sn+  1

(n+1)!
=sn æ
è
1+  1

(n+1)!sn
ö
ø
=sn æ
è
1+  1

un+1
ö
ø

clearly un+1 is an integer and

un=n!sn-1
un+1=(n+1)!sn=(n+1)n! æ
è
sn-1+  1

n!
ö
ø
=(n+1)(n!sn-1+1)

therefore, if we define the sequence un as following

u1=1
un+1=(n+1)(un+1)

we find
e= ¥
Õ
n=1 

æ
è
 un+1

un
ö
ø
= æ
è
 2

1
ö
ø

æ
è
 5

4
ö
ø

æ
è
 16

15
ö
ø

æ
è
 65

64
ö
ø

æ
è
 326

325
ö
ø
...,
with obviously un ³ n!.

The rate of convergence of this product is exactly the same as (3).

5  e and permutations

The number e occurs in the derangements of n objects. For example consider the set of 3 different objects S3=(1,2,3), all the permutations of this set are

(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1)
and it is well known that the number of permutations of Sn is n! (6 in our example). But among some permutations how many don't have a fix point ? In the example the answer is 2, that is
(2,3,1),(3,1,2)
but the general answer for a set of n objects is given by the following theorem.

Theorem 6 (Euler) Let dn be the number of derangements of n different objects (dn is sometime called subfactorial), we have

dn =n! æ
è
1-  1

1!
+  1

2!
-  1

3!
+...+  (-1)n

n!
ö
ø
.

Proof. Left as exercise. Hint: Show the recursion dn =(n-1)(dn-1+dn-2).  

Hence the ratio of the number of derangements to the number of permutations is


 dn

n!
= æ
è
1-  1

1!
+  1

2!
-  1

3!
+...+  (-1)n

n!
ö
ø
,
this ratio is such as


lim
n®¥ 

æ
è
 dn

n!
ö
ø
=  1

e
=0.3678794411714423215955237...
and for example

ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê

 d3

3!
=  2

3!
=0.3333333333...

 d4

4!
=  9

4!
=0.375

 d5

5!
=  44

5!
=0.3666666666...

 d6

6!
=  265

6!
=0.3680555555...

 d7

7!
=  1854

7!
=0.3678571428...

6  Other results

6.1  e and Means

Before the description of efficient ways to compute e, let's consider An, the arithmetic mean of the quantities (1,2,3,...,n) and Gn, the geometric mean of the same quantities, that is

An=  1+2+3+...+n

n
=  n(n+1)

2n
,
and


log(Gn)=  log(1)+log(2)+...+log(n)

n

or
Gn=( 1.2.3...n) 1/n=(n!)1/n.

Thanks to Stirling's formula, when n becomes large :


 An

Gn
=  n+1

2(n!)1/n
~  n+1

2( nne-n(2pn)1/2) 1/n
=  e

2

 n+1

n(2pn)1/(2n)

therefore


lim
n®¥ 

 An

Gn
=  e

2
.

A more general result was stated in [10] and another similar result involving the geometric mean of the prime numbers was given in [12] :



lim
n®¥ 

æ
è

Õ
pi £ n 
pi ö
ø
1/n
 
=e
where the pi is the sequence of primes (2,3,5,7,11,13,...). The rate of converge is extremely slow and for example by computing all primes up to n=106, it only gives 2.7141645...

6.2  e and p in a mirror

As observed by B. Cloitre (from who we borrow the title of this section), the two following symmetrical algorithms converge respectively to the constants e and p. Let u1=v1=0,u2=v2=1 and define the iteration :


ì
ï
ï
í
ï
ï
î
un+2=un+1+  1

n
un
vn+2=  1

n
vn+1+vn
,
then it's easy to show that :


lim
n®¥ 

æ
è
 n

un
ö
ø
=e
and


lim
n®¥ 

æ
è
 2n

vn2
ö
ø
=p.

The first algorithm is just a reformulation of (4) and therefore converges quickly while the second one is derived from Wallis' formula for p and is very slow to converge (it can be a little accelerated by replacing 2n by 2n-1).

7  Direct use of the exponential series

Formula (3) is a good one to compute e while formula (4) may be used to check the result. To have d decimal digits of e, roughly the first K @ dlog(10)/log(d) terms of the series are sufficient.

  • Classical approach. Using the series (3), the classical approach has a complexity of O(n2/log(n)).  The algorithm is simply of the form :

    ì
    ï
    í
    ï
    î
    v=v/n
    u=u+v
    n=n+1
    ,
    (9)

    starting with u=v=n=1 and (9) is repeated until v is smaller than the required precision. A clear and easy C source code which computes e with the classical approach can be found here : Easy programs for constants computation. Note that it may be interesting to reverse and rearrange series (3) in a Horner's way

    e= æ
    è
    æ
    è
    æ
    è
    æ
    è
    ( ...+1)  1

    6
    +1 ö
    ø

     1

    5
    +1 ö
    ø

     1

    4
    +1 ö
    ø

     1

    3
    +1 ö
    ø

     1

    2
    +2

    so that, if we evaluate it by starting with the last terms of series (3), only divisions by small integers will be required (the addition of 1 being of a quasi-null cost when working with a large accuracy) and fewer storage memory is needed. The algorithm, which converges to e-1, now looks like :


    ì
    í
    î
    u=(1+u)/n
    n=n-1

    (10)

    with initial values u=1, n=K and (10) should be repeated until n=1.

  • Binary splitting approach. The series (3) is particularly well suited for the binary splitting method ( the binary splitting page contains a well detailed section dedicated to the computation of the exponential series). It leads to the best results in the practice. Moreover, it is quite easy to implement if one has already a fast multiplication code. A quite clear and short C source code which computes e with the binary splitting approach can be found in Easy programs for constants computation.
  • Powering. For all number x the exponential series (5) is converging. Formula (3) corresponds to the particular case x=1.

    Instead of computing e=exp(1) from its series, one can compute exp(1/2) and square the result. The advantage is that the series of exp(1/2) has a quicker convergence than the one of exp(1). The process can also be used with exp(1)=exp(1/4)4, and more generally

    e=exp(1)=exp( 1/2N) 2N,
    which needs N consecutive squares after the computation of exp(1/2N). This technique can be interesting if one has a fast squaring procedure compared to the series computation. The choice of N can be optimized with respect to the relative timings of series computation and squaring. However, it seems that if one use the binary splitting approach on (3), the powering optimization does not permits to obtain better global timings.

8  AGM based approach

The log(2) page presents an AGM based iteration which converges quadratically to log(x) for any real number x. Thus the computation of y=exp(a) can be achieved by using a simple Newton iteration on the function f(y)=log(y)-a, which leads to the iteration

yn+1=yn( 1+a-log(yn)) .

The convergence of the process is quadratic, that is the number of correct digits is doubled at each iteration. Notice that each iteration requires the computation of log(yn) to the current precision required, which can be done using the AGM based algorithm to compute logarithmic values.

This process permits to compute exp(a) for any real number a, which for the (very) particular case a = 1 gives the value of e.

It must be pointed out that using the AGM based process to compute e is quite complicated and less efficient than the binary splitting approach from the exponential series.

9  A tiny program to compute e

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 p and for some other constants defined by mean of hypergeometric series.

main(){int N=9009,n=N,a[9009],x;while(--n)a[n]=1+1/n; for(;N>9;printf("%d",x)) for(n=N--;--n;a[n]=x%n,x=10*a[n-1]+x/n);}  

This program has 117 characters (try to do better !). 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.

10  Records of computation

The following table illustrates the evolution of the number D of computed digits for the constant e :

D When Who Notes
23 1748 L. Euler (3) was used [5]
137 1853 W. Shanks [19]
205 1871 W. Shanks [20]
346 1884 Boorman [2]
808 1946 ?
2 010 1949 John von Neumann He and his team used the ENIAC [14]
100 265 1961 D. Shanks & J.W. Wrench 2.5 hours on a IBM 7090. (3) was used [17]
10 000 000 1994 R. Nemiroff & J. Bonnell
18 199 978 1997 May Patrick Demichel
20 000 000 1997 Aug Birger Seifert 182.5 hours on a Pentium (133Mhz)
50 000 817 1997 Sep P. Demichel 714 hours on a HP 9000/778 (160Mhz)
200 000 579 1999 S. Wedeniwski Binary splitting
869 894 101 1999 Oct S. Wedeniwski Continued fraction expansion of e
1 250 000 000 1999 Nov X. Gourdon (3) was used, the verification was made with (4). A binary splitting technique was used. The computation took  40 hours on a PentiumII 350 with 320 Mo. 4 Go of disk space was used during the process. (Note : a previous 1.7 billion digits computation by Patrick Demichel was made before, but no verification was performed. It appears that its first 1.25 billion digits are OK).
2 147 483 648 2000 Jul X. Gourdon & S. Kondo The computation was launched by Shigeru Kondo with the program PiFast33 written by X. Gourdon. The same technique as in the previous record was used. The computation took 21 hours on a PentiumIII 933 with 512 Mo. 8 Go of disk space were used during the process.
3 221 225 472 2000 Jul X. Gourdon & C. Martin PiFast33 on an Athlon 650 with just 128 Mo of physical memory and 17.2 GB of disk space. The initial calculation took just over 77 hours and completed on July 13th and the verification took 80 hours.
6 442 450 944 2000 Aug X. Gourdon & S. Kondo PiFast33, the computation took 70 hours on a PentiumIII 800 with 1024 Mo. 24 Go of disk space were necessary. The verification took 71 hours on the same machine
12 884 901 000 2000 Aug X. Gourdon & S. Kondo PiFast33, the computation, ended on Aug 04 2000, was done in 167 hours on a Pentium III 800 with 1024 Mo of physical memory. 48 GB of disk space was needed.

References

[1]
J.M. Borwein and P.B. Borwein, Pi and the AGM - A study in Analytic Number Theory and Computational Complexity, A Wiley-Interscience Publication, New York, (1987)

[2]
Boorman, Computation of the Naperian Base, The Mathematical Magazine, (1884), vol. 1, p. 204

[3]
Brothers, Harlan J. and J.A. Knox, New closed-form approximations to the logarithmic constant e, Math. Intelligencer, (1998)

[4]
J.L. Coolidge, The number e, Amer. Math. Monthly, (1950), vol. 57, p. 591-602

[5]
L. Euler, Introduction à l'analyse infinitésimale (french traduction by Labey), Barrois, ainé, Librairie, (original 1748, traduction 1796), vol. 1, p. 89-90

[6]
S. Finch, Favorite mathematical constants. World Wide Web site at http://pauillac.inria.fr/algo/bsolve/constant/constant.html, (1995)

[7]
X. Gourdon and P. Sebah, Numbers, Constants and Computation, World Wide Web site at the adress : http://numbers.computation.free.fr/Constants/constants.html, (1999)

[8]
G.H. Hardy and E.M. Wright, An Introduction to the Theory of Numbers, Oxford Science Publications, (1979)

[9]
C. Hermite, Sur la fonction exponentielle, C. R. Académie des Sciences, (1873), vol. 77, p. 18-24, 74-79, 226-233, 285-293

[10]
K. Knopp, Theory and application of infinite series, Blackie & Son, London, (1951)

[11]
J.A. Knox and Harlan J. Brothers, Novel series-based approximations to e, College Math. J., (1998)

[12]
F. Le Lionnais, Les nombres remarquables, Paris, Hermann, (1983), p. 47

[13]
E. Maor, e: The Story of a Number, Princeton University Press, (1994)

[14]
G.W. Reitwiesner , An ENIAC Determination of p and e to more than 2000 Decimal Places, Mathematical Tables and other Aids to Computation, (1950), vol. 4, p. 11-15

[15]
L.W. Richardson, The deferred Approach to the Limit, Philosophical Transactions of the Royal Society of London, (1927), serie A, vol. 226

[16]
A. Sale, The Calculation of e to many Significant Digits, Computing Journal, (1968), vol. 11, p. 229-230

[17]
D. Shanks and J.W. Wrench, Jr., Calculation of p to 100,000 Decimals, Math. Comput., (1962), vol. 16, p. 76-79

[18]
D. Shanks and J.W. Wrench, Jr., Calculation of e to 100,000 Decimals, Math. Comput., (1969), vol. 23, p. 679-680

[19]
W. Shanks, Contributions to Mathematics Comprising Chiefly the Rectification of the Circle to 607 Places of Decimals, G. Bell, London, (1853), p. vii.

[20]
W. Shanks, Second paper on the numerical values of e,loge2,loge3 and loge10, also on the numerical value of M the modulus of the common system of logarithms, all to 205 decimals, Proc. of London, (1871), vol. 19, p. 27-29

[21]
D.E. Smith, A Source Book in Mathematics, Dover Publications, New York, (1959, first edition 1929)

[22]
E.W. Weisstein, CRC Concise Encyclopedia of Mathematics, CRC Press, (1999), p. 503-504

 


 

Back to Numbers, Constants and Computation


File translated from TEX by TTH, version 3.01.
On 14 Dec 2001, 15:46.

No comments: