Вычисление с большой точностью числа е

Путь: Математика » Быстрые вычисления » Число е - мощные методы

Вычисление с большой точностью числа е

Если Вы найдете статью полезной и интересной - не сочтите за труд, переведите материал или хотя бы часть его и отправьте на адрес algolist@manual.ru. Задать вопросы или просто написать письмо можно также на странице контактов.

Статья любезно предоставлена:
© Xavier Gordon, numbers.computation.free.fr
The constant e

e = 2.71828182845904523536028747135266249775724709369996ј

" e^ip+1 = 0 , all analysis lies here"

Felix Klein (1849-1925).

1 Introduction

A function equal to it's own derivative

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

(a^x)^ў =

lim
h® 0

ж
з
и

a^x+h-a^x

ц
ч
ш

lim
h® 0

ж
з
и

a^h-1

ц
ч
ш

a^x = Ca^x.

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

C = 1 =

lim
h® 0

ж
з
и

a^h-1

ц
ч
ш

and can be written as

a =

lim
h® 0

( 1+h) ^1/h.

This limit exists and it was denoted by the letter e by Euler, first in a letter to Goldbach in 1731, later in 1736 in his work. The origin of this choice is maybe due to the first letter of the word exponential. Hence the constant e is defined by the monotone increasing sequence

e =

lim
n® Ґ

e_n =

lim
n® Ґ

ж
з
и

ц
ч
ш

but the convergence is very slow

e₁

e₁₀

2.(59374246010...),

e₁₀₀

2.7(0481382942...),

e₁₀₀₀

2.71(692393223...),

e₁₀₀₀₀

2.718(14592682...),...

A quicker convergence is easily obtained by (see [5], [12], [13])

e =

lim
n® Ґ

ж
з
и

2n+1

2n-1

ц
ч
ш

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 a similar second order differential equation.

Links between the previous definition and the law of compound interest in financial calculations are given in [11].

A famous sequence

Using the Newton's binomial theorem

(1+x)ⁿ = 1+nx+

n(n-1)

x²+...+xⁿ,

for x = 1/n and after some manipulations gives the famous relation

1×2

1×2×3

1×2×3×4

+ј

(1)

lim
n® Ґ

s_n =

lim
n® Ґ

ж
з
и

n
е
k = 0

ц
ч
ш

(2)

given by Euler in 1748 [1]. This relation is very efficient to compute e because the factorial of a number grows very quickly and it's easy to show that

s_n < e < s_n+

n.n!

Euler used this relation to find 23 digits of e.

The constant e is also known as the natural base of logarithm, that is

log(e) =

у
х

e

1

= 1,

and is equal to the value of the exponential function at 1 :

e = exp(1).

The number e is irrational (Euler 1744) and transcendental (Hermite 1873, [3]).

2 Continued fraction

In 1737 and 1748 [1], Euler published three important developments related to the constant e. It is more convenient to use the following notation:

x = a₀+

a₁+

a₂+...

= [a₀;a₁,a₂,a₃,...]

The first development is

[2;1,2,1,1,4,1,1,6,1,1,8,1,1,10,1,1,12,...]

ж
з
и

2,3,

106

193

1264

465

1457

536

2721

1001

,ј

ц
ч
ш

the second one

Цe = [1;1,1,1,5,1,1,9,1,1,13,1,1,17,...]

and the last one is

e-1

[0;1,6,10,14,18,22,26,30,34,38,42,...]

ж
з
и

0,3,

193

2721

1001

49171

18089

1084483

398959

,ј

ц
ч
ш

Notice that, unlike p, those developments are dramatically regular.

2.1 An algorithm using continued fraction

Let p_k/q_k be the convergent of order k of a simple continued fraction, and let x = [a₀;a₁,a₂,...,a_k], we have the following recurrence relations

p_k

a_kp_k-1+p_k-2

q_k

a_kq_k-1+q_k-2

with initial conditions

p_-1

1,p₀ = a₀

q_-1

0,q₀ = 1.

If we use the continued fraction for (e-1)/2, then a_k = 4(k-1)+2 for k > 1 and it's easy to build an algorithm to compute a fast converging sequence to e. With k = 1500, e is given to 10⁴ decimal places, with k = 12000, e is given to 10⁵digits.

3 Using Padé approximant

If we consider the series expansion

s(t) =

Ґ
е
k = 0

s_kt^k,s_k О 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

m_kt^k

n
е
k = 0

n_kt^k

ц
ч
ч
ч
ч
ш

= O(t^m+n+1), t® 0

When such an approximation exists, it's called a Padé approximant of degree (m,n) of the series s. There are many properties associated with those approximants, they are used to compute series with bad convergence properties. Here we will just give the Padé approximant of the e^t function:

m
е
k = 0

(m+n-k)!m!

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

t^k

n
е
k = 0

(m+n-k)!n!

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

(-t)^k

For example here are some approximant of increasing degree for e^t:

2+t

2-t

12+6t+t²

12-6t+t²

120+60t+12t²+t³

120-60t+12t²-t³

...

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

ж
з
и

2+t

2-t

ц
ч
ш

1/t

ж
з
и

12+6t+t²

12-6t+t²

ц
ч
ш

1/t

ж
з
и

120+60t+12t²+t³

120-60t+12t²-t³

ц
ч
ш

1/t

...

The error being respectively O(t²),O(t⁴),O(t⁶),... The famous formula (1+t)^1/t converges with an error O(t). As an example of the efficiency of the other formulas, the last formula for different values of t gives

t = 1,

e » 2.718(3098591549295...)

t = 1/2,

e » 2.71828(22539303576...)

t = 1/4,

e » 2.7182818(350588212...)

t = 1/32,

e » 2.7182818284590(703...)

4 Infinite products

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

e =

ж
з
и

ц
ч
ш

1/2

ж
з
и

6·8

5·7

ц
ч
ш

1/4

ж
з
и

10·12·14·16

9·11·13·15

ц
ч
ш

1/8

ј,

and another similar relation is [10]:

e = 2

ж
з
и

ц
ч
ш

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 (1); let s_n be the partial sum of this sequence

s_n = 1+

+...+

then

s_n+1 = s_n+

(n+1)!

= s_n

ж
з
и

(n+1)!s_n

ц
ч
ш

= s_n

ж
з
и

u_n+1

ц
ч
ш

Clearly u_n+1 is an integer and

u_n

n!s_n-1

u_n+1

(n+1)!s_n = (n+1)n!(s_n-1+

) = (n+1)(n!s_n-1+1)

therefore, if we define the sequence u_n as following

u₁

u_n+1

(n+1)(u_n+1)

we find

e =

Ґ
Х
n = 1

ж
з
и

u_n+1

u_n

ц
ч
ш

= 2.(

)(

326

325

)...

The rate of convergence is exactly the same as (1).

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! (in our example 6). 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 1 (Euler) Let d_nbe the number of derangements of n different objects (d_nis sometime called subfactorial), we have

d_n = n!

ж
з
и

+...+

(-1)ⁿ

ц
ч
ш

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

d_n

ж
з
и

+...+

(-1)ⁿ

ц
ч
ш

this ratio is such as

lim
n® Ґ

ж
з
и

d_n

ц
ч
ш

= 0.3678794411714423215955237...

and for example

d₃

= 0.3333333333...

d₄

= 0.375

d₅

= 0.3666666666...

d₆

265

= 0.3680555555...

d₇

1854

= 0.3678571428...

6 Other limits

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

A_n

1+2+3+...+n

n(n+1)

log(G_n)

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

G_n = ( 1.2.3...n) ^1/n = (n!)^1/n,

then thanks to Stirling's formula, when n becomes large :

A_n

G_n

n+1

2(n!)^1/n

n+1

2(nⁿe^-n(2pn)^1/2) ^1/n

n+1

n(2pn)^1/(2n)

therefore

lim
n® Ґ

A_n

G_n

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

lim
n® Ґ

ж
и

Х
p_i Ј n

p_i

ц
ш

1/n

= e

where the p_i is the sequence of primes (2,3,5,7,11,13,...). The rate of converge is very slow and for example with n = 10⁶, it gives 2.7141645...

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

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.

8 Records of computation

Number of digits	When	Who	Notes

23	1748	L. Euler	The exponential series e = 1+1/1!+1/2!+... was used [1].
137	1853	W. Shanks	See [2].
205	1871	W. Shanks
346	1884	Boorman	See [4].
2010	1949	John von Neumann	John von Neumann and his team [6] used the ENIAC. The result confirmed a previous computation to 808 digits published in 1946.
100,265	1961	D. Shanks and J.W. Wrench, Jr.	Euler's formula was used on a IBM 7090. The computation took 2.5 hours [8]. (See Math. Comp. 23, 679-680 (1969))
10,000,000	1994	R. Nemiroff and J. Bonnell
18,199,978	1997, May	Patrick Demichel
20,000,000	1997, Aug	Birger Seifert	The computation took 182.5 hours on a Pentium at 133 Mhz
50,000,817	1997, Sep	Patrick Demichel	The computation took 714 hours of a HP 9000/778 (160Mhz).
200,000,579	1999	S. Wedeniwski	The method was based on binary splitting.
869,894,101	1999, Oct	S. Wedeniwski	The method was based on a continued fraction expansion of e.
1,250,000,000	1999, Nov	X. Gourdon	The formula used was the exponential series e = е1/n!, the verification was made with e = 1/(е(-1)ⁿ/n!). 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 10	S. Kondo and X. Gourdon	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 16	C. Martin and X. Gourdon	Colin Martin used the same PiFast33 program 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 02	S. Kondo and X. Gourdon	Shigeru Kondo used PiFast33 again. 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.

References

[1]: L. Euler, Introduction à l'analyse infinitésimale (french traduction by Labey), Barrois, ainé, Librairie, (original 1748, traduction 1796), vol. 1, p. 89-90
[2]: W. Shanks, Contributions to Mathematics Comprising Chiefly the Rectification of the Circle to 607 Places of Decimals, G. Bell, London, (1853), p. vii.
[3]: C. Hermite, Sur la fonction exponentielle, C. R. Académie des Sciences, (1873), vol. 77, p. 18-24, 74-79, 226-233, 285-293
[4]: Boorman, (On the value of e), The Mathematical Magazine, (1884), vol. 1, p. 204
[5]: L W. Richardson, The deferred Approach to the Limit, Philosophical Transactions of the Royal Society of London, (1927), serie A, vol. 226
[6]: 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
[7]: K. Knopp, Theory and application of infinite series, Blackie & Son, London, (1951)
[8]: D. Shanks and J.W. Wrench, Jr., Calculation of p to 100,000 Decimals, Math. Comput., (1962), vol. 16, p. 76-79
[9]: F. Le Lionnais, Les nombres remarquables, Paris, Hermann, (1983), p. 47
[10]: 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)
[11]: E. Maor, e: The Story of a Number, Princeton University Press, (1994)
[12]: Brothers, Harlan J. and J.A. Knox, New closed-form approximations to the logarithmic constant e, Math. Intelligencer, (1998)
[13]: J. A. Knox and Harlan J. Brothers, Novel series-based approximations to e, College Math. J., (1998)
[14]: E. W. Weisstein, CRC Concise Encyclopedia of Mathematics, CRC Press, (1999), p. 503-504