\magnification=1200
\vsize=20cm\hsize=13.5cm\hoffset=-0.1mm
\parindent=0mm \parskip=6pt plus 2pt minus 2pt

\input macros_ocr.tex

%% Further macros
\let\leq=\leqslant
\let\geq=\geqslant

\def\argtanh{\mathop{{\rm arg}\,{\rm tanh}}}

\def\rl#1{\hbox{\rlap{\hbox{$#1$}}}}
\def\vertspace{{\vrule width 0.4pt height 12pt depth 6pt}}
\def\Vertspace{{\vrule width 0.4pt height 28pt depth 6pt}}
\def\Raise#1{\raise8pt\hbox{$#1$}}
\def\table#1{\baselineskip=0pt\lineskip=0pt
    \halign{$\displaystyle{}##\hfil{}$\hfil&&
    $\displaystyle{}##\hfil{}$\hfil\crcr#1\crcr}}

\font\bigbfgreek=cmmib10 at 12pt

\def\statement{\bf}
\def\frac#1#2{{#1\over #2}}
\def\card{\mathop{\rm card}}

%% main text

\footline{\hfill}
\headline{\ifnum \pageno=1 \hfil\else
\ifodd \pageno {\sevenbf \hfil
Sur le calcul numérique de la constante d'Euler\hfil\folio}%
\else {\sevenbf \folio\hfil Gazette des Mathématiciens, vol.\ {\sevenbf 27}, (1985), p.~113--126\hfil} \fi\fi}

{\sevenrm\baselineskip=9pt
Sur le calcul numérique de la constante d'Euler\\
Gazette des Mathématiciens, vol.\ {\sevenbf 27}, (1985), p.~113--126.
\vskip1.5cm}

\centerline{\hugebf SUR LE CALCUL NUMÉRIQUE}
\bigskip
\centerline{\hugebf DE LA CONSTANTE D'EULER 
*\footnote{}{\sevenrm\baselineskip=9pt
* Le présent article est une version étendue du texte originel paru dans 
la Gazette en 1985, où, pour des contraintes de longueur, la plupart
de calculs concernant l'algorithme de Brent-McMillan avaient été omis.
\vskip-\parskip}}
\bigskip
\centerline{\bf par Jean-Pierre DEMAILLY}
\smallskip
\centerline{\vbox{\hrule width 2cm}}
\vskip1.5cm




\section{0}{Un peu d'histoire}

La constante d'Euler    
$$
\gamma =  \lim_{n\to+\infty} 1+{1\over 2} +\cdots+{1\over n} - \log n,
$$
dite parfois constante d'Euler-Mascheroni, a fait l'objet d'un grand nombre
de travaux et d'évaluations numériques depuis le XVIII${}^{\rm e}$
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 [14] 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
C.F.~Gauss, avec 22 décimales exactes, puis par un certain nombre de
mathématiciens anglais du XIX${}^{\rm e}$ siècle. Le lecteur pourra consulter
J.W.L.~Glaisher [15] pour l'historique détaillé des calculs
antérieurs à 1870.  W.~Shanks [19] publia 110 décimales, dont 101
exactes, en 1867-71 ; peu après, le célèbre mathématicien-astronome
J.C.~Adams [1]   calcula laborieusement 263 décimales, record qui 
devait tenir depuis sa publication en 1878 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 [17]
en 1962 avec 1$\,$271~D,   reposaient  sur le développement asymptotique
de $1+{1\over 2} +\cdots+{1\over n}-\log n$ par la formule 
d'Euler-Maclaurin. Le temps de calcul prohibitif des nombres de 
Bernoulli requis dans cette méthode conduisit Dura W.~Sweeney [21]
à introduire un nouvel algorithme plus efficace, basé sur la formule
$\gamma=-\int_0^{+\infty}\log x\;e^{-x}\,dx$. Sweeney obtint ainsi
3$\,$566~D en 1963, et sa méthode fut reprise  successivement par
Beyer-Waterman [3],~[4] en 1974 (7$\,$114~D, dont 4$\,$879 exactes) et
par R.P.~Brent [7],~[8] (20$\,$700~D en 1977). Enfin en 1980, R.P.~Brent
et E.~Mc~Millan [10] découvrirent un nouvel algorithme plus performant,
utilisant les fonctions de Bessel, et calculèrent 30$\,$100~D [9]. Voici un
bref aperçu des algorithmes évoqués plus haut, avec analyse comparée
des temps de calcul.
\bigskip

\section{1}{Formule d'Euler-Maclaurin}

Cette formule sera utilisée sous la forme suivante (cf.\ par exemple
Bourbaki [5])~:
$$
\sum_{i=1}^nf(i)=\int_1^nf(x)\,dx+{1\over 2}\big(f(n)+f(1)\big)
+\sum_{j=1}^k{b_{2j}\over (2j)!}\Big[f^{(2j-1)}(n)-f^{(2j-1)}(1)\Big]+R_k
$$
où $b_j$ est la suite des nombres de Bernoulli, définie par
$$
{z\over e^z-1}=\sum_{j=0}^{+\infty}~{b_j\over j!}\,z^j,
$$
sorte que
$$
b_0 = 1,\quad b_1 = -{1\over 2},\quad b_2={1\over 6},\quad b_3=0,\quad 
b_4=-{1\over 30},\quad b_{2k+1}=0,~~\ldots~.
$$
Le reste $R_k$ est donné par
$$
R_k={1\over (2k+1)!}\int_1^nB_{2k+1}(\{x\})\;f^{(2k+1)}(x)\,dx,
$$
où $B_m(x)=\sum_{j=0}^m{m\choose j }b_jx^{m-j}$ est le $m$-ième polynôme de
Bernoulli et $\{x\}$ la partie fractionnaire de~$x$. Si nous posons 
$f(x) = {1\over x}$, la formule ci-dessus entraîne (cf.\ par 
exemple D.~Knuth [17])~:
$$
\eqalign{
1+{1\over 2}&+\cdots+{1\over n}\cr
&=\log n +{1\over 2}+{1\over 2n}+{b_2\over 2}\Big(1-{1\over n^2}\Big)+\cdots
+{b_{2k}\over 2k}\Big(1-{1\over n^{2k}}\Big)
-\int_1^n{B_{2k+1}(\{x\})\over x^{(2k+2)}}\,dx.\cr}
$$
Quand $n\to+\infty$, il vient
$$
\gamma ={1\over 2}+{b_2\over 2}+\cdots+{b_{2k}\over 2k}
-\int_1^{+\infty}{B_{2k+1}(\{x\})\over x^{2k+2}}\,dx.
$$
Par soustraction,  nous obtenons donc :
$$
\gamma = 1+{1\over 2}+\cdots+{1\over n} -\log n -{1\over 2n}
+{b_2\over 2n^2}+\cdots+{b_{2k}\over 2k\,n^{2k}}
-\int_n^{+\infty}{B_{2k+1}(\{x\})\over x^{2k+2}}\,dx.
$$
Ce développement asymptotique diverge quand $k\to +\infty$, mais il donne 
cependant de bonnes valeurs approchées de $\gamma$ lorsque $n$ et $k$ 
sont bien choisis. En effet, l'identité classique
$$
B_{2k+1}(\{x\})  = 2(-1)^{k-1}(2k+1)!\;\sum_{r=1}^{+\infty}{\sin 2r\pi x
\over (2r\pi)^{2k+1}}
$$
entraîne
$$
|B_{2k+1}(\{x\})|\leq 4\;{(2k+1)!\over (2\pi)^{2k+1}}\;,
$$
et grâce à la formule de  Stirling,  nous en déduisons
$$
\Bigg|\int_n^{+\infty}{B_{2k+1}(\{x\})\over x^{(2k+2)}}\,dx\Bigg|\leq
{4\over\pi}\Big({k\over n\pi e}\Big)^k.
$$
Le  reste est donc très petit tant que  $k$ demeure  sensiblement 
inférieur à~$n$.
\medskip

{\bf Analyse du temps de calcul}.

Désignons par $E_1$ cet algorithme et par $E_1(d)$ le temps de calcul de
$\gamma$ à la précision $10^{-d}$. 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 
multiprécision par des \lguil\?petits\rguil\ nombres ($1$~mot-machine).
L'évaluation la plus brutale de $1+{1\over 2}+\cdots+{1\over n}$
réclame alors $dn$ unités de temps. Dans la pratique, on choisira pour
$n$ une puissance de $2$ ou $10$,  et $\log n$ sera évalué en $d^2$ 
unités de temps à partir des séries entières $\log 2 = 
2 \argtanh {1\over 3}$, $\log {10\over 8}= 2\argtanh {1\over 9}$
exponentiellement convergentes.

D'autre part,   les nombres $\beta_{2k}= {b_{2k}\over n^{2k}}$ 
peuvent être calculés au moyen de la formule de récurrence des nombres 
de Bernoulli~:
$$
\beta_{2k}={1\over 2k+1}\Bigg[{k-{1\over 2}\over n^{2k}}-\sum_{j=1}^{k-1}
{2k+1\choose 2j}{\beta_{2j}\over n^{2(k-j)}}\Bigg].
$$
L'étape de récurrence exige environ $kd$ unités de temps, à condition
de  stocker en mémoire les résultats partiels ${2k+1\choose 2j}
{\beta_{2j}\over n^{2(k-j)}}$, soit au total $k^2d$ pour 
l'évaluation de la somme
$$
{b_2\over 2n^2}+\cdots + {b_{2k}\over 2kn^{2k}}.
$$
On aboutit donc à~:
$$
E_1(d) \sim nd+d^2+k^2d.
$$
Bien entendu, $n$ et $k$ doivent être choisis de manière que le terme
d'erreur soit${}<10^{-d}$,  ce qui impose
$$
k \log{n\pi e\over k} \simeq d\log 10.
$$
Le choix optimal est obtenu pour $n\sim k^2$, d'où $k\log k\sim d\log 10$,
$k\sim {d\log 10\over \log d}$, et
$$
E_1(d)\sim Cd^3(\log d)^{-2}.
$$
Signalons qu'il existe une formule  (moins connue) permettant de contourner
le calcul des~$b_{2k}$. Cette formule exprime le reste  sous forme de série 
double rapidement convergente~:
$$
\gamma=1+{1\over 2}+\cdots+{1\over n} -\log n-{1\over 2n}
+\sum_{k=1}^{+\infty}\sum_{\ell=1}^{+\infty}
{\ell!\over 2^{\ell +1}2^kn(2^kn+1)\cdots(2^kn+\ell)}.
$$
Pour vérifier cette identité,  on part de l'intégrale convergente
$$
I(p,q) =\int_0^1 x^{p-1}\Big({q\over 1-x^q}-{1\over 1-x}\Big)\;dx,
\qquad\hbox{où $p,q > 0$}.
$$
Pour tout entier $n > 1$, un calcul aisé donne
$$
I(n,n)=\int_0^1(1+x+\cdots+x^{n-1})\;dx - d\bigg[\log{1-x^n\over 1-x}\bigg]\,,
$$
de  sorte que $I(n,n) = 1+{1\over 2}+\cdots+{1\over n} - \log n \to \gamma$
quand $n \to +\infty$. Le changement de variable $x = t^r$ 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
$$
I(n,n) + I(2n,2)  +\cdots + I(2 n,2) = I(2^kn,2^kn),
$$
d'où à la limite quand $k\to+\infty$~:
$$
\eqalign{
\gamma &= I(n,n)  + \sum_{k=1}^{+\infty} I(2^kn,2)\cr
&=1+{1\over 2}+\cdots+{1\over n-1}-\log n+\sum_{k=1}^{+\infty}
\int_0^1 x^{2^kn-1}\;{dx\over 1+x}\,.\cr}
$$
La formule annoncée s'en déduit maintenant par développement en série 
de $1/(1+ x)$~:
$$
{1\over 1 + x} = {1\over 2}\Big(1 - {1-x\over 2}\Big)^{-1} =
\sum_{\ell =0}^{+\infty} {(1-x)^\ell\over 2^{\ell+1}}.
$$
En particulier, pour $n = 1$,  on obtient la formule  simple
$$
\gamma={1\over 2}
+\sum_{k=1}^{+\infty}\sum_{\ell=1}^{+\infty}
{\ell!\over 2^{\ell +1}2^k(2^k+1)\cdots(2^k+\ell)}.
\leqno(E_2)
$$
Le terme général de la série est majoré par 
$\min\big(2^{-\ell-k-1},(\ell 2^{-k})^\ell\big)$~;
la précision $10^{-d}$ est donc atteinte en  sommant pour
$$
k,\ell\leq d \log 10/\log2\leq 4d,\qquad 
\ell\leq {d\log 10\over k\log 2-\log\ell},
$$
ce qui (en majorant $\ell$ par $4d$ pour $k\leq 2\log(4d)/\log 2$
et par $2d\log 10/k\log 2$ sinon) laisse 
$O(d\log d)$ termes à calculer.
Le temps de calcul requis par $(E_2)$ admet donc l'estimation
$$
E_2(d) \sim C d^2 \log d,
$$
asymptotiquement inférieure  à celle de l'algorithme $(E_1)$.
Nous allons voir qu'il existe en fait des algorithmes simples en $O(d^2)$,
mais cela n'empêche pas l'algorithme $(E_1)$ d'être  sans doute le plus 
efficace pour les calculs manuels  (disons $d\leq 100$).
\bigskip

\section{2}{Algorithme de Sweeney}

Le point de départ en est l'égalité
$$
\Gamma'(1) = \int_0^{+\infty} \log t\;e^t\,dt = -\gamma,
$$
qui découle par exemple de l'identité des définitions d'Euler et de 
Weierstrass de la fonction~$\Gamma$~:
$$
\Gamma(1+x) = \int_0^{+\infty} t^x\;e^{-t}\;dt = e^{-\gamma x}
\prod_{n=1}^{+\infty}\Big(1+{x\over n}\Big)^{-1}.
$$
Grâce à des intégrations par parties, on obtient
$$
\gamma = F(x) - \log x - R(x)
$$
avec
$$
\eqalign{
F(x)&=\int_0^x{1-e^{-t}\over t}\,dt\;,\cr
R(x)&=\int_x^{+\infty}{e^{-t}\over t}\,dt = {e^{-x}\over x}\bigg(1-{1!\over x}+
\cdots+{(-1)^kk!\over x^k}\bigg)+(-1)^{k+1}(k+1)!\int_x^{+\infty}
{e^{-t}\;dt\over t^{k+2}}.\cr}
$$
La méthode $(S_1)$ la plus simple, utilisée par Sweeney~[21] et 
Beyer-Waterman~[3],  consiste à choisir un entier $x$ assez grand 
de manière que $R(x)$ soit négligeable et à calculer la valeur approchée 
$\gamma \simeq F(x) - \log x$. Comme $0<R(x)<{e^{-x}\over x}$, la précision
$10^{-d}$ est obtenue pour $x\simeq d\log 10$.

Etant donné $p\geq 0$, soit $a_p$ l'unique racine réelle${}> 0$ de
l'équation
$$
a_p(\log a_p - 1) = p,
$$
de  sorte que
$$
a_0 = e \simeq 2.718,\quad
a_1\simeq 3.591,\quad
a_2\simeq 4.319,\quad
a_3\simeq 4.971.
$$
La formule de Stirling montre que ${x^n\over n!}\simeq \big({ex\over n}\big)^n$ 
est de l'ordre de $e^{-px}$ pour $n\simeq a_px$. La  série $F(x)$ doit donc 
être  sommée ici jusqu'à $n=\lceil a_1x\rceil$.

Le calcul de chaque terme de la  série, après factorisation suivant la règle
de Hörner 
$$
F(x)={x\over 1}\bigg({1\over 1}-{x\over 2}\bigg({1\over 2}-\cdots
{x\over n-1}\bigg({1\over n-1}-{x\over n}\bigg({1\over n}\bigg)\bigg)\cdots
\bigg)\bigg),
$$
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 de ${x^n\over n!}\simeq e^x$ pour $n = x$,   ceci amène à 
travailler avec des nombres à $2d$ chiffres décimaux au lieu de~$d$.
On obtient en définitive le temps de calcul suivant (en négligeant 
le calcul de $\log x$)~:
$$
S_1(d) = a_1\times d \log 10 \times 3 \times 2d   =   6a_1 \log 10\;d^2 
\simeq 49.6\,d^2.
$$
Une méthode $(S'_1)$ plus élaborée, suggérée par Sweeney [21] et mise en 
{\oe}uvre par Brent~[7], consiste  à évaluer le  reste $R(x)$ par son 
développement asymptotique limité à l'ordre $k = x \in\bN$. Compte tenu que
$$
\bigg| R(x) - {e^{-x}\over x}\bigg(1-{1!\over x}+
\cdots+{(-1)^kk!\over x^k}\bigg)\bigg|
\leq {e^{-x}\;x! \over x^{x+1}}
\leq \sqrt{{2\pi\over x}}\;e^{-2x},
$$
on est amené à choisir $x = -{1\over 2}d\log 10$, à sommer la série 
$F(x)$ jusqu'à $n = a_2x$ (\hbox{$a_2\simeq 4.319$}),   et à travailler avec 
des nombres ayant ${3d\over 2}$ chiffres décimaux. Le calcul de $R(x)$ 
ou de $e^x$ nécessite 2 opérations pour chaque terme, effectuées à 
la précision~$10^{-d/2}$.  On en déduit
$$
\textstyle
S_1'(d) = \big({9\over 4}\,a_2 +{1\over 2}a_0 + {1\over 2}\big)
\log 10\;d^2 \simeq 26.7\,d^2.
$$
Il existe deux autres alternatives pour le calcul de $F(x)$ qui évitent la 
perte de précision par compensation intervenant dans l'algorithme $(S_1)$.
Elles reposent sur les développements en séries à termes positifs ci-dessous~:
$$
\leqalignno{
F(x)&=e^{-x}\int_0^x{e^x-e^t\over x-t}\;dt=e^{-x}\sum_{n=1}^{+\infty}
{1\over n!}\int_0^x {x^n-t^n\over x-t}\;dt\cr
&=e^{-x}\sum_{n=1}^{+\infty}\bigg(1+{1\over 2}+\cdots+{1\over n}\bigg)\,
{x^n\over n!}&(S_2)\cr
F(x)&=e^{-x/2}\int_{-x/2}^{x/2}{e^{x/2}-e^t\over x/2-t}\;dt 
= e^{-x/2}\sum_{n=1}^{+\infty}
{1\over n!}\int_{-x/2}^{x/2}{(x/2)^n-t^n\over x/2-t}\;dt\cr
&=e^{-x/2}\sum_{p=1}^{+\infty}2\bigg(1+{1\over 3}+\cdots+{1\over 2p-1}\bigg)\,
\bigg({(x/2)^{2p-1}\over (2p-1)!}+{(x/2)^{2p}\over (2p)!}\bigg).
&(S_3)\cr}
$$
[On utilise le développement ${x^n-t^n\over x-t}=x^{n-1}+x^{n-2}t+\cdots+t^{n-1}$, et pour la dernière sommation, on regarde séparément les termes $n=2p$ et
$n=2p-1$, $p \geq 1$]. La série $(S_2)$, et de façon analogue $(S_3)$,  
peut être évaluée par factorisation suivant la règle de Hörner 
(nous notons $H_n=1+{1\over 2}+\cdots+{1\over n}$ les sommes partielles
de la série harmonique, avec par convention $H_0=0$, et 
$H_{n,N}=H_N-H_n={1\over n+1}+\cdots+{1\over N}$)~:
$$
\eqalign{
e^{-x}\sum_{n=1}^{+\infty} H_n{x^n\over n!}
&=H_N - e^{-x}\sum_{n=0}^{+\infty}H_{n,N}{x^n\over n!}
\simeq H_N-e^{-x}\sum_{n=0}^{N-1}H_{n,N}{x^n\over n!}\cr
&\simeq H_{0,N}-e^{-x}
{x\over 1}\bigg(H_{1,N}+{x\over 2}\bigg(H_{2,N}+{x\over 3}\bigg(\cdots
+{x\over N-1}H_{N-1,N}\bigg)\cdots\bigg)\bigg).\cr}
$$
Le calcul est effectué de manière descendante en prenant $H_{N-1,N}={1\over N}$
et $H_{n-1,N}={1\over n}+H_{n,N}$. Ceci 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'_2)$ et 
$(S'_3)$), ces méthodes conduisent aux temps de calcul
$$
\eqalign{
S_2(d) &= 6a_0\log 10\;d^2 \simeq 37.6\,d^2,\cr
S_3(d) &= {11\over 4}a_1\log 10\;d^2 \simeq 22.7\,d^2,\cr
S'_2(d) &= \Big(3a_1+{1\over 2}\Big)\log 10\;d^2 \simeq 26.0\,d^2,\cr
S'_3(d) &= \Big({11\over 8}a_3+{1\over 2}\Big)\log 10\;d^2 \simeq 16.9\,d^2.\cr}
$$

\section{3}{Algorithme de Brent-McMillan}

L'algorithme introduit par Brent-McMillan dans [10] repose  sur certaines 
identités vérifiées par les fonctions de  Bessel modifiées 
$I_\alpha(x)$ et $K_0(x)$~:
$$
\eqalign{
I_\alpha(x)&
=\sum_{n=0}^{+\infty} {x^{\alpha+2n}\over n!\,\Gamma(\alpha+n+1)},\cr
K_0(x)&=-{\partial I_\alpha(x)\over\partial\alpha}_{|\alpha=0}~.\cr}
$$
Les spécialistes observeront que $2x$ a été substitué à $x$ dans
les notations classiques du traité de  Watson~[22]. On se restreindra
aux valeurs $x>0$ réelles dans ce qui suit. Par dérivation de
$I_\alpha$ et grâce à la formule $\Gamma'(n+1)=(H_n-\gamma)\,n!$ (elle-même
tirée des égalités ${\Gamma'(x+1)\over\Gamma(x+1)}={1\over x}+
{\Gamma'(x)\over\Gamma(x)}$ et $\Gamma'(1)=-\gamma$), il vient
$$
\eqalign{
K_0(x) &= - (\log x + \gamma)I_0(x) + S_0(x)\qquad\hbox{où}\cr
I_0(x) &= \sum_{n=0}^{+\infty} {x^{2n}\over n!^2},\qquad
S_0(x) = \sum_{n=1}^{+\infty} H_n\;{x^{2n}\over n!^2}.\cr}
$$
L'intégrale de Hankel de la fonction $1/\Gamma$ s'exprime 
cette fonction sous la forme
$$
{1\over\Gamma(z)} = {1\over 2\pi i}\int_{(C)}\zeta^{-z}e^\zeta\;d\zeta
$$
où $(C)$ est le contour formé d'un petit cercle $\zeta=\varepsilon e^{iu}$,
$u\in[-\pi,\pi]$, complété par deux intervalles $]-\infty,-\varepsilon]$
d'arguments $-\pi$ et $+\pi$ parcourus en sens inverse. Cette formule
donne
$$
\eqalign{
I_\alpha(x)&=\sum_{n=0}^{+\infty}{x^{\alpha+2n}\over n!}
{1\over 2\pi i}\int_{(C)}\zeta^{-\alpha-n-1}e^\zeta\;d\zeta=
{1\over 2\pi i}\int_{(C)}x^\alpha \zeta^{-\alpha-1}\exp(x^2/\zeta+\zeta)\;d\zeta\cr
&={1\over 2\pi i}\int_{(C)}\zeta^{-\alpha}\exp(x/\zeta+\zeta x)\;d\zeta\cr
&={1\over\pi}\int_0^\pi e^{2x\cos u}\cos(\alpha u)\;du
-{\sin\alpha\pi\over\pi}\int_0^{+\infty}e^{-2x\cosh v}\,e^{-\alpha v}\,dv.
\cr}
$$
La deuxième ligne ci-dessus exprimant $I_\alpha(x)$ est obtenue à l'aide 
du changement de variable $\zeta\mapsto \zeta x$ (rappelons que $x>0$)~;
la première intégrale de la troisième
ligne provient de la partie du contour $\{\zeta=e^{iu}\}$ formée du cercle 
de centre $0$ et de rayon~$1$, et la dernière intégrale des deux segments 
$]-\infty,-1]$ avec $t=-e^{-v}$, $v\in[0,+\infty[\,$.
En~particulier,  on obtient les expressions intégrales et les équivalents 
suivants de $I_0(x)$, $K_0(x)$ quand $x\to +\infty$~:
$$
\leqalignno{
I_0(x)&={1\over\pi}\int_0^\pi e^{2x\cos u}\;du\kern27.5pt
\hbox{de sorte que}~~
I_0(x) \mathop{\sim}\limits_{x\to+\infty}~{1\over\sqrt{4\pi x}}\,e^{2x},&(1)\cr
K_0(x)&=\int_0^{+\infty}e^{-2x\cosh v}\,dv\qquad
\hbox{de sorte que}~~
K_0(x) \mathop{\sim}\limits_{x\to+\infty}~\sqrt{\pi\over 4x}\,e^{-2x}.&(2)\cr}
$$
On a de plus $I_0(x)>{1\over\sqrt{4\pi x}}\,e^{2x}$ si $x\geq 1$ et
$K_0(x) <\sqrt{\pi\over 4x}\,e^{-2x}$ si $x>0$. Ces estimées peuvent se voir
à l'aide de changements de variables
$$
\eqalign{
I_0(x)&={e^{2x}\over 2\pi\sqrt{x}}\int_0^{4x} {e^{-t}\over\sqrt{t(1-t/4x)}}\;dt,
\qquad t=2x(1-\cos u),\cr
K_0(x)&={e^{-2x}\over 2\sqrt{x}}\int_0^{+\infty}{e^{-t}\over\sqrt{t(1+t/4x)}}\;dt,
\qquad t=2x(\cosh v-1),\cr}
$$
en observant que $\int_0^{+\infty}{1\over \sqrt{t}}\,e^{-t}\,dt=
\Gamma({1\over 2})=\sqrt{\pi}$~; la minoration de $I_0(x)$ résulte du
fait que ${1\over \sqrt{1-t/4x}}\geq 1+t/8x$ par convexité, d'où
(au moyen d'une intégration par parties pour le terme $\sqrt{t}\;e^{-t}$)
$$
\eqalign{
\int_0^{4x} {e^{-t}\over\sqrt{t(1-t/4x)}}\;dt&\geq
\Gamma({\textstyle{1\over 2}})+{1\over 8x}\Gamma({\textstyle{3\over 2}})
-\int_{4x}^{+\infty}\Big({1\over\sqrt{t}}+{\sqrt{t}\over 8x}
\Big)\,e^{-t}\,dt\cr
&\geq \sqrt{\pi}+{\sqrt{\pi}\over 16x}
-e^{-4x}\Big({3\over 4\sqrt{x}}+{1\over 32x\,\sqrt{x}}\Big)
>\sqrt{\pi}\cr}
$$
pour $x\geq 1$. En définitive, la constante d'Euler peut s'écrire sous la forme
$$
\gamma={S_0(x)\over I_0(x)}-\log x - {K_0(x)\over I_0(x)},\leqno(3)
$$
avec
$$
0 < {K_0(x)\over I_0(x)} <\pi\;e^{-4x}\qquad\hbox{si $x\ge 1$}.\leqno(4)
$$
Dans l'algorithme basique $(B)$ proposé par Brent, le reste 
${K_0(x)\over I_0(x)}$ est purement et simplement négligé~; la 
précision $10^{-d}$ est donc atteinte pour
$x = {1\over 4}\,d\,\log 10$, et les séries \hbox{$I_0(x)$, $S_0(x)$} doivent 
être sommées jusqu'à $n = \lceil a_1x\rceil$. Le calcul de
$$
I_0(x)=1+{x^2\over 1^2}\bigg(1+{x^2\over 2^2}\bigg(\cdots
{x^2\over (n-1)^2}\bigg(1+{x^2\over n^2}\bigg(\cdots\bigg)\bigg)\cdots
\bigg)\bigg)
$$
nécessite 2 opérations arithmétiques pour chaque terme, et celui de
$$
S_0(x)\simeq H_{0,N} - {1\over I_0(x)}\;
{x^2\over 1^2}\bigg(H_{1,N}+{x^2\over 2^2}\bigg(\cdots
{x^2\over (n-1)^2}\bigg(H_{n-1,N}+{x^2\over n^2}\bigg(H_{n,N}
+\cdots\bigg)\bigg)\cdots\bigg)\bigg)
$$
en exige 4. Le temps requis par l'algorithme $(B)$ est donc
$$
B(d) = a_1 \times {1\over 4}\,d\;\log 10 \times 6 \times d \simeq 12.4\,d^2.
$$

{\bf Raffinement de l'algorithme.}
Comme dans le cas de la méthode de  Sweeney, le reste $K_0(x)/I_0(x)$
peut être évalué  au moyen d'un développement  asymptotique. On a en effet
$$
I_0(x)K_0(x)={1\over2\pi}
\int_{\{-\pi<u<\pi\,,\,v>0\}}\exp\big(2x(\cos u -\cosh v)\big)\,du\,dv,
$$
et le changement de variable 
$$
r\;e^{i\theta}=\sin^2\Big({u+iv\over 2}\Big) = 
{1\over 2}\big(1-\cos(u+iv)\big)=
{1\over 2}\big(1-\cos u\cosh v+i\sin u\sinh v\big),
$$
donne
$$
\eqalign{
&r={1\over 2}(\cosh v-\cos u),\qquad|1-r\,e^{i\theta}|=
\Big|\cos\Big({u+iv\over 2}\Big)\Big|^2,\cr
&r\,dr\,d\theta=\Big|\sin\Big({u+iv\over 2}\Big)
\cos\Big({u+iv\over 2}\Big)\Big|^2\,du\,dv
=r\;|1-r\,e^{i\theta}|\,du\,dv,\cr}
$$
par conséquent
$$
I_0(x)K_0(x)=
{1\over2\pi}\int_0^{+\infty}\exp(-4xr)\;dr\int_0^{2\pi}{d\theta\over
|1-r\,e^{i\theta}|}.\leqno(5)
$$
Notons
$$
{\alpha\choose k}={\alpha(\alpha-1)\cdots(\alpha-k+1)\over k!},\qquad
\alpha\in\bC
$$
les coefficients binomiaux généralisés.
Pour $z=r\,e^{i\theta}$ et $|z|=r<1$, le développement binomial $(1-z)^{-1/2}=
\sum_{k=0}^{+\infty}{-{1\over 2}\choose k}\,(-z)^k$
combiné à la formule de Parseval-Bessel implique le
développement en série$\phantom{\big|}$
$$
\varphi(r):={1\over2\pi}\int_0^{2\pi}{d\theta\over|1-r\,e^{i\theta}|}=
\sum_{k=0}^{+\infty}w_k\,r^{2k},\qquad 
\hbox{$0\leq r<1$,}\leqno(6)
$$
où le coefficient
$$
w_k:={-1/2\choose k}^2=\bigg({1\cdot 3\cdot 5\cdots(2k-1)\over 
2\cdot 4\cdot 6\cdots 2k}\bigg)^2={(2k)!^2\over2^{4k}\,k!^4}.\leqno(7)
$$
est étroitement lié à l'intégrale de Wallis 
$\smash{W_p=\int_0^{\pi/2}\sin^p x\;dx}$. La relation de récurrence
classique $W_p={p-1\over p}W_{p-2}$ implique en effet 
$$
W_{2k}={1\cdot 3\cdot 5\cdots (2k-1)\over 2\cdot 4\cdot 6\cdots 2k}\;
{\pi\over 2},\qquad
W_{2k+1}={2\cdot 4\cdot 6\cdots 2k\over 3\cdot 5\cdots (2k+1)},
$$
d'où $w_k=({2\over \pi}W_{2k})^2$. Les relations
$W_{2k}W_{2k-1}={\pi\over 4k}$, $W_{2k}W_{2k+1}={\pi\over 2(2k+1)}$ entraînent
par décroissance de la suite $(W_p)$ que $\sqrt{\pi\over 2(2k+1)}<W_{2k}<
\sqrt{\pi\over 4k}$, donc
$$
{2\over \pi(2k+1)}<w_k<{1\over \pi k}.\leqno(8)
$$
Le point de départ de notre analyse pour l'estimation de $I_0(x)K_0(x)$ 
est l'égalité tirée de (5), (6)
$$
I_0(x)K_0(x)=\int_0^{+\infty}e^{-4xr}\,\varphi(r)\,dr,
\leqno(9)
$$
où
$$
\leqalignno{
\varphi(r)&=\sum_{k=0}^{+\infty}w_k\,r^{2k}\kern85pt\hbox{pour $r<1$,}
&(9')\cr
\varphi(r)&={1\over r}\varphi\bigg({1\over r}\bigg)=
\sum_{k=0}^{+\infty}w_k\,r^{-2k-1}\qquad\hbox{pour $r>1$.}
&(9'')\cr}
$$
(Cette dernière identité s'obtient en faisant le changement de $\theta$ en 
$-\theta$ dans~(6)). D'autre part, on voit facilement à l'aide de (8) 
que l'on a un équivalent
$$
\varphi(r)\sim\sum_{k=1}^{+\infty}{r^{2k}\over\pi k}=
{1\over \pi}\,\log{1\over 1-r^2}\qquad\hbox{quand $r\to 1-0$},
$$
en particulier l'intégrale (9) converge autour de $r=1$
(il nous sera nécessaire plus loin de préciser davantage
cette approximation, mais il faut des arguments bien plus sophistiqués
pour y parvenir).
Par intégration terme à terme de la série donnant $\varphi(r)$ sur
$[0,+\infty[\,$, en ignorant le fait que la série diverge pour $r\geq 1$, 
on obtient formellement un développement asymptotique, lui aussi divergent~:
$$
I_0(x)K_0(x)\sim \sum_{k\in\bN}w_k\,{(2k)!\over (4x)^{2k+1}} \sim
{1\over 4x}
\sum_{k\in\bN}{(2k)!^3\over k!^4\,(16x)^{2k}}.
\leqno(10)
$$
Si $x$ est entier, le terme général de ce développement passe par un 
minimum exactement pour $k = 2x$, car le rapport entre le $k$-ième et
le $(k-1)$-ième terme vaut
$$
{(2k(2k-1))^3\over k^4\,(16x)^2}=\bigg({k\over 2x}\bigg)^2
\bigg(1-{1\over 2k}\bigg)^3<1\quad\hbox{ssi $k\le 2x$}.
$$
L'idée est de tronquer le développement asymptotique précisément à $k=2x$.
On va vérifier, comme conjecturé par Brent-McMillan [10] et en partie
démontré par Brent-Johansson [$10'$], que 
\lguil\?l'erreur de sommation\?\rguil\ correspondante 
est alors d'un ordre de grandeur comparable à celui du dernier 
terme $k = 2x$, soit $\smash{e^{-4x}\over 2\,\sqrt{2\pi}\;x^{3/2}}$
par la formule de Stirling.
\medskip

{\statement Théorème.} {\it La différence
$$
\Delta(x):=
I_0(x)K_0(x)-{1\over 4x}\sum_{k=0}^{2x}{(2k)!^3\over k!^4\,(16x)^{2k}}
\leqno(11)
$$
admet quand $x\to+\infty$ l'équivalent
$$
\Delta(x) \sim -{5\,e^{-4x}\over 24\sqrt{2\pi}\,x^{3/2}},
\leqno(12)
$$
et on a plus précisément
$$
\Delta(x)=-e^{-4x}\bigg(
{5\over 24\sqrt{2\pi}\,x^{3/2}}+\varepsilon(x)\bigg),
\qquad|\varepsilon(x)|<{0.863\over x^2}.
\leqno(13)
$$
La valeur approchée
$$
{K_0(x)\over I_0(x)}\simeq {1\over 4x\,I_0(x)^2}
\sum_{k=0}^{2x}{(2k)!^3\over k!^4\,(16x)^{2k}}
$$
est donc affectée d'une erreur dont l'ordre de grandeur est 
$$
{\Delta(x)\over I_0(x)^2}\sim \smash{-{5\sqrt{2\pi}\over 12\,x^{1{/}2}}\,e^{-8x}}.
\leqno(13')
$$}%

La vérification de ce résultat demande de nombreux calculs. Comme on le 
verra plus loin, on pourrait probablement obtenir en fait un développement
asymptotique de~$\Delta(x)$, au moins pour les premiers termes, mais les 
calculs nécessaires seraient sans doute assez longs. D'après (5) et la 
définition de~$\Delta(x)$, nous avons
$$
\Delta(x)=\int_0^{+\infty}e^{-4xr}\,\delta(r)\,dr
\leqno(14)
$$ 
où
$$
\delta(r):=\varphi(r)-\sum_{k=0}^{2x}w_k\,r^{2k},\qquad\hbox{de sorte que}~~
\delta(r)=\sum_{k=2x+1}^{+\infty}w_k\,r^{2k}\quad\hbox{pour $r<1$}.
\leqno(15)
$$
Pour $r<1$, on notera que $\varphi(r)$ coïncide avec l'intégrale elliptique
de première espèce
${2\over\pi}\int_0^{\pi/2}(1-r^2\sin^2\theta)^{-1/2}\,d\theta$, comme il 
résulte à nouveau de la série binomiale et de l'expression de~$W_{2k}$.
Nous aurons besoin de connaître le comportement asymptotique
précis de $\varphi(r)$ quand~$r\to 1$. Celui-ci résulte d'une identité 
classique que nous
rappelons ci-dessous. En posant $t^2=1-r^2$, le changement de 
variable $u=\tan\theta$ donne
$$
\leqalignno{
\varphi(r)&={2\over\pi}\int_0^{\pi/2}(1-r^2\cos^2\theta)^{-1/2}\,d\theta
={2\over\pi}\int_0^{+\infty}{du\over\sqrt{(1+u^2)(t^2+u^2)}}\,du\cr
&={4\over\pi}\int_0^1{dv\over \sqrt{(1+v^2)(1+t^2v^2)}}
+{2\over\pi}\int_t^1{dv\over \sqrt{(1+v^2)(t^2+v^2)}}&(16)\cr}
$$
où la dernière ligne est obtenue en découpant l'intégrale
$\int_0^{+\infty}\ldots\,du$ sur les 3 intervalles $[0,t]$, $[t,1]$,
$[1,+\infty[$ et en y effectuant les changements de variables respectifs 
$u=vt$, \hbox{$u=v$}, $u=1/v$ (le premier et le troisième morceaux étant égaux).
Grâce à la série binomiale, la première intégrale de la ligne (16) admet 
un développement en série convergent
$$
{4\over\pi}\int_0^1{dv\over \sqrt{(1+v^2)(1+t^2v^2)}}={4\over\pi}
\sum_{k=0}^{+\infty}
c'_kt^{2k},\qquad 
c'_k={-1/2\choose k}\int_0^1{v^{2k}\,dv\over\sqrt{1+v^2}}.
$$
La deuxième intégrale s'obtient de même à l'aide d'une série double, en
développant simultanément les deux racines carrées~:
$$
{2\over\pi}\int_t^1{dv\over v\sqrt{1+v^2}\,\sqrt{(1+t^2/v^2)}}
={2\over\pi}\int_t^1\sum_{k,\ell\geq 0}{-1/2\choose\ell}v^{2\ell}\,
{-1/2\choose k}(t^2/v^2)^k\;
{dv\over v}.
$$
Les termes diagonaux $k=\ell$ fournissent un terme logarithmique
$$
{2\over\pi}\sum_{k=0}^{+\infty}
{-1/2\choose k}^2\;t^{2k}\;
\log{1\over t}={1\over\pi}\,\varphi(t)\,\log{1\over t^2},
$$
et les autres s'intègrent sous la forme d'une série double absolument 
convergente
$$
{2\over\pi}\sum_{k\neq\ell\geq 0}{-1/2\choose k}{-1/2\choose\ell}t^{2k}
\Bigg[{v^{2\ell-2k}\over 2\ell-2k}\Bigg]_t^1
={2\over\pi}\sum_{k\neq\ell\geq 0}{-1/2\choose k}{-1/2\choose\ell}\;
{t^{2k}-t^{2\ell}\over 2(\ell-k)}.
$$
Après regroupement des puissances  de $t$, celle-ci se réduit à une 
série entière ${4\over\pi}\sum c''_kt^{2k}$ de rayon de convergence $1$, avec
(par raison de symétrie en $k,\ell$)
$$
c''_k=\sum_{0\leq\ell<+\infty,\,\ell\neq k}~{1\over2(\ell-k)}
{-1/2\choose k}
{-1/2\choose\ell}.
$$
D'après (8), on voit en effet a priori que
$$
|c'_k|\leq {1\over\sqrt{\pi k}}\,{1\over 2k+1}=O(k^{-3/2}),
\kern90pt
$$
et
$$
|c''_k|\leq{1\over 2\sqrt{\pi k}}\Bigg({1\over k}+
\sum_{0<\ell\neq k}{1\over|\ell-k|\,\sqrt{\pi\ell}}\Bigg)
=O\bigg({\log k\over k}\bigg).
$$
Au total, en posant $t^2=1-r^2$, on obtient la relation
$$
\varphi(r)={1\over\pi}\bigg(\varphi(t)\,\log{1\over t^2}+
4\sum_{k=0}^{+\infty}c_k\,t^{2k}\bigg),\qquad
c_k=c'_k+c''_k,\leqno(17)
$$
qui permettra de calculer un équivalent arbitrairement précis de 
$\varphi(r)$ quand $r\to 1$. Observons que
$$
c_k=c'_k+c''_k={-1/2\choose k}\,\alpha_k
$$
avec
$$
\eqalign{
\alpha_k&=\int_0^1{v^{2k}\,dv\over\sqrt{1{+}v^2}}
+\int_1^{+\infty}\Bigg({v^{2k}\over\sqrt{1{+}v^2}}-\sum_{\ell=0}^k
{-1/2\choose\ell}\,v^{2k-2\ell-1}\Bigg)\,dv
+\sum_{\ell=0}^{k-1}{1\over2(\ell{-}k)}{-1/2\choose\ell}.\cr}
$$
Un calcul direct donne 
$$
c_0=\alpha_0=\int_0^1{dv\over\sqrt{1{+}v^2}}+
\int_1^{+\infty}\bigg({1\over\sqrt{1{+}v^2}}-{1\over v}\bigg)dv=\log 2.
$$
Écrivons d'autre part
$$
{v^{2k}\over\sqrt{1+v^2}}=v^{2k-1}\cdot {v\over\sqrt{1+v^2}}\quad\hbox{et}\quad
(\sqrt{1+v^2})'={v\over\sqrt{1+v^2}},
$$
et intégrons par parties après avoir factorisé $v^{2k-1}$. Il vient
$$
\eqalign{
\alpha_k&=\sum_{\ell=0}^{k-1}{1\over2(\ell-k)}{-1/2\choose\ell}
+\Big[v^{2k-1}\,\sqrt{1+v^2}\,\Big]_0^1-
\int_0^1(2k-1)\,v^{2k-2}\,\sqrt{1+v^2}\,dv\cr
&\qquad{}+\Bigg[v^{2k-1}\Bigg(\sqrt{1+v^2}-\sum_{\ell=0}^k{-1/2\choose\ell}\,
{v^{1-2\ell}\over 1-2\ell}\Bigg)\Bigg]_1^{+\infty}\cr
&\qquad{}-\int_1^{+\infty}(2k-1)\,v^{2k-2}\Bigg(\sqrt{1+v^2}-\sum_{\ell=0}^k
{-1/2\choose\ell}\,{v^{1-2\ell}\over 1-2\ell}\Bigg)\,dv.\cr}
$$
Ceci suggère de calculer $\alpha_k+(2k-1)\alpha_{k-1}$, en utilisant la
simplification
$$
v^{2k-2}\,\sqrt{1+v^2}-{v^{2k-2}\over\sqrt{1+v^2}}={v^{2k}\over\sqrt{1+v^2}}.
$$
Nous en déduisons
$$
\eqalign{
\alpha_k&+(2k-1)\alpha_{k-1}=-(2k-1)\alpha_k+\sum_{\ell=0}^k
{-1/2\choose\ell}\,{1\over 1-2\ell}\cr
&+\int_1^{+\infty}\kern-10pt(2k-1)\,v^{2k-2}\Bigg(\!\sum_{\ell=0}^k
{-1/2\choose\ell}\,
\bigg({v^{1-2\ell}\over 1-2\ell}-v^{1-2\ell}\bigg)-\sum_{\ell=0}^{k-1}
{-1/2\choose\ell}\,
v^{-1-2\ell}\!\Bigg)dv\cr
&+2k\sum_{\ell=0}^{k-1}{1\over2(\ell-k)}{-1/2\choose\ell}
+(2k-1)\sum_{\ell=0}^{k-2}{1\over 2(\ell-(k-1))}{-1/2\choose\ell}.\cr}
$$
Un changement d'indice $\ell=\ell'-1$ dans les sommes associées à $k-1$
permet d'éliminer la quasi totalité des termes. Il reste seulement
le terme $\ell=k$ de la première sommation, d'où la relation de récurrence
$$
2k\,\alpha_k+(2k-1)\alpha_{k-1}=-{-1/2\choose k}\,{1\over 2k-1},\quad
\hbox{i.e.}\quad
{\alpha_k\over{-1/2\choose k}}
-{\alpha_{k-1}\over{-1/2\choose k-1}}=-{1\over 2k(2k-1)}.
$$
Nous obtenons par conséquent
$$
{c_k\over{-1/2\choose k}^2}={\alpha_k\over{-1/2\choose k}}
={\alpha_0\over 1}-\sum_{\ell=1}^k{1\over 2\ell(2\ell-1)}
=\log 2-\sum_{\ell=1}^{2k}{(-1)^{\ell-1}\over \ell}
$$
et l'expression explicite
$$
c_k=w_k\,\Bigg(\log 2-
\sum_{\ell=1}^{2k}{(-1)^{\ell-1}\over \ell}\Bigg).
\leqno(18)
$$
Le reste de la série alternée donnant $\log 2$ est majoré par la moitié du 
dernier terme calculé, soit $1/4k$, donc d'après (8) nous avons
$0<c_k<{1\over \pi^2k^2}$ si $k\geq 1$, de sorte que le rayon de convergence
vaut~$1$. De (9) et (17) nous tirons 
pour $r\to 1-0$ le développement bien connu de l'intégrale elliptique
$$
\varphi(r)={1\over\pi}\Bigg(
\sum_{k=0}^{+\infty}w_kt^{2k}\log{1\over t^2}
+4\sum_{k=0}^{+\infty}c_kt^{2k}\Bigg),\qquad t^2=1-r^2,\leqno(19)
$$
avec
$$
w_0=1,~~w_1={1\over 4},~~w_2={9\over 64},~~
c_0=\log 2,~~c_1={1\over 4}\bigg(\log 2-{1\over 2}\bigg),
~~c_2={9\over 64}\bigg(\log 2-{7\over 12}\bigg).
$$
Explicitons le début du développement limité en $r=1$ en posant $r=1+h$.
Pour $r=1+h<1$ ($h<0$), nous avons $t^2=1-r^2=-2h-h^2=2|h|(1+h/2)$, d'où
$$
\eqalign{
&\log{1\over t^2}=\log{1\over 2|h|(1+h/2)}=\log{1\over |h|}-\log 2-{1\over 2}h
+{1\over 8}h^2+O(h^2),\cr
&\sum_{k=0}^{+\infty}w_kt^{2k}=1+{1\over 4}(-2h-h^2)+{9\over 64}(2h)^2+O(h^3),\cr
&4\sum_{k=0}^{+\infty}c_kt^{2k}=4\log 2+
\bigg(\log 2 -{1\over 2}\bigg)(-2h-h^2)+
{9\over 16}\bigg(\log 2 -{7\over 12}\bigg)(2h)^2+O(h^3),\cr}
$$
et
$$
\eqalign{
\varphi(1+h)&={1\over\pi}\Bigg(
\bigg(1-{1\over 2}h+{5\over 16}h^2+O(h^3)
\bigg)\bigg(\log{1\over |h|}-\log 2-{1\over 2}h+
{1\over 8}h^2+O(h^3)\bigg)\cr
&\kern67pt{}+4\log 2-\big(2\log 2-1\big)h+
\bigg({5\over 4}\log 2-{13\over 16}\bigg)h^2+O(h^3)\Bigg),\cr}
$$
soit encore (par ordre de grandeur décroissant des termes)~:
$$
\leqalignno{
\varphi(1+h)&={1\over\pi}\Bigg(
\log{1\over |h|}+3\log 2-{1\over 2}h\log{1\over |h|}-
\bigg({3\over 2}\log 2-{1\over 2}\bigg)h\cr
&\kern67pt{}+{5\over 16}h^2\log{1\over |h|}+
\bigg({15\over 16}\log 2-{7\over 16}\bigg)h^2
+O\bigg(h^3\log{1\over|h|}\bigg)\Bigg).&(20)\cr}
$$
Pour $r=1+h>1$, la relation $\varphi(r)={1\over r}\varphi({1\over r})$ donne
de même
$$
\varphi(r)=
{1\over 1+h}\Bigg({1\over\pi}\sum_{k=0}^{+\infty}w_kt^{2k}\log{1\over t^2}
+\sum_{k=0}^{+\infty}c_kt^{2k}\Bigg),\quad t^2=1-{1\over r^2}=2h-3h^2+O(h^3),
$$
et après simplifications on voit que le développement (20) est encore 
valable pour $h>0$. Un passage à la limite $r\to 0$, $t\to 1-0$ dans (19) 
fournit en particulier la relation $\sum_{k\geq 0}c_k={\pi\over 4}$. 
Le lemme suivant nous sera utile.
\medskip

{\statement Lemme 1.} {\it Pour $h>0$, l'écart
$$
\leqalignno{
\rho(h)&=\varphi(1+h)-{1\over\pi}\Bigg(\log{1\over h}+3\log 2
-{1\over 2}h\,\log{1\over h}-\bigg({3\over 2}\log 2-{1\over 2}\bigg)h\Bigg)
&(21)\cr
&=\varphi(1+h)-{1\over 2\pi}\bigg((h-2)\log{h\over 8}+h\bigg)&(21')\cr}
$$
admet la majoration
$$
|\rho(h)|\leq h^2\bigg(2+\log\bigg(1+{1\over h}\bigg)\bigg).\leqno(22)
$$}%

{\it Preuve.} L'utilisation de la formule de Taylor-Lagrange donne
$(1+h)^{-1}=1-h+\theta_1h^2$, $t^2=1-{1\over r^2}=2h-3\theta_2h^2$, 
$\theta_i\in{}]0,1[$, ainsi que $t^2\leq 2h$ et
$$
\eqalign{
\log{1\over t^2}&=\log {r^2\over (r-1)(r+1)}=\log{1\over h}+2\log(1+h)
-\log\bigg(1+{h\over 2}\bigg)-\log 2\cr
&=\log{1\over h}-\log 2+{3\over 2}h-{7\over 8}\theta_3h^2,\quad\theta_2\in
{}]0,1[,\cr}
$$
tandis que les restes $\sum_{k\ge 2}w_kt^{2k}$ et $\sum_{k\ge 2}c_kt^{2k}$
sont majorés respectivement par 
$$
{w_2t^4\over 1-t^2}\leq 4w_2r^2h^2\leq{225\over 256}h^2\quad\hbox{et}\quad
{c_2t^4\over 1-t^2}\leq 4c_2r^2h^2<{1\over 10}h^2\quad\hbox{si}~~
h\leq{1\over 4},~r=1+h\leq{5\over 4}.
$$
Pour $h\leq{1\over 4}$, on a donc une égalité
$$
\eqalign{
\varphi(1+h)&={1\over\pi}(1-h+\theta_1h^2)\times\Bigg(\cr
&\qquad{}
\bigg(1+{1\over 4}(2h-3\theta_2h^2)+{225\over 256}\theta_4h^2\bigg)
\bigg(\log{1\over |h|}-\log 2+{3\over 2}h-{7\over 8}\theta_3h^2\bigg)\cr
&\qquad{}+4\log 2+\bigg(\log 2-{1\over 2}\bigg)(2h-3\theta_2h^2)
+{4\over 10}\theta_5h^2\Bigg)\cr}
$$
avec $\theta_i\in{}]0,1[\,$. Un développement de $\rho(h)$ à l'aide de 
cette expression, suivi d'une\break majoration de chaque terme par sa valeur 
absolue, montre que pour $h\leq{1\over 4}$ on a\break
$|\rho(h)|\leq h^2(0.885\,\log{1\over h}+2.11)$, 
de sorte que (22) est satisfait. Lorsque $h\geq{1\over 4}$,
on~a
$$
\rho'(h)=\varphi'(1+h)-{1\over2\pi}\bigg(\log{h\over 8}+2-{2\over h}\bigg),
\qquad
\varphi'(r)=-\sum_{k=0}^{+\infty}(2k+1)w_k\,r^{-2k-2},
$$
et d'après (8)
$$
\sum_{k=0}^{+\infty}{2\over\pi}r^{-2k-2}<-\varphi'(r)<{1\over r^2}+\sum_{k=1}^{+\infty}{3k\over\pi k}r^{-2k-2}
<\sum_{k=0}^{+\infty}r^{-2k-2}={1\over r^2-1},
$$
d'où
$$
\eqalign{
&{2\over\pi}\,{1\over h(h+2)}<-\varphi'(1+h)<{1\over h(h+2)},\cr
&{1\over2\pi}\bigg(\log{8\over h}-2+{2\over h}-{2\pi\over h(h+2)}\bigg)<
\rho'(h)<{1\over2\pi}\bigg(\log{8\over h}-2+{2\over h+2}\bigg).\cr}
$$
On obtient ainsi
$$
\eqalign{
&-1.72<{1\over 2\pi}\bigg(\log 4{-}2{+}{1\over 4}{-}{32\pi\over 9}\bigg)<
\rho'(h)<{1\over 2\pi}\bigg(\log 32{-}2{+}{8\over 9}\bigg)<1.51
~~\hbox{sur}~\bigg[{1\over 4},2\bigg],\cr
&-{1\over 2\pi}\bigg(\log{h\over 8}+2\bigg)<
\rho'(h)<{1\over 2\pi}\bigg(\log 4 -{3\over 2}\bigg)<0\qquad\hbox{
sur~~$[2,+\infty[\,$},\cr}
$$
par suite $|\rho'(h)|\leq{1\over 2\pi}(h-1-\log 8+2)\leq{1\over 2\pi}h$
pour $h\in[2,+\infty[$. Comme $\rho(2)\simeq 0.00249<{1\over\pi}$, on voit
que $|\rho(h)|\leq {1\over 4\pi}h^2$, ce qui montre que (22) est encore 
vérifié sur $[2,+\infty[\,$. Un calcul numérique de $\rho(h)$ en des points 
suffisamment rapprochés de $[{1\over 4},2]$ permet enfin de vérifier (22)
sur ce dernier intervalle.\qed
\medskip

On décompose maintenant l'intégrale (14) sur les intervalles $[0,1]$ et 
$[1,+\infty[\,$, en commençant par l'intégrale de~$\varphi$ sur 
$[1,+\infty[\,$. Le changement de variable $r=1+t/4x$ fournit
$$
\int_1^{+\infty}e^{-4xr}\,\varphi(r)\,dr = {e^{-4x}\over 4x}
\int_0^{+\infty}e^{-t}\,\varphi\Big(1+{t\over 4x}\Big)\,dt,\leqno(23)
$$
et le lemme 1~$(21')$ donne de cette dernière intégrale une approximation
$$
\eqalign{
{e^{-4x}\over 8\pi x}&
\int_0^{+\infty}e^{-t}\bigg(\Big({t\over 4x}-2\Big)
\log{t\over 32x}+{t\over 4x}\bigg)\,dt\cr
&={e^{-4x}\over 8\pi x}\bigg(
\log(32x)\Big(2-{1\over 4x}\Big)+2\gamma
+{1\over 4x}\int_0^{+\infty}e^{-t}(t\log t+t)\,dt\bigg)\cr
&={e^{-4x}\over 4\pi x}\bigg(\log x+\gamma+5\log 2-
{\log x\over 8x}-{\gamma+5\log 2-2\over 8x}\bigg),\cr}
$$
avec une erreur majorée par
$$
\eqalign{
{e^{-4x}\over 4x}\int_0^{+\infty}e^{-t}\bigg({t\over 4x}\bigg)^2&\bigg(
2+\log\bigg(1+{4x\over t}\bigg)\bigg)dt\cr
&={e^{-4x}\over 4x}\bigg({1\over 4x^2}+{1\over 16x^2}
\int_0^{+\infty}t^2\,e^{-t}\,\log{t+4x\over t}\,dt\bigg).\cr}
$$
Si on écrit
$$
0<\log{t+4x\over t}=\log{4x\over t} +
\log\bigg(1+{t\over 4x}\bigg)
\leq \log{4x\over t} + {t\over 4x},
$$
on voit que
$$
\eqalign{
\int_0^{+\infty}t^2\,e^{-t}\,\log{t+4x\over t}\,dt&\leq
\int_0^{+\infty}t^2\,e^{-t}\,\bigg(\log{4x\over t}+{t\over 4x}\bigg)\,dt
=2\log 4x+{3\over 2x}+2\gamma-3.\cr}
$$
On en déduit 
$$
\int_1^{+\infty}e^{-4xr}\,\varphi(r)\,dr =
{e^{-4x}\over 4\pi x}\bigg(\log x+\gamma+5\log 2-
{\log x\over 8x}\bigg)+{e^{-4x}\over 4x}\,R_1(x),
\leqno(24)
$$
et au moyen d'un calcul numérique il vient
$$
|R_1(x)|<{\gamma+5\log 2-2\over 8\pi x}+{1\over 4x^2}+
{2\log 4x+{3\over 2x}+2\gamma-3\over 16x^2}<{0.483\over x}.\leqno(25)
$$
Il faut maintenant estimer les deux intégrales
$\int_0^1e^{-4xr}\,\sum_{k\geq 2x+1}w_k\,r^{2k}\,dr$
et\break $\int_1^{+\infty}e^{-4xr}\,\sum_{k\leq 2x}w_k\,r^{2k}\,dr$.
Des intégrations par parties successives donnent
$$
\leqalignno{
\int_0^1e^{-4xr}\,r^{2k}\,dr&=e^{-4x}\sum_{\ell=1}^{+\infty}
{(4x)^{\ell-1}\over (2k+1)\cdots(2k+\ell)},&(26)\cr
\int_1^{+\infty}e^{-4xr}\,r^{2k}\,dr&=
{e^{-4x}\over 4x}\bigg(1+
\sum_{\ell=1}^{2k}{2k(2k-1)\cdots(2k-\ell+1)\over(4x)^\ell}\bigg).&(27)\cr}
$$
En combinant les identités (14), (15), (24), (26), (27) il vient
$$\quad
\Delta(x)={e^{-4x}\over 4x}\Bigg({1\over\pi}\Big(\log x+\gamma+5\log 2\Big)
-{\log x\over 8\pi x}-\sum_{k=0}^{2x}w_k+S(x)+R_1(x)+R_2(x)\Bigg)\leqno(28)
$$
avec
$$
~~S(x)=\sum_{k=2x+1}^{+\infty}
\sum_{\ell=1}^{2x-1}{w_k\,(4x)^\ell\over (2k+1)\cdots(2k+\ell)}
-\sum_{k=1}^{2x}
\sum_{\ell=1}^{2x-1}w_k\,{2k(2k-1)\cdots(2k-\ell+1)
\over(4x)^\ell},\leqno(29)
$$
et
$$
~~R_2(x)=
\sum_{k=2x+1}^{+\infty}\sum_{\ell=2x}^{+\infty}
{w_k\,(4x)^\ell\over (2k+1)\cdots(2k+\ell)}
-\sum_{k=1}^{2x}\sum_{\ell=2x}^{+\infty}
w_k\,{2k(2k-1)\cdots(2k-\ell+1)\over(4x)^\ell}
\leqno(30)
$$
(dans la dernière sommation, seuls les termes $\ell\leq 2k$ sont non nuls).
La formule (28) nous amène maintenant à étudier le développement asymptotique
de $\sum_{k=0}^{2x}w_k$. Ce développement est facile à établir à partir de (19), 
on pourrait même le calculer à un ordre arbitrairement grand.
\medskip

{\statement Lemme 2.} {\it On a
$$
\leqalignno{
\qquad&w_k={1\over\pi k}\bigg(1-{1\over 2(2k-1)}+\varepsilon_k\bigg)
\quad\hbox{où}~~
{1\over 12k(2k-1)}<\varepsilon_k<{5\over 16k(2k-1)},~~k\geq 1,
&(31)\cr
&\sum_{k=0}^{2x}w_k={1\over\pi}\Big(\log x+5\log 2+\gamma\Big)
+R_3(x),\qquad {1\over 4\pi x}<R_3(x)<{19\over 48\pi x}.
&(31')\cr}
$$}%

{\it Preuve.} La minoration (31) va résulter de la formule
d'Euler-Maclaurin du \S1 appliquée à la fonction $f(x)=\log{2x-1\over 2x}$.
Celle-ci fournit
$$
{1\over 2}\,\log w_k=\sum_{i=1}^kf(i)=C+\int_1^kf(x)\,dx+{1\over 2}f(k)+
\sum_{j=1}^p{b_{2j}\over (2j)!}\,f^{(2j-1)}(k)+\tilde R_p
$$
où $C$ est une constante et où le reste $\tilde R_p$ est le produit 
par un facteur dans $[0,1]$ du terme suivant, à savoir
$$
{b_{2p+2}\over (2p+2)!}\,f^{(2p+1)}(k)={2^{2p+1}\,b_{2p+2}\over(2p+1)(2p+2)}
\bigg({1\over (2k-1)^{2p+1}}-{1\over (2k)^{2p+1}}\bigg).
$$
On a ici 
$$
\eqalign{
\int_1^kf(x)\,dx&={1\over 2}(2k-1)\log(2k-1)-k\log k-(k-1)\,\log 2\cr
&=\bigg(k-{1\over 2}\bigg)\log\bigg(1-{1\over 2k}\bigg)
-{1\over 2}\,\log k+{1\over 2}\,\log 2\cr}
$$
et la constante $C$ est ajustée par la formule de Wallis, d'où,
avec $b_2={1\over 6}$,
$$
\eqalign{
\log w_k&=\log{1\over\pi k} + 2k\,\log\bigg(1-{1\over 2k}\bigg)+1
+2\theta\, b_2\bigg({1\over (2k-1)}-{1\over 2k}\bigg)\cr
&\geq \log{1\over\pi k}-{1\over 4k}-\sum_{\ell=3}^{+\infty}
{1\over \ell (2k)^{\ell-1}}
> \log{1\over\pi k}-{1\over 4k}-{1\over 3}\,{1\over (2k)^2}\,
{1\over 1-{1\over 2k}}.\cr}
$$
L'inégalité $e^{-x}\geq 1-x$ donne alors
$$
w_k>{1\over\pi k}\bigg(1-{1\over 4k}-{1\over 6k(2k-1)}\bigg)
={1\over\pi k}\bigg(1-{1\over 2(2k-1)}+{1\over 12k(2k-1)}\bigg)
$$
et la minoration (31) s'ensuit pour tout $k\geq 1$. Dans l'autre sens,
nous obtenons
$$
\log w_k<\log{1\over\pi k} - {1\over 4k }- {1\over 12k^2} - {1\over 32k^3}
+ {1\over 6k(2k-1)}
=\log{1\over\pi k} - {1\over 4k} + {1\over 12k^2(2k-1)} - {1\over 32k^3}
$$
et l'inégalité $e^{-x}\leq 1-x+{1\over 2}x^2$ implique
$$
\leqalignno{
w_k&<{1\over\pi k}\Bigg(1- \bigg({1\over 4k} - {1\over 12k^2(2k-1)} + {1\over 32k^3}\bigg) + {1\over 2}\bigg({1\over 4k}\bigg)^2\Bigg)\cr
&&\hbox{d'où (par différence polynomiale et réduction au même dénominateur)}\cr
w_k&<{1\over\pi k}\Bigg(1- {1\over 2(2k-1)}+{5\over 16k(2k-1)}\Bigg)\qquad
\hbox{si $k\geq 3$.}\cr}
$$
On peut vérifier que l'inégalité finale a bien lieu aussi si $k=1,2$, ce qui
entraîne l'estimation (31). La formule (19) implique d'autre part
$$
\eqalign{
w_0+\sum_{k=1}^{+\infty}\Big(w_k-{1\over \pi k}\Big)r^{2k}
&=\varphi(r)-{1\over\pi}\,\log{1\over 1-r^2}\cr
&={1\over\pi}\,\big(\varphi(t)-1\big)\,\log{1\over 1-r^2}
+{4\over\pi}\,\log2+\sum_{k\geq 1}
c_k\,t^{2k}\cr}
$$
avec $t=\sqrt{1-r^2}$ et $\varphi(t)=1+O(1-r^2)$. Un passage à la limite
quand $r\to 1-0$ et $t\to 0$ donne donc
$$
w_0+\sum_{k=1}^{+\infty}\Big(w_k-{1\over \pi k}\Big)={4\over\pi}\,\log2.
$$
On en déduit
$$
w_0+\sum_{k=1}^{2x}\Big(w_k-{1\over \pi k}\Big)-{4\over\pi}\,\log2 
=\sum_{2x+1}^{+\infty}\Big({1\over \pi k}-w_k\Big)
$$
et l'encadrement (31) donne
$$
0<\sum_{2x+1}^{+\infty}\Big({1\over \pi k}-w_k\Big)\leq
\sum_{2x+1}^{+\infty}{1\over 2\pi\,k(2k-1)}
<\sum_{2x+1}^{+\infty}{1\over 4\pi}\,{1\over k(k-1)}
={1\over 8\pi x}.
$$
L'estimation d'Euler-Maclaurin
$$
\sum_{k=1}^{2x}{1\over k}=\log(2x)+\gamma+{1\over 4x}+{b_2\over 2(2x)^2}
- {b_4\over4(2x)^4}+\cdots\leqno(32)
$$
implique alors $(31')$.\qed

Il reste à évaluer la somme $S(x)$, ce qui est beaucoup plus difficile 
du fait d'une compensation partielle entre les termes positifs et négatifs. 
L'approximation (31) du lemme~2 implique
$$
S(x)={2\over\pi}\bigg(T(x)-{1\over 2}U(x)+{5\over 8}R_4(x)\bigg),\leqno(33)
$$
où, en convenant qu'un produit vide vaut~$1$ et que
$(2k{-}2)\cdots(2k{-}\ell{+}1)={1\over 2k-1}$ si $\ell=1$,
$$
\leqalignno{
\qquad T(x)&=\sum_{\ell=1}^{2x-1}
\sum_{k=2x+1}^{+\infty}{(4x)^\ell\over 2k(2k+1)\cdots(2k+\ell)}
-\sum_{\ell=1}^{2x-1}
\sum_{k=1}^{2x}{(2k-1)\cdots(2k-\ell+1)\over(4x)^\ell},&(34)\cr
\qquad U(x)&=\sum_{\ell=1}^{2x-1}
\sum_{k=2x+1}^{+\infty}{(4x)^\ell\over (2k-1)\cdots(2k+\ell)}
-\sum_{\ell=1}^{2x-1}
\sum_{k=1}^{2x}{(2k-2)\cdots(2k-\ell+1)\over (4x)^\ell},&(35)\cr}
$$
et où la nouvelle erreur $R_4(x)$ admet la majoration
$$
~~|R_4(x)|\leq
\sum_{\ell=1}^{2x-1}
\sum_{k=2x+1}^{+\infty}{(4x)^\ell/2k\over (2k-1)\cdots(2k+\ell)}
+\sum_{\ell=1}^{2x-1}
\sum_{k=1}^{2x}{(2k-2)\cdots(2k-\ell+1)\over 2k\,(4x)^\ell}.\leqno(36)
$$
La~méthode consiste à sommer d'abord par rapport à l'indice~$k$, et à utiliser
pour cela des ``intégrations par parties discrètes''. Si on pose 
$$
u^{a,b}_k:={1\over (2k+a)(2k+a+1)\cdots(2k+b-1)},\qquad a\leq b\leqno(37)
$$
(en convenant que le dénominateur vaut~$1$ si $a=b$), alors
$$
\eqalign{
u^{a,b}_k-u^{a,b}_{k+1}&= {(2k+b)(2k+b+1)-(2k+a)(2k+a+1)
\over (2k+a)(2k+a+1)\cdots(2k+b+1)}\cr
&={(b-a)(4k+a+b+1)\over (2k+a)(2k+a+1)\cdots(2k+b+1)}.\cr}
$$
L'encadrement $2(2k+a)\leq 4k+a+b+1\leq 2(2k+b+1)$ implique
$$
{1\over (2k+a+1)\cdots(2k+b+1)}\leq {u^{a,b}_k-u^{a,b}_{k+1}\over 2(b-a)}
\leq{1\over (2k+a)(2k+a+1)\cdots(2k+b)}
$$
avec une erreur par excès et une erreur par défaut toutes deux égales à
$$
{b-a+1\over 2}\,{1\over (2k+a)(2k+a+1)\cdots(2k+b+1)}.
$$
Ces inégalités donnent en particulier par sommation de
$\sum_{k=2x+1}^{+\infty}{u^{a-1,b-1}_k-u^{a-1,b-1}_{k+1}\over 2(b-a)}$\break
l'inégalité
$$
\sum_{k=2x+1}^{+\infty}
{1\over (2k+a)\cdots(2k+b)}\leq 
{u^{a-1,b-1}_{2x+1}\over 2(b-a)}
={1\over 2(b-a)}\,{1\over (4x+a+1)\cdots(4x+b)},
$$
avec une erreur par excès égale à
$$
\eqalign{{b-a+1\over 2}
\sum_{k=2x+1}^{+\infty}
{1\over (2k+a-1)\cdots(2k+b)}
&\leq{1\over 4}\,{1\over (4x+a)\cdots(4x+b)}\cr}
$$
et une ``erreur sur l'erreur'' (de nouveau par excès) égale à
$$
{(b{-}a{+}1)(b{-}a{+}2)\over 4}
\sum_{k=2x+1}^{+\infty}
{1\over (2k+a-2)\cdots(2k+b)}
\leq{b{-}a{+}1\over 8}\,{1\over (4x+a-1)\cdots(4x+b)},
$$
autrement dit
$$
\leqalignno{
\sum_{k=2x+1}^{+\infty}{1\over (2k+a)\cdots(2k+b)}&=
{1\over 2(b-a)}{1\over (4x+a+1)\cdots(4x+b)}
-{1\over 4}{1\over (4x+a)\cdots(4x+b)}\cr
&\qquad\quad
{}+\theta\,{b-a+1\over 8}\,{1\over (4x+a-1)\cdots(4x+b)},\qquad
\theta\in[0,1].&(38^{a,b}_3)\cr}
$$
On pourrait bien sûr pousser ce développement à 3 termes à un
nombre quelconque $p$ de termes si nécessaire, ce que nous noterons
$(37^{a,b}_p)$, et on utilisera ici les cas $p=2,3$. Pour les sommes 
$\sum_{k=1}^{2x}\ldots~$, on pose de manière similaire
$$
v^{a,b}_k=(2k-a)(2k-a-1)\cdots(2k-b+1),\qquad a\leq b,\leqno(39)
$$
et on obtient
$$
\eqalign{
v^{a,b}_k-v^{a,b}_{k-1}
&=(2k-a-2)\cdots(2k-b+1)\big((2k-a)(2k-a-1)-(2k-b)(2k-b-1)\big)\cr
&=(2k-a-2)\cdots(2k-b+1)\big((b-a)(4k-a-b-1)\big).\cr}
$$
Pour $a<b$, l'encadrement $2(2k-b)\leq(4k-a-b-1)\leq 2(2k-a-1)$ implique
$$
(2k-a-2)\cdots(2k-b)\leq
{v^{a,b}_k-v^{a,b}_{k-1}\over 2(b-a)}\leq (2k-a-1)\cdots(2k-b+1)
$$
avec une erreur par excès et une erreur par défaut toutes deux égales à
$$
{1\over 2}(b-a-1)\,(2k-a-2)\cdots(2k-b+1).
$$
Par sommation de 
$\sum_{k=1}^{2x}{v^{a,b}_k-v^{a,b}_{k-1}\over 2(b-a)}$, on obtient donc
$$
\sum_{k=1}^{2x}(2k-a-1)\cdots(2k-b+1)\geq {v^{a,b}_{2x}-
v^{a,b}_0\over 2(b-a)}
$$
avec une erreur par défaut
$$
{b-a-1\over 2}
\sum_{k=1}^{2x}(2k-a-2)\cdots(2k-b+1)\leq{v^{a,b-1}_{2x}-v^{a,b-1}_0\over 4}
$$
et une erreur (par excès) sur l'erreur égale à
$$
{(b-a-1)(b-a-2)\over 4}
\sum_{k=1}^{2x}(2k-a-2)\cdots(2k-b+2)
\leq{b-a-1\over 8}\,\Big(v^{a,b-2}_{2x}-v^{a,b-2}_0\Big),
$$
autrement dit, il existe $\theta\in[0,1]$ tel que
$$
\leqalignno{
\sum_{k=1}^{2x}(2k-&a-1)\cdots(2k-b+1)\cr
\noalign{\vskip-8pt}\cr
&={1\over 2(b-a)}\Big(v_{2x}^{a,b}-v^{a,b}_0\Big)
+{1\over 4}\Big(v^{a,b-1}_{2x}-v^{a,b-1}_0\Big)
-\theta\,{b-a-1\over 8}\,\Big(v^{a,b-2}_{2x}-v^{a,b-2}_0\Big),\cr
&={1\over 2(b-a)}\,v_{2x}^{a,b}
+{1\over 4}\,v^{a,b-1}_{2x}-\theta\,{b-a-1\over 8}\,v^{a,b-2}_{2x}
+C^{a,b}_3,&(40^{a,b}_3)\cr}
$$
avec
$$
|C^{a,b}_3|\leq {1\over 2(b-a)}|v^{a,b}_0|+{1\over 4}|v^{a,b-1}_0|+
{b-a-1\over 8}|v^{a,b-2}_0|,\leqno(41^{a,b}_3)
$$
et en particulier $C^{a,b}_3=0$ si $a=0$. La version à l'ordre 2
avec erreur initiale par excès donne plus simplement
$$
\leqalignno{
\sum_{k=1}^{2x}(2k-&a-2)\cdots(2k-b)
={1\over 2(b-a)}\Big(v_{2x}^{a,b}-v^{a,b}_0\Big)
-\theta\,{1\over 4}\Big(v^{a,b-1}_{2x}-v^{a,b-1}_0\Big)\cr
&={1\over 2(b-a)}\,(4x-a)\cdots(4x-b+1)
-\theta\,{1\over 4}\,(4x-a)\cdots(4x-b+2)+C^{a,b}_2.&(40^{a,b}_2)\cr}
$$
Pour le développement à l'ordre 3, il sera commode d'utiliser une 
transformation supplémentaire
$$v^{a,b}_k-v^{a+1,b+1}_k=(2k-a-1)\cdots(2k-b+1)\Big((2k-a)-(2k-b)\Big)=
(b-a)v^{a+1,b}_k,
$$
et en appliquant ceci pour $(a,b)$, $(a,b-1)$ et $k=2x$ on voit
que le développement $(40^{a,b}_3)$ peut s'écrire sous la forme équivalente
$$
\leqalignno{
&\sum_{k=1}^{2x}(2k-a-1)\cdots(2k-b+1)-C^{a,b}_3\cr
\noalign{\vskip-8pt}\cr
&={1\over 2(b-a)}v_{2x}^{a+1,b+1}+{3\over 4}v^{a+1,b}_{2x}
+{b-a-1\over 8}\Big(2\,v^{a+1,b-1}_{2x}
-\theta\,v^{a,b-2}_{2x}\Big),\cr
&={1\over 2(b-a)}\,(4x-a-1)\cdots(4x-b)
+{3\over 4}\,(4x-a-1)\cdots(4x-b+1)\cr
&\kern1.1cm{}+{b-a-1\over 8} \Big(2(4x-a-1)\cdots(4x-b+2)
-\theta\,(4x-a)\cdots(4x-b+3)\Big)&(\widetilde{40}^{\,a,b}_3)\cr}
$$
D'après (34), $(38^{0,\ell}_3)$, $(\widetilde{40}^{\,0,\ell}_3)$, nous 
obtenons 
$$
T(x)=T'(x)-T''(x)+R_5(x)\leqno(42)
$$
avec
$$
\leqalignno{
T'(x)&=
\sum_{\ell=1}^{2x-1}
{1\over 2\ell}\bigg({(4x)^\ell\over (4x+1)\cdots(4x+\ell)}-
{(4x-1)\cdots(4x-\ell)\over(4x)^\ell}\bigg),&(43)\cr
T''(x)&=
\sum_{\ell=1}^{2x-1}{1\over 4}\,{(4x)^\ell\over 4x(4x+1)\cdots(4x+\ell)}+
{3\over 4}\,{(4x-1)\cdots(4x-\ell+1)\over(4x)^\ell},&(43')\cr
\qquad |R_5(x)|&\leq{1\over 8}
\sum_{\ell=1}^{2x-1}\bigg({(\ell+1)\,(4x)^\ell\over (4x-1)4x\cdots(4x+\ell)}+
{2(\ell-1)\,(4x-1)\cdots(4x-\ell+2)\over(4x)^\ell}\bigg).&(43'')\cr}
$$
Le dernier terme de la dernière ligne provient de 
$(\widetilde{40}^{\,0,\ell}_3)$
si on observe que l'inégalité $4x\leq 2(4x-\ell+2)$ pour $\ell\leq 2x-1$
entraîne
$$
4x(4x-1) \cdots(4x-\ell+3)\leq 2(4x-1)\cdots(4x-\ell+2).
$$
De manière analogue, grâce à $(35)$, $(38^{-1,\ell}_2)$ et $(40^{0,\ell-1}_2)$ 
on obtient une décomposition 
$$
U(x)=U'(x)-U''(x)+R_6(x)\leqno(44)
$$
avec
$$
\leqalignno{
\kern36pt U'(x)&=\sum_{\ell=1}^{2x-1}
{1\over 2(\ell+1)}\,{(4x)^\ell\over 4x\cdots(4x+\ell)}-
\sum_{\ell=2}^{2x-1}{1\over 2(\ell-1)}\,
{4x(4x-1)\cdots(4x-\ell+2)\over(4x)^\ell},&(45)\cr
U''(x)&={1\over 4x}\sum_{k=1}^{2x}{1\over 2k-1}\qquad\hbox{(terme négatif
$\ell=1$ issu de $U(x)$)},&(45')\cr
\quad|R_6(x)|&\leq
{1\over 4}\sum_{\ell=1}^{2x-1}
{(4x)^\ell\over (4x-1)\cdots(4x+\ell)}+
{1\over 4}\sum_{\ell=2}^{2x-1}{4x(4x-1)\cdots(4x-\ell+3)\over(4x)^\ell}.
&(45'')\cr}
$$
Les restes $R_2(x)$ $[\,$resp.\ $R_4(x)\,$] se majorent de même à l'aide de
$(38^{0,\ell}_2)$ et $(40^{-1,\ell-1}_2)$ $[\,$resp.\ 
$(38^{-1,\ell}_2)$ et $(40^{0,\ell-2}_2)\,$]
et (8), (30), (36) entraînent ainsi
$$
\leqalignno{
|R_2(x)|&\leq{2\over \pi}\Bigg(
\sum_{\ell=2x}^{+\infty}\sum_{k=2x+1}^{+\infty}
{(4x)^\ell\over (2k)\cdots(2k+\ell)}
+\sum_{\ell=2x}^{+\infty}\sum_{k=1}^{2x}
{(2k-1)\cdots(2k-\ell+1)\over(4x)^\ell}\Bigg)\cr
&\leq{2\over \pi}\sum_{\ell=2x}^{+\infty}{1\over 2\ell}\bigg(
{(4x)^\ell\over (4x+1)\cdots(4x+\ell)}+
{(4x+1)\cdots(4x-\ell+2)\over(4x)^\ell}\bigg),&(46)\cr
|R_4(x)|&\leq
\sum_{\ell=1}^{2x-1}
\sum_{k=2x+1}^{+\infty}{(4x)^{\ell-1}\over (2k-1)\cdots(2k+\ell)}
+\sum_{\ell=1}^{2x-1}
\sum_{k=1}^{2x}{(2k-2)\cdots(2k-\ell+2)\over (4x)^\ell}\cr
&\leq
\sum_{\ell=1}^{2x-1}{1\over 2(\ell+1)}\,{(4x)^{\ell-1}\over
4x\cdots(4x+\ell)}+\sum_{\ell=3}^{2x-1}{1\over 2(\ell-2)}
{4x(4x-1)\cdots(4x-\ell+3)\over(4x)^\ell}&(46')\cr
&\qquad{}+\sum_{k=1}^{2x}{1\over 2k(2k-1)}\,{1\over 4x}
+\sum_{k=1}^{2x}{1\over (2k-1)}\,{1\over (4x)^2}\quad
\hbox{[termes $\ell=1,2$ de la somme]}.&(46'')\cr}
$$
En conclusion, d'après (28), $(31')$, (33) et (42), (44) nous obtenons 
la décomposition
$$
\leqalignno{
\Delta(x)={e^{-4x}\over 4\pi x}\Bigg(
2\,T'(x)&-2\,T''(x)-U'(x)+U''(x)-{\log x\over 8x}\cr
&+\pi\Big(R_1(x)+R_2(x)-R_3(x)\Big)-{5\over 4}R_4(x)+2\,R_5(x)-R_6(x)\Bigg).
&(47)\cr}
$$
{\statement Lemme 3.} {\it On a les inégalités
$$
\leqalignno{
&\log2-{1\over 8x}<
\sum_{k=1}^{2x}{1\over 2k(2k-1)}<\log2-{1\over 2(4x+1)},\kern90pt&(48)\cr
&\sum_{k=1}^{2x}{1\over 2k-1}
<{3\over 2}\log 2+{1\over 2}\Big(\log x+\gamma\Big)+{1\over 24x^2},&(48')\cr
&U''(x)={\log x\over 8x}+R_7(x),\qquad
0<R_7(x)<{1.37\over x}.&(48'')\cr}
$$}

{\it Preuve.} Pour (48), la somme de la série vaut $\log 2$ et on observe que 
le reste admet l'encadrement
$$\eqalign{
{1\over 2(4x+1)}&=\sum_{k=2x+1}^{+\infty}{1\over 4}
\bigg({1\over k-1/2}-{1\over k+1/2}\bigg)\cr
&<\sum_{k=2x+1}^{+\infty}{1\over 2k(2k-1)}<\sum_{k=2x+1}{1\over 4}
\bigg({1\over k-1}-{1\over k}\bigg)={1\over 8x}.\cr}
$$
D'après l'estimation d'Euler-Maclaurin (32), on en déduit d'une part
$$
\eqalign{
\sum_{k=1}^{2x}{1\over 2k-1}&=\sum_{\ell=1}^{4x}{(-1)^{\ell-1}\over \ell}
+\sum_{\ell=1}^{2x}{1\over 2\ell}
=\sum_{k=1}^{2x}{1\over 2k(2k-1)}+{1\over 2}\sum_{k=1}^{2x}{1\over k}\cr
&<\log 2-{1\over 2(4x+1)}+
{1\over 2}\bigg(\log(2x)+\gamma+{1\over 4x}+{1\over 12(2x)^2}\bigg)\cr
&={3\over 2}\log 2+{1\over 2}\Big(\log x+\gamma\Big)+{1\over 8x(4x+1)}+
{1\over 96x^2},\cr}
$$
ce qui donne $(48')$, et d'autre part
$$
\eqalign{
\sum_{k=1}^{2x}{1\over 2k-1}
&>\log 2+{1\over 2}\bigg(\log(2x)+\gamma+{1\over 12(2x)^2}-
{1\over 120(2x)^4}\bigg)\cr
&>{3\over 2}\log 2+{1\over 2}\Big(\log x+\gamma\Big)+{1\over 96x^2}
-{1\over 1920x^4}.\cr}
$$
Une estimation numérique élémentaire 
${3\over 2}\log 2+{1\over 2}\gamma+{1\over 24}<1.37$
fournit alors $(48'')$.\qed
\medskip

On va maintenant vérifier que tous les restes $R_i(x)$ sont d'ordre 
inférieur à celui des termes principaux et plus précisément $O(1/x)$. 
Le plus facile à estimer est $R_6(x)$. On~peut
en effet utiliser la majoration brutale
$$
|R_6(x)|\leq {1\over 4}\sum_{\ell=1}^{2x-1}{1\over 4x(4x-1)}+
{1\over 4}\sum_{\ell=2}^{2x-1}{1\over (4x)^2}\leq 
{1\over 4}\,{2x-1\over 4x(4x-1)}+
{1\over 4}\,{2x-2\over (4x)^2}<{1\over 16x}.\leqno(49)
$$
Passons maintenant à $R_4(x)$, et utilisons le lemme 3 pour majorer les deux 
sommations de $(46'')$~: on obtient
$$
[[(46'')]]\leq {\log 2-{1\over 2(4x+1)}\over 4x}
+{{3\over 2}\log 2+{1\over 2}(\log x+\gamma)+{1\over 24x^2}\over(4x)^2}<{0.234\over x}
$$
à l'aide d'une étude numérique pour les valeurs $x=1,2\ldots$~. Par des 
majorations plus brutales, on a par ailleurs 
$$
\eqalign{
\sum_{\ell=1}^{2x-1}{1\over 2(\ell{+}1)}\,{(4x)^{\ell-1}\over
4x\cdots(4x+\ell)}&\leq
\sum_{\ell=1}^{2x-1}{1\over 2(\ell{+}1)}\,{1\over (4x)^2}
\leq{\log2x{+}\gamma{+}{1\over 4x}{+}{1\over 12(2x)^2}{-}1\over 32x^2}
<{0.025\over x},\cr
&\kern-122pt\sum_{\ell=3}^{2x-1}{1\over 2(\ell-2)}
{4x(4x-1)\cdots(4x-\ell+1)\over(4x)^{\ell+2}}\leq
\sum_{\ell=1}^{2x-3}{1\over\ell}\,{1\over 32x^2}
\leq{\log2x+\gamma\over 32x^2}<{0.040\over x},\cr}
$$
ce qui donne au total
$$
|R_4(x)|\leq {0.299\over x}.\leqno(50)
$$
Pour borner optimalement les autres termes, il nous faut estimer 
plus précisément les produits partiels
$\prod(4x\pm j)$ qui apparaissent, et pour cela on utilise des encadrements
des logarithmes. Pour $t>0$, on a $t-{1\over 2}\,t^2<\log(1+t)<t$. Pour
$t={j\over 4x}$, ceci donne
$$
-{\sum\limits_{1\leq j\leq\ell}~j\over 4x}<
\log{(4x)^\ell\over(4x+1)\cdots(4x+\ell)}=
\sum_{1\leq j\leq\ell}\log{1\over 1+{j\over 4x}}<
-{\sum\limits_{1\leq j\leq\ell}~j\over 4x}+
{\sum\limits_{1\leq j\leq\ell}~j^2\over 2(4x)^2}
$$
où $\sum_{1\leq j\leq\ell}~j={\ell(\ell+1)\over 2}$ et
$\sum_{1\leq j\leq\ell}~j^2={\ell(\ell+1)(2\ell+1)\over 6}$,
donc
$$
-{\ell(\ell{+}1)\over 8x}<\log{(4x)^\ell\over(4x{+}1)\cdots(4x{+}\ell)}
<-{\ell(\ell{+}1)\over 8x}+{\ell(\ell{+}1)(2\ell{+}1)\over 12\,(4x)^2},
$$
soit encore
$$
\exp\bigg({1\over 32x}{-}{(\ell{+}1/2)^2\over 8x}\bigg)<
{(4x)^\ell\over(4x{+}1)\cdots(4x{+}\ell)}<\exp\bigg({1\over 32x}
{-}{(\ell{+}1/2)^2\over 8x}{+}{(\ell{+}1/2)^3\over 96x^2}\bigg).\kern-8pt
\leqno(51)
$$
Pour $\ell\leq 2x-1$, nous avons
$$
{(\ell+1/2)^2\over 8x}-{(\ell+1/2)^3\over 96x^2}
={(\ell+1/2)^2\over 8x}\bigg(1-{(\ell+1/2)\over 12x}\bigg)\geq
{5\over 6}\,{(\ell+1/2)^2\over 8x}
$$
de sorte que (à l'aide d'un nouveau calcul numérique)
$$
\leqalignno{
\qquad&{(4x)^\ell\over(4x+1)\cdots(4x+\ell)}<
\exp\bigg({1\over 32x}-{5\over 6}\,{(\ell+1/2)^2\over 8x}\bigg)
\kern53pt\hbox{pour $\ell\leq 2x-1$,}
&(51')\cr
&{(4x)^\ell\over(4x+1)\cdots(4x+\ell)}<\exp\bigg({1\over 32x}-{5\over 6}
{(2x-1/2)^2\over 12x}\bigg)<{1.52\over x}\qquad\hbox{pour $\ell\geq 2x-1$.}\cr}
$$
Pour $\ell\geq 2x$, chaque nouveau facteur ${4x\over 4x+\ell}$ est inférieur
où égal à ${2\over 3}$, donc
$$
\sum_{\ell=2x}^{+\infty}
{(4x)^\ell\over(4x+1)\cdots(4x+\ell)}<{1.52\over x}\sum_{p=1}^{+\infty}
\bigg({2\over 3}\bigg)^p<{3.04\over x}.
\leqno(51'')
$$
D'autre part, l'inégalité analogue $-t-16\,t^2/26<\log(1-t)<-t$ pour $t={j\over 4x}\leq 1/4$ implique
$$
-{\ell(\ell+1)\over 8x}-{16\,\ell(\ell+1)(2\ell+1)\over 6\cdot 26\,(4x)^2}
<\log{(4x-1)\cdots(4x-\ell)\over (4x)^\ell}<-{\ell(\ell+1)\over 8x}.
\leqno(52)
$$
Comme $\exp(1/4x)>1+1/4x$, on en tire
$$
~~{(4x+1)\cdots(4x-\ell+2)\over(4x)^\ell}\leq
\bigg(1+{1\over 4x}\bigg)\,\exp\bigg(-{(\ell-1)(\ell-2)\over 8x}\bigg)
<\exp\bigg(-{\ell(\ell-3)\over 8x}\bigg)\kern-4pt\leqno(52')
$$
avec un rapport entre deux majorants consécutifs d'indices $\ell$, $\ell+1$
inférieur à\break $\exp(-(2\ell-2)/8x)\leq e^{-1/4}$ si $\ell=2x$ 
et${}\leq e^{-1/2}$ si $\ell\geq 2x+1$, donc
$$
\sum_{\ell=2x}^{+\infty}{(4x+1)\cdots(4x-\ell+2)\over(4x)^\ell}
\leq \exp\bigg({3\over 4}-{x\over 2}\bigg)\Bigg(1+e^{-1/4}\sum_{p=0}^{+\infty}
e^{-p/2}\Bigg)<{4.65\over x},
$$
Comme $2\ell\geq 4x$, on en déduit d'après (46)
$$
|R_2(x)|\leq{2\over\pi}\,{1\over 4x}\,{7.69\over x}<{1.224\over x^2}\leqno(53)
$$
(la décroissance de $R_2(x)$ étant même exponentielle en fait).
Les inégalités $(43'')$, $(51')$ et~$(52')$ donnent aussi, au moyen d'une
comparaison série-intégrale
$$
\leqalignno{
|R_5(x)|&\leq{1\over 8}\sum_{\ell=1}^{2x-1}{\ell+1\over 4x(4x-1)}\,
\exp\bigg({1\over 32x}-{5\over 6}\,{(\ell+1/2)^2\over 8x}\bigg)
+2\,{\ell-1\over (4x)^2}\,\exp\bigg({3\ell\over 8x}-{\ell^2\over 8x}\bigg)\cr
&\leq{1\over 8(4x)(3x)}\Bigg(e^{1\over 32}\!\!
\int_0^{+\infty}\!\bigg(t\,{+}{3\over 2}\bigg)
\exp\bigg(\!-{5\over 6}{t^2\over 8x}\bigg)dt+{3\,e^{3\over 4}\over 2}\!
\int_0^{+\infty}t\,\exp\bigg(\!-{t^2\over 8x}\bigg)dt\!\Bigg)\cr
&={1\over 96x^2}\Bigg(e^{1\over 32}\bigg({24\over 5}x+{3\over 2}
\sqrt{48x\over 5}\,{1\over 2}\,\sqrt{\pi}\bigg)+6\,e^{3\over 4}\,x\Bigg)
<{0.229\over x}\qquad
\hbox{pour $x\geq 1$}.&(54)\cr}
$$
Il vient ensuite d'après (43) et (51)
$$
\eqalign{
T'(x)=\sum_{\ell=1}^{2x-1}{1\over 2\ell}\bigg(
&{(4x)^\ell\over(4x+1)\cdots(4x+\ell)}-
{(4x-1)\cdots(4x-\ell)\over(4x)^\ell}\bigg)\cr
&=\sum_{\ell=1}^{2x-1}{1\over2\ell}\,{(4x)^\ell\over(4x+1)\cdots(4x+\ell)}\,
\bigg(1-\prod_{j=1}^\ell\Big(1-{j\over 4x}\Big)\Big(1+{j\over 4x}\Big)\bigg)\cr
&\leq\sum_{\ell=1}^{2x-1}
\exp\bigg({1\over 32x}-{(\ell+1/2)^2\over 8x}+{(\ell+1/2)^3\over 96\,x^2}
\bigg)\,{(\ell+1)^2\over 96x^2}\cr}
$$
grâce aux inégalités $1-\prod(1-a_j)\leq\sum a_j$ avec
$a_j={j^2\over (4x)^2}<1$, et 
$\sum_{j\leq \ell}j^2={\ell(\ell+1)(2\ell+1)\over 6}$.
Dans l'autre sens, on a $\prod(1-a_j)^{-1}-1\geq\sum a_j$, donc
d'après (52)
$$
\eqalign{
T'(x)=\sum_{\ell=1}^{2x-1}{1\over 2\ell}\,&
{(4x-1)\cdots(4x-\ell)\over(4x)^\ell}\Bigg(\prod_{j=1}^\ell
\bigg(1-\Big({j\over 4x}\Big)^2\bigg)^{-1}-1\Bigg)\cr
&\qquad{}\geq\sum_{\ell=1}^{2x-1}
\exp\bigg(-{\ell(\ell+1)\over 8x}-{(\ell+1/2)^3\over 78\,x^2}
\bigg)\,{(\ell+1)(2\ell+1)\over 12\,(4x)^2}\cr
&\qquad{}\geq\sum_{\ell=1}^{2x-1}
\exp\bigg(-{(\ell+1/2)^2\over 8x}-{(\ell+1/2)^3\over 78\,x^2}
\bigg)\,{(\ell+1)(\ell+1/2)\over 96x^2}\cr
&\qquad{}\geq\sum_{\ell=1}^{2x-1}
\exp\bigg(-{(\ell+1/2)^2\over 8x}\bigg)\bigg(1-{(\ell+1/2)^3\over 78\,x^2}
\bigg)\,{(\ell+1)(\ell+1/2)\over 96x^2}.\cr}
$$
Nous évaluons maintenant ces sommes en les remplaçant par des intégrales.
Il vient
$$
T'(x)\leq e^{1\over 32x}\int_0^{2x}
\exp\bigg(-{t^2\over 8x}+{t^3\over 96x^2}\bigg)\,{(t+3/2)^2\over 96x^2}\,dt
$$
en majorant le terme d'indice $\ell$ par l'intégrale sur $[\ell-1/2,\ell+1/2]$.
Le changement de variable
$$
u={t^2\over 8x}-{t^3\over 96x^2}={t^2\over 8x}\bigg(1-{t\over 12x}\bigg),
\qquad du={t\over 4x}\bigg(1-{t\over 8x}\bigg)dt
$$
implique $u\geq{5\over 48x}\,t^2$, donc
$t\leq\sqrt{48x\over 5}\,\sqrt{u}$. De plus, on a par convexité 
$(1-{v\over p})^{-1}\leq 1+{1\over p-1}v$ si $v\leq 1\;$; 
on prend ici $v={t\over 2x}$ et $p=6$, resp.\ $p=3$) pour voir que
$$
\eqalign{
&t=\sqrt{8xu}\,\bigg(1-{t\over 12x}\bigg)^{-1/2}
\leq\sqrt{8xu}\,\bigg(1+{t\over 20x}\bigg)
\leq\sqrt{8xu}\,\bigg(1+{\sqrt{3\over 125x}}\,\sqrt{u}\bigg),\cr
&dt={4x\over t}\bigg(1-{t\over 8x}\bigg)^{-1}du\leq
{4x\over t}\bigg(1+{t\over 6x}\bigg)du
\leq{4x\over \sqrt{8x u}\,
}\bigg(1+{2\over\sqrt{15x}}\,\sqrt{u}\bigg)du,\cr}
$$
par conséquent
$$
T'(x)\leq
{e^{1\over 32 x}\over 96x^2}\int_0^{+\infty}
e^{-u}\bigg({3\over 2}+\sqrt{8xu}\bigg(1+\sqrt{3\over 125x}\sqrt{u}\bigg)\bigg)^2
\bigg(1+{2\over\sqrt{15x}}\,\sqrt{u}\bigg)
{\sqrt{2x}\,du\over\sqrt{u}}.
$$
Cette intégrale s'évalue explicitement, avec un terme dominant qui est
$$
{e^{1\over 32 x}\over 96x^2}\int_0^{+\infty}
e^{-u}(\sqrt{8xu})^2{\sqrt{2x}\,du\over\sqrt{u}}
\sim{\sqrt {2}\over 12\sqrt{x}}\int_0^{+\infty}e^{-u}\,\sqrt{u}\,du
={\sqrt{2\pi}\over 24\,x^{1/2}}.
$$
Le facteur $e^{1\over 32x}$ est majoré par $1+{1\over 31.5\,x}$, ce qui amène 
une erreur
$$
{\sqrt{2\pi}\over 24\,x^{1/2}}\cdot {1\over 31.5\,x}<{0.004\over x}.
$$
Tous les autres termes intervenant dans l'intégrale mettent en jeu des termes
$O({1\over x})$, avec des coefficients qui sont des produits de facteurs
$\Gamma(a)$, ${1\over 2}\leq a\le 2$, par des coefficients dont la somme 
est majorée par
$$
{e^{1\over 32}\over 96}\bigg[
\bigg({3\over 2}+\sqrt{8}\bigg(1+\sqrt{3\over 125}\,\bigg)\bigg)^2
\bigg(1+{2\over\sqrt{15}}\bigg)\sqrt{2}-8\sqrt{2}\,\bigg]
<0.4021.
$$
Comme $\Gamma(a)\leq\sqrt{\pi}$, on obtient
$$
T'(x)<{\sqrt{2\pi}\over 24\,x^{1/2}}+{0.717\over x}.
$$
De même, la minoration de $T'(x)$ donne
$$
\eqalign{
T'(x)&\geq\sum_{\ell=1}^{2x-1}
\exp\bigg(-{(\ell+1/2)^2\over 8x}\bigg)\,
\bigg(1-{(\ell+1/2)^3\over 78\,x^2}\bigg){(\ell+1)(\ell+1/2)\over 96x^2}\cr
&\geq\int_{3/2}^{2x+1/2}
\exp\bigg(-{t^2\over 8x}\bigg)\bigg(1-{t^3\over 78\,x^2}
\bigg)\,{(t-1)(t-1/2)\over 96x^2}\,dt\cr
&\geq\int_{2}^{2x}
\exp\bigg(-{t^2\over 8x}\bigg)\bigg(1-{t^3\over 78\,x^2}
\bigg)\,{t^2-3t/2\over 96x^2}\,dt\cr
&=\int_{1/2x}^{x/2}
e^{-u}\bigg(1-{8\sqrt{8}\,u^{3/2}\over 78\,x^{1/2}}
\bigg)\,{8xu-3\sqrt{8}\,x^{1/2}u^{1/2}/2\over 96x^2}\,{
\sqrt{8}\,x^{1/2}\,du\over 2\,u^{1/2}}\cr
&\geq\int_{1/2x}^{x/2}
e^{-u}\bigg(1-{8\sqrt{8}\,u^{3/2}\over 78\,x^{1/2}}
\bigg)\,{\sqrt{8}\,u-3\,x^{-1/2}u^{1/2}/2\over 24\,x^{1/2}}\,{du\over u^{1/2}}\cr
&\geq\int_{1/2x}^{x/2}
e^{-u}\bigg({\sqrt{2}\,u^{1/2}\over 12\,x^{1/2}}-
{8\,u^2\over 3\cdot 78\,x}-{1\over 16x}\bigg)\,du\cr
&\geq\int_{0}^{+\infty}
e^{-u}\bigg({\sqrt{2}\,u^{1/2}\over 12\,x^{1/2}}-
{4\,u^2\over 117\,x}-{1\over 16x}\bigg)\,du-
\int_{\complement}
e^{-u}{\sqrt{2}\,u^{1/2}\over 12\,x^{1/2}}\,du.\cr}
$$
L'intégrale complémentaire $\int_{\complement}...$ est majorée 
sur $[0,1/2x]$ par
$$
\int_{0}^{1/2x}{\sqrt{2}\,u^{1/2}\over 12\,x^{1/2}}\,du
={1\over 36\,x^2},
$$
tandis que sur $[A,+\infty[{}=[x/2,+\infty[$ on a l'inégalité
$$
\int_A^{+\infty}u^\alpha\,e^{-u}\,du=A^\alpha\,e^{-A}+
\int_A^{+\infty}\alpha\,u^{\alpha-1}\,e^{-u}\,du
\leq e^{-A}(A^\alpha+\alpha A^{\alpha-1}),\quad\alpha\in{}]0,1],
$$
ce qui fournit la majoration
$$
\int_{{x\over 2}}^{+\infty}
e^{-u}{\sqrt{2}\,u^{1/2}\over 12\,x^{1/2}}\,du
\leq\exp\bigg(-{x\over 2}\bigg)\bigg({1\over 12}+
{1\over 12x}\bigg)
\leq{{1\over 6}\,e^{-1/2}\over x}.
$$
On obtient ainsi la borne inférieure explicite
$$
T'(x)>{\sqrt{2\pi}\over 24\,x^{1/2}}-\bigg({8\over 117}+{1\over 16}+{1\over 36}
+{1\over 6}e^{-1/2}\bigg){1\over x}>{\sqrt{2\pi}\over 24\,x^{1/2}}-
{0.260\over x}.
$$
De la même manière, mais cette fois, sans compensation de termes, avec
des calculs beaucoup plus simples, et grâce à $(43'')$, (51), (52),
on trouve une majoration
$$
T''(x)\leq{1\over 4x}\sum_{\ell=1}^{2x-1}
{1\over 4}\exp\bigg({32\over x}-{(\ell{+}1/2)^2\over 8x}+{(\ell{+}1/2)^3\over 
96x^2}\bigg)+{3\over 4}\exp\bigg({32\over x}-{(\ell{-}1/2)^2\over 8x}\bigg).
$$
Par des majorations intégrales analogues à celles déjà utilisées, ceci donne
$$
\eqalign{
T''(x)&\leq{e^{32\over x}\over 4x}\Bigg(
{1\over 4}\int_0^{2x}\exp\bigg(-{t^2\over 8x}+{t^3\over 96x^2}\bigg)dt
+{3\over 4}\int_0^{2x}
\exp\bigg(-{t^2\over 8x}\bigg)dt\Bigg)+{3\over 16x}\cr
&\leq{e^{32\over x}\over 4x}\Bigg(
{1\over 4}\int_0^{+\infty}e^{-u}\bigg(1+{2\over\sqrt{15x}}\,\sqrt{u}\bigg)
{\sqrt{2x}\,du\over \sqrt{u}}
+{3\over 4}\int_0^{+\infty}e^{-u}\,{\sqrt{2x}\,du\over \sqrt{u}}
\Bigg)+{3\over 16x}\cr
&\leq{e^{32\over x}\over 4x}\int_0^{+\infty}e^{-u}\,{\sqrt{2x}\,du\over \sqrt{u}}
+{e^{32\over x}\over 4x}\,{1\over\sqrt{30}}+{3\over 16x}
<{\sqrt{2\pi}\over 4x^{1/2}}+{0.255\over x},\cr}
$$
et une minoration
$$
\eqalign{
T''(x)&\geq{1\over 4x}\sum_{\ell=1}^{2x-1}
{1\over 4}\exp\bigg(-{(\ell+1/2)^2\over 8x}\bigg)+
{3\over 4}\exp\bigg(-{(\ell-1/2)^2\over 8x}
-{(\ell-1/2)^3\over 78x^2}
\bigg)\cr
&\geq{1\over 4x}\Bigg({1\over 4}\int_{3/2}^{2x+1/2}
\exp\bigg(-{t^2\over 8x}\bigg)dt
+{3\over 4}\int_{1/2}^{2x-1/2}
\exp\bigg(-{t^2\over 8x}\bigg)\bigg(1-{t^3\over 78x^2}\bigg)dt\Bigg)\cr
&\geq{1\over 4x}\Bigg({1\over 4}\int_0^{2x}
\exp\bigg(-{t^2\over 8x}\bigg)dt
+{3\over 4}\int_{0}^{2x}
\exp\bigg(-{t^2\over 8x}\bigg)\bigg(1-{t^3\over 78x^2}\bigg)dt-{9\over 8}\Bigg)\cr
&\geq{1\over 4x}\Bigg(\int_{0}^{2x}
\exp\bigg(-{t^2\over 8x}\bigg)dt
-{3\over 4}\int_0^{+\infty}
\exp\bigg(-{t^2\over 8x}\bigg)\,{t^3\over 78x^2}\,dt-{9\over 8}\Bigg)\cr
&={1\over 4x}\Bigg(\int_0^{x/2}e^{-u}\,{\sqrt{2x}\,du\over\sqrt{u}}
-{1\over 104}\int_0^{+\infty}e^{-u}{u\,du\over 2}-{9\over 8}\Bigg)\cr
&\geq{1\over 4x}\Bigg(\int_0^{+\infty}e^{-u}\,{\sqrt{2x}\,du\over\sqrt{u}}
-{235\over 208}-2\,e^{-x/2}\Bigg)>{\sqrt{2\pi}\over 4x^{1/2}}-{0.586\over x}.
\cr}
$$
Nous obtenons finalement l'estimation
$$
T'(x)-T''(x)=-{5\over 24}\,{\sqrt{2\pi}\over x^{1/2}}+R_8(x),\qquad
-{0.515\over x}<R_8(x)<{1.303\over x}.
\leqno(55)
$$
Il reste seulement $U'(x)$ à évaluer. Or, d'après (45),
un changement de variable $\ell=\ell'+1$ suivi d'une
décomposition $4x=(4x-\ell)+\ell$ permet de transformer la deuxième
sommation figurant dans $U'(x)$~:
$$
\eqalign{
U'(x)&=\sum_{\ell=1}^{2x-1}
{1\over 2(\ell+1)}\,{(4x)^\ell\over 4x\cdots(4x+\ell)}-
\sum_{\ell=1}^{2x-2}{1\over 2\ell}\,
{4x(4x-1)\cdots(4x-\ell+1)\over(4x)^{\ell+1}}\cr
&=\sum_{\ell=1}^{2x-1}
{1\over 2(\ell+1)}\,{(4x)^{\ell-1}\over (4x+1)\cdots(4x+\ell)}-
\sum_{\ell=1}^{2x-2}{1\over 2\ell}\,
{(4x-1)\cdots(4x-\ell+1)(4x-\ell)\over(4x)^{\ell+1}}\cr
&\kern160pt{}-\sum_{\ell=1}^{2x-2}{1\over 2}\,
{(4x-1)\cdots(4x-\ell+1)\over(4x)^{\ell+1}}.\cr}
$$
En écrivant maintenant ${1\over \ell+1}={1\over \ell}-{1\over \ell(\ell+1)}$,
il vient
$$
U'(x)={1\over 4x}\,T'(x)-R_9(x)
$$
avec
$$
\eqalign{
R_9(x)&=\sum_{\ell=1}^{2x-1}
{1\over 2\ell(\ell+1)}\,{(4x)^{\ell-1}\over (4x+1)\cdots(4x+\ell)}
+\sum_{\ell=1}^{2x-2}{1\over 2}\,
{(4x-1)\cdots(4x-\ell+1)\over(4x)^{\ell+1}}\cr
&\kern100pt{}-\bigg({1\over 2\ell}{(4x-1)\cdots(4x-\ell)
\over(4x)^{\ell+1}}\bigg)_{\ell=2x-1}\cr}
$$
et, pour $x\geq 2$, une majoration
$$
0<R_9(x)<{1\over 4x}\sum_{\ell=1}^{+\infty}{1\over 2\ell(\ell+1)}
+{1\over 2}(2x-2){1\over (4x)^2}<{3\over 16x}.
$$
Il résulte alors de l'encadrement de $T'(x)$ (et d'un calcul explicite 
de $U'(x)$ pour \hbox{$x=1{,}2{,}3$}) que
$$
|U'(x)|<{0.206\over x}.\leqno(56)
$$
En combinant (25), $(31')$, (47), $(48'')$, (49), (50) et (53 -- 56), on 
obtient
$$
\Delta(x)={e^{-4x}\over 4\pi x}\bigg(
-{5\sqrt{2\pi}\over 12\,x^{1/2}}+R(x)\bigg)\leqno(57)
$$
avec
$$
R(x)=-U'(x)
+\pi\Big(R_1(x)+R_2(x)-R_3(x)\Big)-{5\over 4}R_4(x)+2\,R_5(x)-R_6(x)+R_7(x)
+2\,R_8(x),
$$
d'où
$$
|R(x)|<{10.835\over x}.
\leqno(58)
$$
Ces estimées impliquent (12), (13). La preuve du théorème est achevée.\qed
\medskip

{\bf Complexité numérique de la méthode de Brent-McMillan pour le calcul
de la constante d'Euler.} La version raffinée permet d'atteindre une
précision finale meilleure que $e^{-8x}$. Ceci amène à choisir 
$x = {1\over 8}\,d\,\log 10$ et conduit ainsi au temps de calcul
$$
B'(d) = \Big({3\over 4}a_3 + {1\over 2}\Big)\,\log 10\;d^2\simeq 9.7\,d^2,
$$
optimal parmi tous les algorithmes présentés ici.

Nous terminons en donnant le \lguil\?hit-parade\?\rguil\ de ces 
algorithmes, rangés par ordre d'effi\-cacité croissante.

$$\table{\quad
&\hrulefill&\hrulefill&\hrulefill&\hrulefill&\hrulefill&\cr
&\vertspace~\hbox{\rm Algorithmes}~&\vertspace\kern29pt\hbox{\rm Euler}
&\vertspace\kern57pt\hbox{\rm Sweeney}
&\vertspace\kern5pt\hbox{\rm Brent-McM.}~&\vertspace\cr
\noalign{\vskip-1pt}
&\vertspace&
\vertspace\rl{\kern15pt E_1}\kern42pt
\vertspace\rl{\kern15pt E_2}\kern42pt&
\vertspace\rl{\kern8pt S_1}\kern25pt
\vertspace\rl{\kern8pt S_2}\kern25pt
\vertspace\rl{\kern8pt S_1'}\kern25pt
\vertspace\rl{\kern8pt S_2'}\kern25pt
\vertspace\rl{\kern8pt S_3}\kern25pt
\vertspace\rl{\kern8pt S_3'}\kern25pt
&\vertspace\rl{\kern11pt B}\kern29pt\vertspace\rl{\kern11pt B'}\kern29pt
&\vertspace&\cr
&\hrulefill&\hrulefill&\hrulefill&\hrulefill&\hrulefill&\cr
&\Vertspace
\Raise{\smash{~{\hbox{\rm Temps de}~~ \atop ~\hbox{\rm calcul / $d^{\,2}$}}}}
&\Vertspace\rl{\Raise{~\displaystyle{Cd\over (\log d)^2}~}}\kern42pt
\Vertspace\rl{\Raise{~C\,\log d}}\kern42pt&
\Vertspace\rl{\Raise{~49.6}}\kern25pt
\Vertspace\rl{\Raise{~37.6}}\kern25pt
\Vertspace\rl{\Raise{~26.7}}\kern25pt
\Vertspace\rl{\Raise{~26.0}}\kern25pt
\Vertspace\rl{\Raise{~22.7}}\kern25pt
\Vertspace\rl{\Raise{~16.9}}\kern25pt
&\Vertspace\rl{\Raise{~~12.4}}\kern29pt
\Vertspace\rl{\Raise{~~9.7}}\kern29pt&\Vertspace\cr
&\hrulefill&\hrulefill&\hrulefill&\hrulefill&\hrulefill&\cr}$$

Signalons qu'il existe des algorithmes théoriquement encore plus
rapides, permettant d'évaluer $\gamma$ à $10^{-d}$ près en
$O(d(\log d)^3\log\log d)$ 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^x$, $I_0(x)$, $S_0(x)$)~; voir Brent~[6]. De tels 
algorithmes \lguil\?rapides\?\rguil\ sont toutefois difficiles 
à programmer et ne l'emportent sur les algorithmes 
\lguil\?classiques\?\rguil\ présentés ici que lorsque
$d$ est très grand.
\bigskip

\section{4}{Développement en fraction continue de 
$\hbox{\bigbfgreek\char'015}$}

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'à
l'ordre 29$\,$000 (Brent - Mac Millan [11]).  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 [16]) :
$$
f_n= \log_2\bigg(1+{1\over n}\bigg) - \log_2\bigg(1+{1\over n+1}\bigg)
$$
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  
suivant, qui rend extrêmement  improbable la rationalité de $\gamma$ ou
$e^\gamma$~:
\medskip

{\statement Théorème.\pointir} {\it Si l'un des nombres $\gamma$ ou $e^\gamma$ 
est un rationnel $p/q$ pour des entiers $p,q$ positifs, alors 
$q > 10^{15\,000}$.}
\bigskip 

{\bigbf Bibliographie}
{
\parindent=7.5mm
\medskip

\item{[1]} {\petcap J.C.\ Adams} -- {\it On the value of Euler's constant}~;
Proc.\ Roy.\ Soc.\ London, {\bf 27} (1878) p.~88--94. Voir aussi vol.~{\bf 42}
(1887), p.~22--25.

\item{[2]} {\petcap P.\ Appell} -- {\it Sur la nature arithmétique de la 
constante d'Euler}~; Comptes Rend.\ Ac.\ Sc.\ Paris, v.~{\bf 182}, 
12 avril 1926, p.~897--899 ; 19 avril 1926, p.~949.

\item{[3]} {\petcap W.A.\ Beyer \&\ M.\ S.\ Waterman} -- {\it Error analysis
of a computation of Euler's constant}~; Math, of Comp., v.~{\bf 28} 
(1974) p.~599--604.

\item{[4]} {\petcap W.A.\ Beyer \&\ M.S.\ Waterman} -- {\it Decimals and
partial quotients of Euler's constant and $\ln 2$}~; UMT 19, Math, of Comp., 
v.~{\bf 28}, 1974, p.~667.\ Errata : Math, of Comp., MTE 549, v.~{\bf 32}
(1978), p.~317--318.

\item{[5]} {\petcap N.\ Bourbaki} -- {\it Fonctions d'une variable réelle}~;
chap.~VI, \S1, n${}^\circ${\bf 7}; Hermann, Paris, 1951.

\item{[6]} {\petcap R.P.\ Brent} -- {\it The complexity of multiple-precision
arithmetic}~; Complexity of Computational Problem Solving (R.S.\ 
Anderssen and R.P.\ Brent, Editors), Univ.\ of Queensland Press, 
Brisbane, 1976, p.~126--165.

\item{[7]} {\petcap R.P.\ Brent} -- {\it Computation of the regular continued
fraction for Euler's constant}~; Math, of Comp., v.~{\bf 31}, July 1977, 
p.~771--777.

\item{[8]} {\petcap R.P.\ Brent} -- {\it $\gamma$ and $\exp(\gamma)$ to 
20,700 D and their regular continued fractions to 20,000 partial quotients}~;
UMT 1, Math.\ of Comp., v.~{\bf 32}, 1978, p.~311.

\item{[9]} {\petcap R.P.\ Brent} -- {\it Euler's constant and its exponential
to 30,100 decimals}~; Math.\ of Comp., UMT File.

\item{[10]} {\petcap R.P.\ Brent \&\ E.M.\ McMillan} -- {\it Some new 
algorithms for high-precision computation of Euler's constant}~;
Math.\ of Comp., v.~{\bf 34}, January 1980, p.~305--312.

\item{[$10'$]} {\petcap R.P.\ Brent \&\ F.\ Johansson} -- {\it  A bound
for the error term in the Brent-McMillan algorithm}~;
Math.\ of Comp., v.~{\bf 84}, 2015, p.~2351--2359.

\item{[11]} {\petcap R.P.\ Brent \&\ E.M.\ McMillan} -- {\it The first 
29,000 partial quotients in the regular continued fraction for Euler's 
constant and its exponential}~; Math, of Comp., UMT File.

\item{[12]} {\petcap L.\ Euler} -- {\it De progressionibus harmonicis
observationes}~; Euleri Opera Omnia, Ser.\ 1, v.~{\bf 14},
Teubner, Leipzig and Berlin, 1925, p.~93--100.

\item{[13]} {\petcap L.\ Euler} -- {\it De summis serierum numeros 
Bernoullianos involventium}~; Euleri Opera 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.

\item{[14]} {\petcap A.\ Froda} -- {\it La constante d'Euler est 
irrationnelle}~; Atti Accad.\ Naz.\ Lincei Rend.\ Cl.\ Sci.\ Fis.\ 
Mat.\ Natur.\ (8) {\bf 38} (1965), p.~338--344.

\item{[15]} {\petcap J.W.L.\ Glaisher} -- {\it History of Euler's constant}~; 
Messenger of Math., v.~{\bf 1} (1872), p.~25--30.

\item{[16]} {\petcap A.\ Ya.\ Khintchine (A.\ Ja.\ Hin\v{c}in)} -- {\it
Continued Fractions}~; 3${}^{\rm e}$ édition, traduction anglaise par P.\ Wynn, 
Noordhoff, Groningen, 1963.

\item{[17]} {\petcap D.E.\ Knuth} -- {\it Euler's constant to 1271 places}~;
Math, of Comp., v.~{\bf 16} (1962), p.~275--281.

\item{[18]} {\petcap K.\ Mahler} -- {\it Applications of a theorem by 
A.B.\ Shidlovskii}~; Proc.\ Roy.\ Soc.\ London, A~{\bf 305} (1968), p.~149--173.

\item{[19]} {\petcap W.\ Shanks} -- {\it On the numerical value of Euler's
constant}~; Proc.\ Roy.\ Soc.\ London, v.~{\bf 15} (1867), p.~429--432 ; 
v.~{\bf 20} (1871), p.~29--34.

\item{[20]} {\petcap A.\ Schönhage \&\ V.\ Strassen} -- {\it Schnelle
Multiplikation grosser Zahlen}~; Computing, v.~{\bf 7} (1971), p.~281--292.

\item{[21]} {\petcap D.W.\ Sweeney} -- {\it On the computation of Euler's 
constant}~; Math, of Comp., v.~{\bf 17} (1963), p.~170--178.

\item{[22]} {\petcap G.N.\ Watson} -- {\it A Treatise on the Theory of
Bessel Functions}~; 2${}^{\rm nd}$  edition, Cambridge Univ.\ Press,
London, 1944.

\item{[23]} {\petcap J.W.\ Wrench, Jr} -- {\it A new calculation of Euler's
constant}~; Math.\ Tables \&\ other Aids to Comp., v.~{\bf 6} (1952), p.~255.
\bigskip

\noindent
{\bf Ouvrages de référence sur la fonction $\Gamma$.}
\medskip

\item{[24]} {\petcap E.\ Artin} -- {\it The Gamma Function}~;
Holt, Rinehart and Winston, 1964, traduit de~: {\it Einführung in die 
Theorie der Gammafunktion}, Teubner, 1931.

\item{[25]} {\petcap R.\ Campbell} -- {\it Les intégrales eulériennes et 
leurs applications}~; Dunod, Paris, 1966.
\bigskip\bigskip}

Jean-Pierre Demailly\\
Université de Grenoble I, Institut Fourier\\
Laboratoire de Mathématiques Pures associé au CNRS\\
 B.P.\ 74 -- Tél.\ (76) 51-46-00\\
38402 Saint-Martin-d'Hères Cedex
\bigskip\bigskip

\centerline{(juin 1984~; octobre 2016~: ajout des calculs non initialement 
publiés)}

\end

% Local Variables:
% TeX-command-default: "UpTeX"
% End:
