Если Вы найдете статью полезной и интересной - не сочтите за труд,
переведите материал или хотя бы часть его и отправьте на адрес
algolist@manual.ru.
Задать вопросы или просто написать письмо можно
также на странице контактов.
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]).
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 :
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.
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.
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.
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.
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.
More information about this technique can be found in [12].
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
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.
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.
