Путь: Математика » Быстрые вычисления » The Euler constant gamma
The Euler constant gamma

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

Статья любезно предоставлена:
web: numbers.computation.free.fr

g = 0.57721566490153286060651209008240243104215933593992ј

1.  Introduction

Euler's Constant was first introduced by Leonhard Euler (1707-1783) in 1734 as the limit

 g = lim n® Ґ жз и 1+ 12 + 13 +ј+ 1n -log(n) цч ш .
(1)

It is also known as the Euler-Mascheroni constant (according to Glaisher [4], the symbol g is probably due to the geometer Lorenzo Mascheroni (1750-1800) who used it in 1790 while Euler used the letter C).

g is deeply related to the Gamma function G(x) thanks to the Weierstrass formula

 1G(x) = xexp(gx) Х n > 0 йк л жз и 1+ xn цч ш exp жз и - xn цч ш щъ ы ,

and from this formula follows the relation

 Gў(1) = -g.
(2)

We don't know if g is an irrational or a transcendental number. The question of it's irrationality (or not ?) has challenged mathematicians since Euler, it remains a famous unresolved problem. By computing sufficiently digits of g it has been shown that if g is rational ( = p/q) then the denominator q must have at least 242080 digits.

Even if g is less famous than the constants p and e, it deserves a great attention since it plays an important role in Analysis (Gamma function, Bessel functions, exponential-integral, ...) and occurs frequently in Number Theory (order of magnitude of arithmetical functions for instance [11]).

2.  Computation of the Euler constant

### 2.1  Basic considerations

Direct use of formula (1) to compute Euler constant is of poor interest since the convergence is very slow. In fact, using the harmonic number notation

 Hn = 1+ 12 + 13 +ј+ 1n ,
we have the estimation
 Hn-log(n)-g ~ 12n .
This estimation is the first term of an asymptotic expansion which can be used to compute effectively g, as shown in next section. Nevertheless, other formulae for g (see next sections) provide a simpler and more efficient way to compute it at a large accuracy. The estimation can be refined as :

 12(n+1)
 <
 Hn-log(n)-g < 12n Young [13],
 0
 <
 Hn- log(n)+log(n+1)2 -g < 16n(n+1) Cesaro

and

 -148n3 < Hn-log(n+ 12 + 124n )-g < -148(n+1)3 Negoi

for example if we apply the last relation with n = 100

 -0.6127.10-9 < 0.577215664432730-g < 0

and with n = 1000

 -0.6238.10-13 < 0.577215664901484-g < 0.

Another similar estimation is given in [6].

### 2.2  Asymptotic expansion of the harmonic numbers

The Euler-Maclaurin summation can be used to have a complete asymptotic expansion of the harmonic numbers. We have (see the essay on Bernoulli's numbers)

 Hn-log(n) » g+ 12n - е k і 1 B2k2k 1n2k ,
where the B2k are the Bernoulli numbers. Since B2k grows like 2(2k)!/(2p)2k, the asymptotic expansion should be stopped at a given k. For example, the first terms are given by

 g = Hn-log(n)- 12n + 112n2 - 1120n4 + 1252n6 - 1240n8 + 1132n10 - 69132760n12 + 112n14 .

This technique, directly inherited from the definition, can be employed to compute g at a high precision but suffers from two major drawbacks :

• It requires the computation of the B2k, which is not so easy ;

• the rate of convergence is not so good compared to other formulas with g.

#### 2.2.1  Euler's computation

In 1736, Euler used the previous relation to compute g to 16 decimal places. He went up to k = 7 and n = 10, and he wrote

 g = H10-log(10)- 120 + 11200 - 11,200,000 + 1252,000,000 - 124,000,000,000 +...

with

 H10
 =
 2.92897
 log(10)
 =
 2.30259

giving the approximation

 g » 0.5772156649015329.

#### 2.2.2  Mascheroni's mistake

During the year 1790, in ''Adnotationes ad calculum integrale Euleri'', Mascheroni made a similar calculation up to 32 decimal places. But, a few years later, in 1809, Johann von Soldner (1766-1833) found a value of the constant which was in agreement only with the first 19 decimal places of Mascheroni's calculation ... Embarrassing !

It was in 1812, supervised by the famous Mathematician Gauss, that a young calculating prodigy Nicolai (1793-1846) evaluated g up to 40 correct decimal places, in agreement with Soldner's value [4].

In order to avoid such miscalculations (see also William Shanks famous error on his determination of the value of p), digits hunters are using different calculations to verify the result.

#### 2.2.3  Stieltjes approach

In 1887, Stieltjes computed z(2),z(3),...,z(70) to 32 decimal places and extended a previous calculation done by Legendre up to z(35) with 16 digits. He was then able to compute Euler's constant to 32 decimal places thanks to the fast converging sequence

 g = 1-log( 32 )- Ґ е k = 1 (z(2k+1)-1)4k(2k+1) ,

for large values of k

 z(2k+1)-1
 =
 122k+1 + 132k+1 +ј ~ 122k+1 hence
 z(2k+1)-14k
 ~
 12.16k .

This relation is issued from properties of the Gamma function and a proof is given in the Gamma function essay.

The first iterates are

 x0
 =
 0.5(9453489189183561...) = 1-log( 32 )
 x1
 =
 0.577(69681662853609...) = 1312 -log( 32 )- z(3)12
 x5
 =
 0.57721566(733782033...)
 x10
 =
 0.57721566490153(417...)

#### 2.2.4  Knuth's computation

With a computer, in 1962, Knuth used the Euler-Maclaurin summation with k = 250 and n = 104. The error was about

 ek,n = B(2k+2)(2k+2) 1n(2k+2) » 2(2k+2)!(2k+2)(2pn)2k+2 » 10-1272

In fact the exact value for g was given to 1271 decimal places [8].

#### 2.2.5  Some numerical results on the error function

To appreciate the rate of convergence of this algorithm we give a table of the approximative number of digits one can find with different values for k and n. This number is given by -log10(ek,n).

 k = 10
 k = 100
 k = 250
 k = 500
 n = 103
 63
 390
 769
 1235
 n = 104
 85
 592
 1272
 2237
 n = 105
 107
 794
 1773
 3239
 n = 106
 129
 996
 2275
 4241

From this table, we see that the Euler-Maclaurin summation is limited to a few thousands decimal places for g.

### 2.3  Exponential integral methods

The relation (2) entails, for all integer N

 g+log(N) = IN-RN,       IN = ух N 0 1-e-tt dt,    RN = ух Ґ N e-tt dt.

The series development of (1-e-t)/t gives

 IN = Ґ е n = 1 (-1)n-1 Nnn·n! .
The bound RN = O(e-N) gives another method to compute g :

 g = aN е n = 1 (-1)n-1 Nnn·n! -log(N)+O(e-N),       a @ 3.59      (A2)

The constant a is such that NaN/(aN)! is of order e-N. To obtain d decimal places of g with (A2), the formula should be used with N @ dlog(10) and computations should be done with a precision of 2d decimal places to compensate for cancellation in the sum for IN. This method was used by Sweeney to compute 3566 decimal places of g [9].

A refinement is obtained by approximating RN by its asymptotic expansion, leading to the formula

 g = bN е n = 1 (-1)n-1 Nnn·n! -log(N)- e-NN N-2 е n = 0 n!(-N)n +O(e-2N),    b @ 4.32.     (A3)
This improvement, also due to Sweeney [9], permits to take N @ d/2log(10) and to work with a precision of 3d/2 decimal places to obtain d decimal places of g.

Notice that RN can be approximated as accurately as desired by using Euler's continued fraction

 eNRN = 1/n+1/1+1/n+2/1+2/n+3/1+3/n+ј
This can be used to improve the efficiency of the technique, but leads to a much more complicated algorithm.

### 2.4  Bessel function method

A better method (see also [12]) is based on the modified Bessel functions and leads to the formula

 g = ANBN - log(N) + O(e-4N),
with
 AN = bN е n = 0 жз и Nnn! цч ш 2 Hn,       BN = bN е n = 0 жз и Nnn! цч ш 2 ,
where b = 4.970625759ј satisfies b(log(b)-1) = 3.

This technique is quite easy, fast and it has a great advantage compared to Exponential integral techniques : to obtain d decimal places of g, the intermediate computations can be done with d decimal places.

A refinement can be obtained from an asymptotic series of the error term. It consists in computing

 CN = 14N N е n = 0 [(2n)!]3(n!)4(16N)2k .
Brent and McMillan in [12] suggest that
 g = ANBN - CNBN2 - log(N) + O(e-8N).
(3)
The error O(e-8N) followed an empirical evidence but the result had not been proved by Brent and McMillan. Paul Zimmermann recently afforded a proof of this bound.

Formula (3) has been used by Xavier Gourdon with a binary splitting process to obtain more than 100 millions decimal digits of g in 1999.

Unlike the constant p with the AGM iteration for instance, no quadratically (or more) convergent algorithm is known for g.

3.   Records of computation

 Number of digits When Who Notes 5 1734 L. Euler He found g = 0.577218. 15 1736 L. Euler The Euler-Maclaurin summation was used [1]. 19 1790 L. Mascheroni Mascheroni computed 32 decimal places, but only 19 were correct. 24 1809 J. von Soldner In a work on the logarithm-integral function. 40 1812 F.B.G. Nicolai In agreement with Soldner's calculation. 19 1825 A.M. Legendre Euler-Maclaurin summation was used with n = 10 [2]. 34 1857 Lindman Euler-Maclaurin summation was used with n = 100. 41 1861 Oettinger Euler-Maclaurin summation was used with n = 100. 59 1869 W. Shanks Euler-Maclaurin summation was used with n = 1000. 110 1871 W. Shanks 263 1878 J.C. Adams Adams also computed the first 62 Bernoullian numbers [5]. 32 1887 T. J. Stieltjes He used a series based on the zeta function. ??? 1952 J.W. Wrench Euler-Maclaurin summation [7]. 1271 1962 D.E. Knuth Euler-Maclaurin summation [8]. 3566 1962 D.W. Sweeney The exponential integral method was used [9]. 20,700 1977 R.P. Brent Brent used Sweeney's approach [10]. 30,100 1980 R.P. Brent and E.M. McMillan The Bessel function method [12] was used 172,000 1993 J. Borwein A variant of Brent's method was used. 1,000,000 1997 T. Papanikolaou This is the first gamma computation based on a binary splitting approach. He used a Sun SPARC Ultra, and the computation took 160 hours. He also proved that if g is rational, its denominator has at least 242080 decimal digits. 7,286,255 1998 Dec. X. Gourdon Sweeney's method (with N = 223 ) with binary splitting was used. The computation took 47 hours on a SGI R10000 (256 Mo). The verification was done with the value N = 223+1. 108*106 1999 Oct. X. Gourdon and P. Demichel Formula ( 3) was used with a binary splitting process. The program was from X. Gourdon and Launched by P. Demichel on a HP J5000, 2 processors PA 8500 (440 Mhz) with 2 Go of memory.

References

[1]
L. Euler, Inventio summae cuiusque seriei ex dato termino generali, St Petersbourg, (1736)

[2]
A.M. Legendre, Traité des Fonctions Elliptiques, Paris, (1825-1828), vol. 2, p. 434

[3]
W. Shanks, (On Euler's constant), Proc. Roy. Soc. London, (1869), vol. 18, p. 49

[4]
J.W.L. Glaisher, History of Euler's constant, Messenger of Math., (1872), vol. 1, p. 25-30

[5]
J.C. Adams, On the value of Euler's constant, Proc. Roy. Soc. London, (1878), vol. 27, p. 88-94

[6]
G. Horton, A note on the calculation of Euler's constant, American Mathematical Monthly, (1916), vol. 23, p. 73

[7]
J.W. Wrench Jr., A new calculation of Euler's constant, MTAC, (1952), vol. 6, p. 255

[8]
D.E. Knuth, Euler's constant to 1271 places, Math. Comput., (1962), vol. 16, p. 275-281

[9]
D.W. Sweeney, On the Computation of Euler's Constant, Mathematics of Computation, (1963), p. 170-178

[10]
R.P. Brent, Computation of the regular continued fraction for Euler's constant, Math. Comp., (1977), vol. 31, p. 771-777

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

[12]
R.P. Brent and E.M. McMillan, Some New Algorithms for High-Precision Computation of Euler's constant, Math. Comput., (1980), vol. 34, p. 305-312

[13]
R.M. Young, Euler's constant, Math. Gazette 75, (1991), vol. 472, p. 187-190