Institut Fourier
Université de Grenoble I
SUR LE CALCUL NUMERIQUE DE LA CONSTANTE D'EULER
par   Jean-Pierre DEMAILLY
0.      UN PEU D'HISTOIRE
La constante d'Euler    y =   lim    1+ - +...+ ~ - Logn  ,   dite
n-+ 00	2	*i
parfois constante de Mascheroni,   a fait l'objet d'un grand nombre de travaux et d'évaluations numériques depuis le  XVIIIe siècle.   Elle garde néanmoins encore beaucoup de son mystère aujourd'hui.   Ainsi,  on ne sait toujours pas si    y    est ou non irrationnel,   en dépit de plusieurs tentatives de démonstration,   dont celle de P. Appell  [2]   en 1926 qui avorta par suite d'une erreur matérielle,   et celle de A.  Froda fl4] en   1965   qui   repose   sur  un   critère   d'irrationnalité   encore    incomplètement démontré.   Signalons toutefois que certains résultats de transcendance pour des expressions faisant intervenir    y   ont été obtenus par K. Manier [18].
Nous nous intéresserons ici surtout  à la mise en oeuvre d'algorithmes rapidement convergents qui,   outre l'intérêt calculatoire,   laissent espérer l'obtention de résultats arithmétiques.
La première évaluation de    y    est due naturellement   à Leonhard Euler,   qui obtint    la valeur    0.577218    en 1735  [12],   bientôt étendue par Mascheroni et quelques autres.   En  1781,   Euler détermina la valeur plus précise    0.577215664901532    [13].   Il fut  suivi notamment par CF. Gauss,   avec    22 décimales   exactes,  puis par un certain nombre de mathématiciens anglais du XtXe  siècle.   Le lecteur pourra consulter J.W.L.   Glaisher  [15]  pour l'historique détaillé des calculs antérieurs à 1870.     YY.  Shanks  fl9]  publia    110 décimales,   dont.    10]    exactes, en 1867-71   ; peu après,   le célèbre mathématicien-astronome
Institut Fourier - B. P. 74  - 3S402 Saint-Martin-d'Hères Cedex
Laboratoire de Mathématiques Pures associé au CNRS - Tél.   (76)   51-46-00
J.C.   Adams fl]   calcula laborieusement    263 décimales,     record    qui devait tenir depuis sa publication en 18 78 jusqu'à l'apparition des premiers ordinateurs et les    328 décimales    de J.W.   Wrench Jr  [23]   en 1952.
Tous ces calculs,   ainsi que celui ultérieur de  D. E.   Knuth fl7]
en  1962 avec    1271  D ,   reposaient  sur le développement asymptotique
de     1 +...+ - -  Logn    par la formule d'Euler - Mac  Laurin.   Le temps
de calcul prohibitif des nombres de Bernoulli requis dans cette méthode
conduisit Dura W.   Sweeney  [21]   à introduire un nouvel algorithme plus
Logx e     dx .     Sweeney obtint 0 ainsi    3566 D    en 1963,   et sa méthode fut reprise  successivement par
Beyer-Waterman  [3],   [4]   en 1974    (7114 D ,   dont    4879 exactes) et
par R.P.   Brent  [7],   [8]     (20700 D    en  1977).   Enfin en 1980, R.P. Brent
et E. Me Millan  [10]   découvrirent un nouvel algorithme plus performant,
utilisant les fonctions de Bessel,   et calculèrent    30100 D    [9] .   Voici un
bref aperçu des algorithmes évoqués plus haut,   avec analyse comparée
des temps de calcul.
I.       FORMULE D'EULER   - MAC LAURIN
Cette formule sera utilisée sous la forme suivante (cf. Bourbaki
t5] )  :
Lf(i> = J f(x)dk +|(f(n)+f(l))
M	"l	2
k     bQi   r /<M
f<2J-D(n).f(2J-D(1)
j=i «m l
+ Rk
où    b.    est la suite des nombres de Bernoulli,   définie par
3    z     -SÎL-J
Z
eZ-l       J-0  * ' de  sorte que
b0- 1 , b = -\ ,  b2=l,  b3= o... .
m	j=0U / J
1 {x}    la partie fractionnaire de    x .   Si nous posons    f(x) = -   ,   la formule ci-dessus entraîne d'après D. Knuth  [17]   :
= Iogn+-+-  +T 11--    +...+ -  1--    -J	2k+2      dx.
n"/	""   v    n
Quand    n - +» ,   il vient
1      b2	b2k       n+»B2k+l<W>
Y	= 2+T   +~+wJ,	2k+2       dX  '
1	X
Par soustraction,  nous obtenons donc :
,i	i     T	i      b2	b2k      r+»-WW>
v = 1 + - +...+ - -  Logn	+	 +...+ 	   -		
Y	2	n	ë	2n      _ 2	Q1   2k      J	2k+2
2n	2kn	n	x
Ce développement asymptotique diverge quand    k -*- +½ ,   mais il donne cependant de bonnes valeurs approchées de    y    lorsque    n    et    k    sont bien choisis.   En effet,  l'identité classique
B_1({x})  = 2(-l)k-1(2k+l"   £ -ta*?
entra me
lB2k+l<W>l   S4f^    '
(2]7)
et grâce à la formule de  Stirling,  nous en déduisons

+. B2k+1(M)dx
n	/k+2
k/   k   >k
Le  reste est donc très petit tant que    k    demeure  sensiblement inférieur à   n  .
Analyse du temps de calcul.
Désignons par    E      cet algorithme et par    E  (d)    le temps de calcul de    y    à la précision    10      .On attribue par convention une
unité de temps à chaque opération arithmétique élémentaire  (portant sur 1 mot de la machine).   La somme de 2 nombres ayant    d    décimales nécessite donc,   à des constantes près que nous négligerons ici,    d   unités de temps ; idem pour les produits ou quotients de nombres en mul-tiprécision par des "petits" nombres (1 mot-machine).
L'évaluation la plus brutale de    1 +-+...+ -    réclame alors    dn
2	n
unités de temps.   Dans la pratique,  on choisira pour    n    une puissance
2 de    2    ou    10  ,  et    Logn    sera évalué en    d      unités de temps à partir
des séries entières    Log2 = 2Argth -   ,     Log -  = 2Argth-     exponen-
3	o	y
tiellement convergentes.
b2k D'autre part,   les nombres    B      = -^r   peuvent être calculés au
n
moyen de la formule de récurrence des nombres de Bernoulli :
1
k-è      k^ f2k+l\    P2j 2k      2k+l[n2k      jti l  2j  ;n2(k-j)J
L'étape de récurrence exige environ    kd    unités de temps,   à condition
3 On
/"2k+l\      2i de  stocker en mémoire les résultats partiels  [?     .   --~rr    ,   soit
V 2J /n2(k"J) 2 total    k d    pour l'évaluation de la somme

b2 2n2
+ .
b2k "+ 2kn2k
E1(d) ~
Bien entendu,    n    et    k    doivent être choisis de manière que le terme d'erreur soit    < 10      ,  ce qui impose
k Log-r-  =*  d LoglO .
2 Le choix optimal est obtenu pour    n ~ k    ,   d'où    kLogk ~ d LoglO ,
^disgio   et
Log d E2(d) ~ Cd3(Logd)'2  .
Signalons qu'il existe une formule  (peu connue) permettant de contourner
le calcul des    b_.   .   Cette formule exprime le reste  sous forme de série 2k
double rapidement convergente :
+ 00        +,
y - X ? ± ?-.+ £ - logn ? i  +   £    S   ^rA(A+i)   >(2kn+{)
Pour vérifier cette identité,  on part de l'intégrale convergente
I(p,q) =   f1xP_1f-1--T2-ldx    ,  où    p.q > 0  .
0	Vi-x^     1 XJ
Pour tout entier    n > 1   ,   un calcul aisé donne
I(n,n) = J (l + x+...+x      )dx - d  Log-p
0
de  sorte que    I(n,n) = 1 + - +...+ -- -  Logn -- y    quand    n -- +½ .
Le changement de variable    x = tr    dans    I(p,q)	fournit d'autre part
l'identité
I(p,q) + I(pr,r) = I(pr,qr)  .
Appliquons cette formule par récurrence avec    p = q  ,     r = 2  .  Il vient
k	k      k
I(n,n)  + I(2n,2)  +...+ 1(2 n,2) = 1(2 n,2 n)   ,
d'où à la limite quand    k -- +o=    :
Y = I(n,n)  +   E I(2kn,2) k=l
TCO	J
= 1 + - +...+ --.-  Logn  +   £
2	n-1	6	k=1J0
La formule annoncée s'en déduit maintenant par développement en série de    1/(1+ x)    :
1 + x	2 l       2
En particulier, pour   n = 1   ,  on obtient la formule  simple
oo	co
<E,>   v = \ + s. s
2      k=l e=l  2f+12k(2k+D...(2k+
6)
Le terme général de la série est majoré par    min(2   "	, ( £2    ) )    ;
la précision    10        est donc atteinte en  sommant pour
k,££dlog10/Log2,        ^k^Vlge     -
ce qui laisse environ	d Logd    termes à calculer.   Le temps de
calcul requis par    (E )    admet donc l'estimation
E2(d) ~ C d2 Logd  ,
asymptotiquement inférieure  à celle de l'algorithme   (E ) .   Nous allons voir qu'il existe en fait des algorithmes simples en    0(d2)   ,   mais cela n'empêche pas l'algorithme      (E )   d'être  sans doute le plus efficace pour les calculs manuels    (d<100) .
2.       ALGORITHME DE SWEENEY
Le point de départ en est l'égalité
r'(l) = J+0° Logt e_tdt = - y  . 0 qui découle de l'identité des définitions d'Euler et de  Weierstrass de
la fonction    T    :
r(l+x) =	t e   dt = e	I   I e      1+-
0	n=l       V    n
Grâce à une intégration par parties,  on obtient Y = F(x) - Logx - R(x)
avec	,
F(x). r i^_dt = z tR-^
0      '	n=l     nn!
x	>?	X	J	X      t
La méthode    (S )    la plus simple,  utilisée par Sweeney [21]   et Beyer-Waterman  f3î,  consiste à choisir un entier    x    assez grand de manière
que    R(x)    soit négligeable et à calculer la valeur approchée Y =*  F(x) - Logx .   Comme    ( obtenue pour    x^ d LoglO  .
Etant donné    p > 0  ,   soit    a      l'unique racine réelle    > 0   de
l'équation
VLog ap "1} = P '
de  sorte que
a    = e- 2.718   ,       a    - 3.591   ,       a    =* 4.319  ,       a   - 4.971   .
U	X	<£	O
x°      /ex\ La formule de Stirling montre que    -- -   -I       est de l'ordre de    e
pour    n = a x .   La  série    F(x)    doit donc être  sommée ici jusqu'à
n  as a x .
Le calcul de chaque terme de la  série,   après factorisation suivant la règle de Horner
réclame 3 opérations arithmétiques  (1 multiplication,   1 division,   1 soustraction).   Une difficulté supplémentaire  se présente ici du fait qu'une partie des chiffres significatifs est perdue par compensation des termes
de signes opposés.   Le terme de valeur absolue maximale étant de l'ordre
xn        x de   -- ** e      pour    n = x ,   ceci amène à travailler avec    2d    décimales n!
au lieu de    d .  On obtient en définitive le temps de calcul suivant  (en négligeant le calcul de    Logx )  :
S (d) = a   x d LoglO x 3 x  2d   =   6a   Log 10 d2 -   49. 6 d2 .
Une méthode    (S1)    plus élaborée,   suggérée par Sweeney f2l]   et mise en oeuvre par Brent  [7],   consiste  à évaluer le  reste    R(x)    par son développement as}rmptotique limité à l'ordre    k = x .   Compte tenu que
««-VK-^
e~Xx!	/2tT   -2x
s --	<; J- e
x+1	v  x
X
on est amené à choisir    x = - d LoglO  ,   à sommer la  série    F(x)
jusqu'à   n = a x   (a   =* 4.319)   ,   et à travailler avec des nombres ayant
3d	x
~    chiffres décimaux.   Le calcul de    R(x)    ou de    e      nécessite 2 opé-
rations pour chaque terme,   effectuées à la précision    fO^2  .  On en
déduit
si(d) =(!a2 + lao + l)Logl0d2 " 26'7d2 -
Il existe deux autres alternatives pour le calcul de     F(x)    qui évitent la perte de précision intervenant dans l'algorithme    (S )  .   Elles reposent sur les développements en séries à termes positifs ci-dessous :
x   x     t	»   .	.   n
<S")        F(x)  = e"Xr   ^-^-dt = e"X Z    l+^+...+ i]^r
2'	'0   x-t	n=1V      2	n/n!
h    x/2      x/2      t
-x/2
= 2e-/2 jLi^J-if^2""1^)
nfA     3	2n-lA (2n-l) !	(2n) !    j
La série    (S )  ,   et de façon analogue    (S )   ,  peut être évaluée par
11 factorisation suivant la règle de Horner  (nous posons    H   = 1+-+...+ - )
?^	n	2	n
N	n	N-l,	.   n
-x ^-       x	"	-x   ^       1	,1    x
= h    -e"XL   Xf-L+...+ i+   x   '
N	V*"  nin+1   '"    N     n+1 V"  N-i{N)'"
ce qui nécessite 1 multiplication,     1 division    et    2 additions    à chaque étape.   Suivant que l'on néglige ou non le reste    R(x)    (les algorithmes avec reste  seront notés    (S')    et    (S') ),   ces méthodes conduisent aux
Z	o
temps de calcul
S (d) = 6a0Logl0d2 -   37. 6d2  , S3(d) =i-aiLogl0d2 -   22. 7 d2  ,
S£(d) =  (3ai+|]lx)gl0d2 -   26. Od2   , S^(d) = (y a3+|)Logl0d2 -   16. 9 d2
3.       ALGORITHME DE BRENT - MAC MILLAN
Cet algorithme repose  sur certaines identités vérifiées par les fonctions de  Bessel modifiées    I (x)    et    Kfi(x)    :
+ 00	"
aW == nfo  n!T(a+n+l)
Kn(x) = - ~^-
a=0    '
les spécialistes observeront que nous avons substitué    2x    à    x   dans
les notations classiques du traité de  Watson  [22].   Par dérivation de    I
a
il vient
KQ(x) = - (Logx + Y)I0(x)  *  S0(x)
+½  v2n	+» /      -,	-, x   2n
où    I (x) -    I £-    ,       S (x) -   £    1+!+...+ £ P«   -
0	n=0n!2	°	n=lV      2	n'n!2
Un calcul faisant intervenir l'intégrale de Hankel de la fonction   -
donne   par ailleurs
T / x      1 r»77       /o           x            j        sinan P+c0       . _     ,    .   -av.
I (x) = - I   exp(2xcosu)cosaudu	exp(-2xchv)e       dv .
a	ttoo	n      _0
En particulier,  on obtient les expressions intégrales et les équivalents suivants de    I (x)   ,     K (x)    quand    x -*? +co    :
T             1 rn   ^xcosu,                       1       zx        .	n
In(x) = -      e	du ~ et >   -	e	si    x ^ 1   ,
0	TT ?'
^ / x       P^00   -*xcav.	,     ^    /TT
KQ(x) = J       e	dv ~   et    <  ,/-
La constante d'Euler peut donc s'écrire  sous la forme
S0(x)	k0(x)
K (x)	_4x
avec     0 <---  < ne	si    x^l  .
I0(x)
0	K0(x)
Dans l'algorithme    (B)    utilisé par Brent,   le reste	est    pure-
-d	°
ment et simplement négligé  ; la précision    10        est donc atteinte pour
x = -d LoglO ,   et les séries    I"(x)   ,     S (x)    doivent être  sommées
jusqu'à   n = a x .   Le calcul de
2 /	2 /	2 ,	2
0	.2 \       021    *    2 V	2
1     v      2   v      n	(n+1)
nécessite    2 opérations arithmétiques    pour chaque terme,   et celui de
S°(X>.H .£Lij.t.,L *2
I0(x)    -     N      x2 f" "2U+1   -   H    (n+1)2 en exige    4  .   Le temps requis par l'algorithme    (B)    est donc
B(d) = a   x i d LoglO x  6 x d -   12.4d2  .
K0(x)
Comme dans le cas de la méthode de  Sweeney,   le reste    -	     peut
I0(x)
être évalué  au moyen d'un développement  asymptotique.  On  a en effet
1  (x)K_(x) = --	exp(2x(cosu - chv))dudv
(r        0V  '        2n -!-TT<U<n
v>0
+ oo	2TT	jq
+ 00
«= - j      exp(-4xr)drf
0	u   j 1 - re    I
grâce au changement de variable    re     = s in    ---   .   Du développement
en série
1   P '        de	t?     (2k)!      2k
7T-        	rr-   =    12    -\,	r        ,     0 £ r < 1   ,
2nl0    |l-reie|	k=0  24kk!4
on déduit alors le développement asymptotique  (divergent)
k!   (16x)
Le terme général de ce développement passe par un minimum pour
k = 2x ,   et on peut vérifier que le reste correspondant est de l'ordre
e-4x
de grandeur du terme    k = 2x ,   soit  	r-y-   .   La valeur approchée
vSïïx372
K0<X>   ^	1	g	(2k)!3
!0(X)	4x1  (x)2   k=0   (k!)4(16x)2k
-8x
est donc affectée d'une erreur ne dépassant pas    e	.   Cette méthode
amène à choisir    x = - d LoglO    et conduit au temps de calcul
o
B'(d) = (ja3 + iJLogl0d2 -   9.7 d2  ,
optimal parmi tous les algorithmes présentés ici.
Nous terminons en donnant le "hit-parade" de ces algorithmes, rangés par ordre d'efficacité croissante.
Algorithmes
Euler El        1      E2
si
S2
Sweeney Si|   S2
S3
S3
Brent B   1   B'
Temps de calcul/d2
Cd
C Logd
49.6
37.6
26.7
26.0
22.7
16.9
12.4
9.7


(Ix>gd)2

















Signalons qu'il existe des algorithmes théoriquement encore plus rapides, permettant d'évaluer    y    à    10      près    en   0(d(Logd)    LogLogd)    unités de temps . Ces derniers reposent  sur l'utilisation de l'algorithme de multiplication rapide de  Schônhage-Strassen   [20]   et  sur une factorisation par blocs de la  série    F(x)  (resp.   e    ,   I (x) ,   S (x))    ; voir Brent  [6].   De tels algorithmes "rapides"  sont toutefois très difficiles à programmer et ne l'emportent sur les algorithmes "classiques" présentés ici que lorsque d    est très grand.
4.  DEVELOPPEMENT EN FRACTION CONTINUE DE y
Les valeurs numériques de    y   obtenues par les auteurs mentionnés ci-dessus ont été utilisées pour déterminer les développements en fraction continue de    y    et    e    ,   qui sont connus maintenant jusqu'à Tordre    29 000    (Brent - Mac Millan  [il]).   La distribution statistique des réduites successives ne fait apparaître aucune différence significative au niveau de    5 %   par rapport à la loi de Gauss - Kusmin (cf. Khintchine  fl6])  :
fa = log2H)   -lo^{1 + lk donnant la fréquence des réduites égales  à   n    dans le développement de presque tout nombre réel.  On obtient d'autre part le résultat  sui-
vant,   qui rend extrêmement  improbable la rationalité de    y    ou    e
Y
Théorème.     Si    y   ou    eY - P/Q   pour des entiers    P,Q   positifs,   alors    Q > io15000 .
BIBLIOGRAPHIE
fl]       J.C.   ADAMS.   - On the value of Euler's constant  ;
Proc.    Roy.   Soc.   London,   v.27,   1878,  p.   88-94. Voir aussi    v.42,   1887,  p.   22-25.
[2]        P.   APPELL.   - Sur la nature  arithmétique de la constante d'Euler ; Comptes Rend.   Ac.   Se.   Paris,   v.182,   12 avril  1926, p.   897-899  ;     19 avril  1926,   p.   949.
[3]      W.A.   BEYER & M. S.   WATERMAN.   - Error analysis of a compu-tation of Euler's constant  ; Math,   of Comp.,   v.28,   1974, p.   599-604.
[4]        W.A.   BEYER & M. S.  WATERMAN.   - Décimais and partial quotients of Euler's constant  and £n 2    ;
UMT 19,   Math,  of Comp.,   v.28,   1974,   p.   667.   Errata : Math,  of Comp.,   MTE  549,   v.32,   1978,  p.   317-318.
[5]       N.   BOURBAKI.   - Fonctions d'une variable réelle  ; chap.   VI,   § 1,  n° 7  ; Hermann,   Paris,   1951.
[6]        R.P.   BRENT.   - The complexity of multiple-précision arithmetic  ;
Complexity of Computational Problem Solving (R. S.   Anderssen and R.P.   Brent,   Editors),   Univ.   of Queensland Press, Brisbane,   1976,  p.   126-165.
[7]       R.P.   BRENT.   - Computation of the regular continued fraction for Euler's constant  ; Math,   of Comp.,   v.31,  July 1977,  p.   771-777.
[8]        R.P.   BRENT.   -    y    and    exp(y)    to    20700 D    and their regular continued fractions to    20 000   partial quotients ; UMT 1,   Math,   of Comp.,   v.32,   1978,  p. 311.
[9]        R.P.   BRENT.   - Euler's constant and its exponential to    30,100 décimais ; Math,   of Comp.,   UMT File.
[10]     R.P.   BRENT  & E.M.   Me MILLAN.   - Sbme new algorithms for high-precision computation of Euler's constant  ; Math,  of Comp.,   v.34,  January  1980,  p. 305-312.
fil]      R.P.   BRENT & E.M.   Me MILLaN.   - The first     29,000    partial
quotients in the regular continued fraction for Euler's constant and its exponential  ; Math,  of Comp.,   UMT File.
r 12]      L.   EULER.   - De progressionibus harmonieis observationes  ;
Euleri Opéra Omnia,   Ser.  1,   v.14,   Teubner,   Leipzig and Berlin,   1925,  p.  93-100.
[13]      L.   EULER.   - De  summis serierum numéros Bernoullianos invol-ventium ;
Euleri Opéra Omnia,   Ser.  1,   v.15,   Teubner,   Leipzig and Berlin,   1927,  p.  91-130.   Voir en particulier p.  115.   Le calcul est donné en détail p. 569-583.
f 14]      A.   FRODA.   -  La constante d'Euler est irrationnelle  ;
Atti Accad.   Naz.   Lincei Rend.   Cl.   Sci.   Fis.   Mat.   Natur. (8)  38   (1965),   p. 338-344.
[15]     J.W. L.   GLAISHER.   - History of Euler's constant ; Messenger of Math.,   v.l,   1872,  p. 25-30.
f 16]      A. Ya.   KHINTCHINE   (A. Ja. HINCIN).   -  Continued Fractions  ; 3e édition,  traduction anglaise par P. Wynn,   Noordhoff, Groningen,   1963.
[17]      D.E.   KNUTH.   - Euler's constant to  1271 places  ; Math,   of Comp.,   v. 16,   1962,  p.  275-281.
[18]     K.   MAHLER.   - Applications of a theorem  by A.B. Shidlovskii ; Proc.   Roy.   Sbe.   London,   A 305  (1968),  p.   149-173.
[19]     W. SHANKS.   - On the numerical value of Euler's constant ; Proc.   Roy.   Soc.   London,   v.15,   1867,  p. 429-432  ; v. 20,   1871,  p.  29-34.
[20]      A.   SCHÔNHAGE  & V. STRASSEN.   - Schnelle    Multiplikation grosser  Zahlen  ; Computing,  v.7,   1971,  p. 281-292.
[21]      D.W.   SWEENEY.   - On the computation of Euler's constant  ; Math,   of Comp.,  v. 17,   1963,  p.  170-178.
[22]      G.N.   WATSON.   - A Treatise on the Theory of Bessel  Functions ; 2nc* édition,   Cambridge Univ.  Press,   London,   1944.
[23]     J.W.   WRENCH,  Jr.   - A new calculation of Euler's constant  ;
Math.   Tables  & other Aids to Comp.,   v.6,   1952,  p. 255.
Ouvrages de référence  sur la fonction    r -
Î24]     E.   ARTIN.   - The  Gamma  Function ;
Holt,   Rinehart and Winston,   1964,  traduit de  :  Einfuhrung in die  Théorie der Gammafunktion,   Teubner,   1931.
[25]     R.   CAMPBELL.   -  Les intégrales eulériennes et leurs applications ; Dunod,   Paris,   1966.
(juin  1984)