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
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 :
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
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
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 :
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
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
the following more compact notation :
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 :
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
with the initial conditions
It's a well known properties on continued fraction that
with the bound
If we use the continued fraction (6) for (e-1)/2, then
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)- |
æ ç
ç è
|
|
ö ÷
÷ ø
|
=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=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
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+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
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
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... |
|
|
|
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
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 :
then it's easy to show that :
and
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 :
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 :
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
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
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.