The Euler constant : g
g = 0.57721566490153286060651209008240243104215933593992¼
(Postcript version of this page : gamma.ps, pdf version : gamma.pdf. Versions of Formulas on the Euler constant are also available : gammaFormulas.ps, gammaFormulas.pdf.)

1  Introduction

Euler's Constant was first introduced by Leonhard Euler (1707-1783) in 1734 as
g =
lim
n® ¥ 
 æ
è 		1+  1
2
 +  1
3
 +¼+  1
n
 -log(n) ö
ø 		.
(1)
It is also known as the Euler-Mascheroni constant. According to Glaisher [4], the use of the symbol g is probably due to the geometer Lorenzo Mascheroni (1750-1800) who used it in 1790 while Euler used the letter C.
The constant g is deeply related to the Gamma function G(x) thanks to the Weierstrass formula
 1
G(x)
 =xexp(gx)
Õ
n > 0 
 é
ë 		æ
è 		1+  x
n
 ö
ø 		exp æ
è 		-  x
n
 ö
ø 		ù
û 		.
This identity entails the relation
G¢(1)=-g.
(2)
It is not known if g is an irrational or a transcendental number. The question of its irrationality has challenged mathematicians since Euler and remains a famous unresolved problem. By computing a large number of digits of g and using continued fraction expansion, it has been shown that if g is a rational number 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+  1
2
 +  1
3
 +¼+  1
n
 ,
we have the estimation
Hn-log(n)-g ~  1
2n
 .
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. Better estimations are :
 1
2(n+1)
<
Hn-log(n)-g
<
 1
2n
       (Young [13])
0
<
Hn-  log(n)+log(n+1)
2
 -g
<
 1
6n(n+1)
       (Cesaro)
 -1
48n3
<
Hn-log(n+  1
2
 +  1
24n
 )-g
<
 -1
48(n+1)3
 .
       (Negoi)
Application of the third relation with n=100 gives
-0.6127.10-9 < 0.577215664432730-g < 0
and n=1000 gives
-0.6238.10-13 < 0.577215664901484-g < 0.
A 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+  1
2n
 -
å
k ³ 1 
  B2k
2k
  1
n2k
 ,
(3)
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)-  1
2n
 +  1
12n2
 -  1
120n4
 +  1
252n6
 -  1
240n8
 +  1
132n10
 -  691
32760n12
 +  1
12n14
 .
This technique, directly inherited from the definition, can be employed to compute g at a high precision but suffers from two major drawbacks :

2.2.1  Euler's computation

In 1736, Euler used the asymptotic expansion 3 to compute the first 16 decimal digits of g. He went up to k=7 and n=10, and wrote
g = H10-log(10)-  1
20
 +  1
1200
 -  1
1,200,000
 +  1
252,000,000
 -  1
24,000,000,000
 +...
with
H10
=
2.9289682539682539
log(10)
=
2.302585092994045684
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 usually doing two different calculations to check 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 series
g = 1-log(  3
2
 )- ¥
å
k=1 
  (z(2k+1)-1)
4k(2k+1)
 .
(4)
For large values of k we have
z(2k+1)-1 =  1
22k+1
 +  1
32k+1
 +¼ ~  1
22k+1
hence the series converges geometrically :
 z(2k+1)-1
4k
 ~  1
2.16k
 ,
This relation is issued from properties of the Gamma function and a proof is given in the Gamma function essay.
The first partial sums of series (4) are
x0
=
0.5(9453489189183561...)=1-log(  3
2
 )
x1
=
0.577(69681662853609...)=  13
12
 -log(  3
2
 )-  z(3)
12
x5
=
0.57721566(733782033...)
x10
=
0.57721566490153(417...)

2.2.4  Knuth's computation

In 1962, Knuth used a computer to approximate g with the Euler-Maclaurin expansion (3), with the parameters k=250 and n=104. The error is about
ek,n=  B(2k+2)
(2k+2)
  1
n(2k+2)
 »  2(2k+2)!
(2k+2)(2pn)2k+2
 » 10-1272
and Knuth gave 1271 decimal places of g [8].

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 integer in the table is the number of digits of 1/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
This table shows that the Euler-Maclaurin summation could not be reasonably used to compute more than a few thousands of decimal places of g.

2.3  Exponential integral methods

An efficient way to compute decimal digits of the Euler constant is to start from the identity g = -G¢(1) (see (2)) which entails for any positive integer n, after integrating by part the formula
g+log(n)=In-Rn,       In= ó
õ 		n

0 
  1-e-t
t
  dt,    Rn= ó
õ 		¥

n 
  e-t
t
  dt.
Plugging the series expansion of (1-e-t)/t in In, we obtain
In= ¥
å
k=1 
 (-1)k-1  nk
k·k!
 .
The value In is an approximation to g with the error bound Rn=O(e-n). By stopping the summation at the right index, we obtain the following formula which provides an efficient way to approximate the Euler constant :
g = an
å
k=1 
 (-1)k-1  nk
k·k!
 -logn+O(e-n),       a = 3.5911¼
(5)
The constant a is such that nan/(an)! is of order e-n, and satisfies a(log(a)-1)=1. To obtain d decimal places of g with (5), the formula should be used with n @ dlog(10) and computations should be done with a precision of 2d decimal places to compensate cancellations in the sum In. This method was used by Sweeney in 1963 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
å
k=1 
 (-1)k-1  nk
k·k!
 -log(n)-  e-n
n
 n-2
å
k=0 
  k!
(-n)k
 +O(e-2n),   b = 4.32¼
(6)
The constant b is such that b(log(b)-1)=2. 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.
More information about this technique can be found in [12].

2.4  Bessel function method

A better method (see also [12]) is based on the modified Bessel functions and leads to the formula
g =  An
Bn
 - log(n) + O(e-4n),
with
An = an
å
k=0 
 æ
è 		 nk
k!
 ö
ø 		2

 
 Hk,       Bn = an
å
k=0 
 æ
è 		 nk
k!
 ö
ø 		2

 
 ,
where a = 3.5911¼ satisfies a(log(a)-1)=1.
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 =  1
4n
 2n
å
k=0 
  [(2k)!]3
(k!)4(16n)2k
 .
Brent and McMillan in [12] suggest that
g =  An
Bn
 -  Cn
Bn2
 - log(n) + O(e-8n).
(7)
This time, the summations in An and Bn should go up to bn where b = 4.970625759¼ satisfies b(log(b)-1)=3. The error O(e-8n) followed an empirical evidence but the result had not been proved by Brent and McMillan. Formula (7) 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  Collection of formulae for the Euler constant

Integral and series formulae for the Euler constant can be found in Collection of formulae for the Euler constant.

4  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,000,000 1999 Oct. X. Gourdon and P. Demichel Formula (7) 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
Back to mathematical constants and computation


File translated from TEX by TTH, version 3.40.
On 14 Apr 2004, 23:26.