Institut Fourier

Université de Grenoble I

SUR LE CALCUL NUMERIQUE DE LA CONSTANTE D'EULER

par   Jean-Pierre DEMAILLY

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

0.      UN PEU D'HISTOIRE

La constante d'Euler    \gamma =   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    \gamma    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    \gamma   ont été obtenus par K. Mahler [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    \gamma    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


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    \gamma    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    \gamma    à 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 -- \gamma    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 = - \gamma  . 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    \gamma    à    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 $\gamma$

Les valeurs numériques de    \gamma   obtenues par les auteurs mentionnés ci-dessus ont été utilisées pour déterminer les développements en fraction continue de    \gamma    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    \gamma    ou    e

Y

Théorème.     Si    \gamma   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.   -    \gamma    and    exp(\gamma)    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)