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

 

e^ip+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 = Öe^p.

- 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 a^x (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 e² 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 e^p/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 ( a_k) 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 (e^x-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=2^m, 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 p_k/q_k be the convergent of order k of a simple continued fraction of a number x and let x=[a₀;a₁,a₂,...,a_k,...], 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 a_k=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 10⁴ decimal places, with k=12000, 10⁵ 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 e^t 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 e^t :

é ë  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=(e^t)^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(t²),O(t⁴),O(t⁶),... 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 s_n 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 u_n+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 u_n 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 u_n ³ 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 S₃=(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 S_n 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 d_n be the number of derangements of n different objects (d_n 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 d_n =(n-1)(d_n-1+d_n-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 A_n, the arithmetic mean of the quantities (1,2,3,...,n) and G_n, 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 p_i 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=10⁶, 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 u₁=v₁=0,u₂=v₂=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(n²/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)⁴, and more generally
  

  e=exp(1)=exp( 1/2N) 2N,

  which needs N consecutive squares after the computation of exp(1/2^N). 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(y_n) 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,loge²,loge³ and loge¹⁰, 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 T_EX by T_TH, version 3.01.
On 14 Dec 2001, 15:46.