# Nombres aléatoires, suites arithmétiques, méthode de Monte-Carlo

J. Bass

► To cite this version:

J. Bass. Nombres aléatoires, suites arithmétiques, méthode de Monte-Carlo. Annales de l'ISUP, 1960, IX (2), pp.289-325. hal-04092917

HAL Id: hal-04092917

https://hal.science/hal-04092917

Submitted on 9 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# NOMBRES ALÉATOIRES, SUITES ARITHMÉTIQUES, MÉTHODE DE MONTE-CARLO

J. BASS

## SOMMAIRE

La méthode de Monte-Carlo est une méthode de calcul numérique faisant appel à des analogies probabilistes. Elle utilise des nombres aléatoires qui sont les résultats d'épreuves successives sur une variable aléatoire donnée. Les nombres aléatoires les plus simples sont ceux dont la distribution est uniforme entre 0 et 1. Il en existe des tables, construites par des procédés difficiles à justifier. On montre, en prenant comme modèle le calcul d'une intégrale définie, qu'il y a intérêt à substituer aux méthodes probabilistes des méthodes arithmétiques. Ces méthodes reposent sur les propriétés des suites de nombres réels uniformément denses dans un intervalle. Les théorèmes de H. Weyl en précisent la structure et en donnent des exemples utilisables. Le plus simple est celui des multiples modulo 1 d'un nombre irrationnel A. La discussion de l'erreur

commise en remplacant l'intégrale simple $\int_{0}^{1} f(x) dx$ par $\frac{1}{N} \sum_{n=0}^{N-1} f(x_n)$,

où $x_n = An$, mod 1, montre que la convergence est plus rapide qu'avec des nombres aléatoires véritables. On établit une formule permettant de calculer dans ce cas une limite supérieure de l'erreur commise. Lorsque A est algébrique, et moyennant des conditions de régularité sur $f(x)$, des

formules de Richtmyer montrent que l'erreur est proportionnelle à $\frac{1}{N}$. Une discussion est amorcée pour discuter l'ordre de grandeur de l'erreur lorsque -

$x_n = An^2$ mod 1. La précision semble moins bonne, mais la suite $An^2$ présente cependant quelques avantages sur la suite An: elle n'a pas de corrélation interne, et elle peut être utilisée pour former des nombres aléatoires prenant les valeurs 0 et 1 avec des probabilités égales.

## I - INTRODUCTION -

Donner une définition générale de ce qu'on appelle méthode de Monte-Carlo est une entreprise difficile. Il semble qu'il s'agisse d'une méthode de calcul analogique, dans laquelle on remplace un phénomène naturel, présentant certains caractères aléatoires ou statistiques, par une image de nature probabiliste et mathématique de telle sorte que, par des opérations statistiques pures, on puisse simuler le phénomène physique aléatoire, et en calculer les paramètres d'une façon plus économique que par des expériences directes.
Une telle définition est très vague, et je n'ai pas l'intention de la discuter ici. Je me limiterai à un point de vue plus restreint. Je considérerai comme "méthode de Monte-Carlo" toute méthode qui permet de résoudre numériquement un problème mathématique par une analogie probabiliste. Le "problème mathématique" consiste à calculer un nombre ou un ensemble de nombres définis par des opérations d'analyse pure, où n'apparaît pas la notion de probabilité: calculs d'intégrales simples ou multiples, inversions de matrices, résolution d'équations intégrales, et cela sans qu'intervienne l'origine physique du problème. Il s'agit de découvrir, entre le formalisme du problème posé, et celui d'un certain problème de calcul des probabilités, une analogie telle que le second puisse servir à résoudre le premier. Or, résoudre un problème de probabilités, cela se traduit en pratique par une expérimentation statistique qui fournit des valeurs approchées du résultat, avec une exactitude d'autant meilleure que le nombre d'épreuves est plus grand. D'où des méthodes pratiques de calcul numérique, dans la mesure où l'on dispose de moyens matériels permettant d'effectuer, avec une précision suffisante, un grand nombre d'expériences successives.

La statistique fournit des fréquences et des valeurs moyennes. La méthode de Monte-Carlo consiste donc à interpréter, après des transformations convenables, le problème posé comme un calcul de fréquences ou de moyennes. Mais on peut toujours traduire une fréquence en langage de valeurs moyennes. D'une façon un peu schématique, on peut donc dire que le but de la méthode de Monte-Carlo est le calcul des valeurs moyennes, c'est-à-dire d'intégrales simples ou multiples.

Bien que cette conclusion soit loin d'être absolue, je vais la choisir comme base de départ, et je prendrai dans ce qui suit comme modèle de problème le calcul numérique d'une intégrale simple ou multiple. Si le domaine d'intégration est borné, on peut toujours le supposer intérieur à l'hypercube $C: 0 \le x < 1, 0 \le y < 1, \dots$, moyennant une transformation linéaire préalable. Si la fonction à intégrer est définie dans un domaine intérieur à cet hypercube, on peut toujours la prolonger à tout l'hypercube en lui donnant la valeur 0 dans la partie de C où elle n'est pas initialement définie. On est donc conduit à l'intégrale

$$I = \int_C f(x, y, \dots) \, dx dy \dots$$

## II - CALCUL D'UNE INTEGRALE -

Soit $f(x)$ une fonction bornée, intégrable, définie sur le segment $[0, 1]$. L'intégrale
$$I = \int_{0}^{1} f(x) dx, \qquad (2.1)$$

qui a priori n'a aucune relation avec le calcul des probabilités, peut être considérée comme la valeur moyenne (espérance mathématique) d'une variable aléatoire X uniformément distribuée entre 0 et 1.

D'après la loi des grands nombres, si $x_1, ..., x_N$ sont les résultats de N épreuves successives sur X, la moyenne arithmétique

$$S_N = \frac{1}{N} [f(x_1) + ... + f(x_N)] \qquad (2.2)$$

tend en probabilité vers $I$ lorsque $N \to \infty$. Il est d'ailleurs commode de considérer $x_1, ..., x_N$ comme les résultats d'épreuves successives sur N variables aléatoires indépendantes uniformément distribuées sur (0,1).

Désignons par $\sigma$ l'écart-type de $f(X)$, défini par

$$\sigma^2 = \int_{0}^{1} [f(x) - I]^2 dx = \int_{0}^{1} [f(x)]^2 dx - I^2. \qquad (2.3)$$

On sait que $S_N$ est une variable aléatoire qui a pour moyenne (stochastique) $I$, pour écart-type $\frac{\sigma}{\sqrt{N}}$, et que $\frac{\sqrt{N}}{\sigma}$ ($S_N - I$) est asymptotiquement distribuée suivant une loi normale réduite.

Par suite, si N est suffisamment grand,

$$\text{Prob.} [|S_N - I| > h] = \frac{2}{\sqrt{2\pi}} \cdot \frac{\int_{h\sqrt{N}}^{\infty} e^{-\frac{t^2}{2}} dt}{\sigma}$$

$|S_N - I|$ est la valeur absolue de l'erreur commise sur l'évaluation de $I$ lorsqu'on prend $S_N$ comme valeur approchée. Si $I$ on fait $4N$ épreuves au lieu de $N$, on a la même probabilité de commettre sur $I$ une erreur deux fois plus petite. C'est ce qu'on exprime en disant que l'erreur décroît comme $\frac{1}{\sqrt{N}}$. Si $I$ on a fait une évaluation grossière de $\sigma$, et qu'on se fixe un seuil de probabilité, on peut calculer le nombre de termes N nécessaire pour obtenir $I$ avec une précision donnée à l'avance. Si par exemple on néglige des éventualités de probabilité 0,05, le nombre de termes N permettant d'obtenir $I$ à 0,001 près est défini par

$$\frac{2}{\sqrt{2\pi}} \cdot \frac{\int_{h\sqrt{N}}^{\infty} e^{-\frac{t^2}{2}} dt}{\sigma} = 0,05 \quad \text{avec } h = 0,001,$$
soit $\frac{h\sqrt{N}}{\sigma} = 1,96$, et approximativement $N = 4.10^6 \sigma^2$.

Mais ce résultat n'a pas grande valeur en soi, puisqu'il dépend du seuil de probabilité accepté, qui est tout à fait imprécis.

Le seul résultat à retenir est que l'erreur décroît comme $\frac{1}{\sqrt{N}}$, c'est-à-dire lentement. La suite des $x_n$ est irrégulière et remplit progressivement l'intervalle (0.1). Mais on ne sait rien d'autre sur elle que son caractère aléatoire. Elle ne possède aucune organisation qui, sans modifier les principes, puisse favoriser la convergence.

Il est évidemment maladroit de calculer une intégrale par la formule (2.1). Les méthodes numériques usuelles sont plus efficaces. La "méthode des rectangles", moins précise que la méthode des traèzes ou les méthodes d'interpolation par polynômes, donne déjà une précision en $\frac{1}{N}$. La méthode statistique consiste à prendre des points au hasard, de sorte que leur répartition, médiocre pour les petites valeurs de N, en ce sens qu'il peut y avoir peu de points dans un intervalle, et beaucoup dans un autre, ne s'uniformise vraiment que lorsque N devient grand. Au contraire, la méthode des rectangles consiste, une fois N fixé, à placer N - 1 points suivant la meilleure répartition possible entre 0 et 1, c'est-à-dire de telle sorte que leurs abscisses soient en progression arithmétique. Une valeur de N relativement petite peut donc déjà donner une précision que la méthode probabiliste n'atteindrait qu'avec une valeur bien plus grande.

On notera cependant que l'évaluation de la précision dans la méthode des rectangles n'est possible que si la fonction $f(x)$ est suffisamment régulière. Si $f(x)$ admet une dérivée continue, telle que, sur l'intervalle (0.1), on ait $|f'| < M$, on démontre à l'aide d'une intégration par parties que

$$ \left| I - \frac{1}{N} \sum_{k=0}^{N-1} f\left(\frac{k}{N}\right) \right| < \frac{M}{2N} \quad (2.5) $$

Cette formule n'est donc pas applicable si $f(x)$ n'a pas une dérivée continue, ou est définie empiriquement. On est alors conduit à choisir N d'une façon plus subjective que rationnelle, et on doit évaluer la précision du calcul par des procédés empiriques. Si la précision semble insuffisante, on doit recommencer l'ensemble des opérations avec une valeur supérieure de N. Il est d'ailleurs difficile d'évaluer à vue la précision.
<figure><img src="image_2.png" /><figcaption>Figure 1</figcaption></figure>

Au contraire, la méthode statistique ne nécessite pas un choix initial sur N. Au fur et à mesure de la croissance de N, la valeur de S, se stabilise et on la voit tendre vers la limite I. On s'arrête quand on a atteint la précision souhaitée.

### III - NOMBRES ALEATOIRES -

Le premier problème qui se pose alors est de construire une suite de "nombres aléatoires". Il s'agit de nombres aléatoires au sens restreint, c'est-à-dire de nombres dont la suite ordonnée puisse être considérée comme le résultat d'une suite d'épreuves sur une variable aléatoire uniformément répartie entre 0 et 1. Je ne m'occuperai pas ici du problème général de la construction des nombres aléatoires ayant une répartition quelconque. Si X a pour fonction de répartition F(x), la variable aléatoire Y = F(x) est uniformément distribuée entre 0 et 1. En effet, puisque F(x) est non décroissante,

$$ \text{Prob} [X < x] = \text{Prob} [Y < F(x)] = F(x), $$

et si l'on pose F(x) = y, on a bien

$$ \text{Prob} [Y < y] = y $$

On peut donc toujours ramener une variable aléatoire quelconque à une variable uniformément distribuée entre 0 et 1 et le problème qui se pose est un problème purement pratique.
Pour construire une succession de valeurs numériques qui soient le résultat d'épreuves indépendantes sur une variable aléatoire uniformément répartie entre 0 et 1, on doit recourir à un jeu statistique convenablement choisi. On obtient ainsi ce qu'on appelle une table de nombres aléatoires. Il n'est d'ailleurs pas facile de réaliser un jeu statistique qui soit à la fois assez précis, assez rapide et assez garanti pour donner avec un nombre de décimales adapté aux machines à calculer modernes plusieurs dizaines de milliers de nombres successifs. C'est ce qui a été fait pour calculer les tables RAND, et les tests statistiques semblent prouver que le résultat est correct.

Cependant d'autres méthodes ont été proposées pour calculer des nombres aléatoires. Ces méthodes reposent sur des calculs de nature arithmétique dans lesquels le hasard ne joue aucun rôle. Ces calculs sont organisés de telle sorte que la suite $x_1, x_2, ...$ qu'ils fournissent soit compliquée et donne l'impression du hasard. Les nombres obtenus sont ensuite testés par des procédés statistiques.

Malheureusement, les méthodes arithmétiques sont souvent difficiles à étudier du point de vue mathématique. Comme les tests statistiques ne sont pas des preuves absolues de la validité des hypothèses, les méthodes arithmétiques généralement employées restent très incertaines. On ne peut pas affirmer qu'elles puissent jouer le même rôle que les méthodes de nature spécifiquement probabilistes.

Il existe cependant certains cas où l'arithmétique est une science assez avancée pour qu'on puisse l'invoquer en toute connaissance de cause lorsqu'il s'agit de construire des nombres aléatoires. C'est à l'examen de ces théories arithmétiques qu'est consacrée la suite de cette étude. On notera que, dès le moment où l'on utilise des théorèmes précis, les méthodes de calcul qu'on obtient n'ont plus rien de statistique, sauf l'apparence. Ce sont des méthodes mathématiques, présentant de simples analogies de forme et d'emploi avec les méthodes statistiques. On notera aussi qu'elles présentent un avantage, c'est que les "rombres aléatoires" arithmétiques n'ont pas besoin d'être rassemblés et prélevés dans des tables. Ils se calculent par des procédés compatibles avec la technique habituelle de programmation des machines à calculer.

## IV - THEOREMES DE H. WEYL(1)

Donnons d'abord une définition. Soit

$$x_1, x_2, ..., x_n, ...$$

-----------------

<footnote>(1) Référence (5).</footnote>
une suite de nombres compris entre 0 et 1. Considérons un intervalle $(\alpha, \beta)$ tel que $0 < \alpha < \beta < 1$. Supposons que, parmi les N premiers termes de la suite, il y en ait $N'$ dans l'intervalle $(\alpha, \beta)$. Si le rapport $\frac{N'}{N}$ tend, lorsque $N \to \infty$, vers une limite égale à $\beta - \alpha$, et cela quel que soit l'intervalle $(\alpha, \beta)$ choisi, on dit que la suite $x_n$ est _uniformément dense_, ou _équirépartie_, sur $(0, 1)$.

Les théorèmes de H. Weyl s'énoncent de la façon suivante :

### THEOREME 1 -

Pour que la suite $x_n$ soit équirépartie sur $(0, 1)$, il faut et il suffit que, pour toute fonction $f(x)$ intégrable au sens de Riemann sur $(0, 1)$, on ait

$$ \int_0^1 f(x) dx = \lim_{N \to \infty} \frac{1}{N} \sum_{n=0}^{N-1} f(x_n) $$

### THEOREME 2 -

Pour que la suite $x_n$ soit équirépartie sur $(0, 1)$, il faut et il suffit que, pour tout entier $l \neq 0$, on ait

$$ \lim_{N \to \infty} \frac{1}{N} \sum_{n=0}^{N-1} e^{2i\pi l x_n} = 0 $$

On définit d'une façon analogue une suite de points

$P_1, P_2, \dots, P_n, \dots$

_équirépartie_ dans l'hypercube fondamental $C$ de l'espace à dimensions $c$ est-à-dire l'hypercube dont les points $P$ ont des coordonnées toutes comprises entre 0 et 1. Soit $\emptyset$ un domaine quarrable intérieur à $C$, et de volume $V$. Supposons que, parmi les N premiers points $P_n$, il y en ait $N'$ à l'intérieur de $\emptyset$. Si le rapport $\frac{N'}{N}$ tend lorsque $N \to \infty$ vers une limite égale à $V$, et cela quel que soit le domaine $\emptyset$, on dit que la suite $P_n$ est équirépartie dans $C$. On démontre alors que :

### THEOREME 3 -

Pour que la suite $P_n$ soit équirépartie dans l'hypercube $C$, il faut et il suffit que, pour toute fonction $f(P)$ intégrable au sens de Riemann dans l'hypercube $C$, on ait

$$ \int_C f(P) d\omega = \lim_{N \to \infty} \frac{1}{N} \sum_{n=0}^{N-1} f(P_n) $$
où dw est l'élément de volume de l'espace.

Les théorèmes 1 et 3 permettent donc de calculer numériquement des intégrales simples et multiples. Le théorème 1 est l'équivalent arithmétique exact du théorème statistique du paragraphe 2. La seule différence, c'est que, en arithmétique, il y a convergence au sens ordinaire du second membre vers le premier, alors qu'en statistique il y a convergence en probabilité (ou convergence presque sûre), ce qui est moins précis.

La raison "physique" de validité de ces théorèmes est la même dans les deux cas. La suite $x_n$ tend à remplir "d'une façon homogène" le segment (0,1). La suite $P_n$ tend à remplir "régulièrement" l'hypercube C. La moyenne des valeurs de $f$ en des points régulièrement répartis est une approximation de la valeur moyenne théorique, c'est-à-dire de l'intégrale. Le sens du mot "remplir" ne doit naturellement pas être interprété d'une façon trop concrète. La suite $x_n$ est dense dans le segment (0,1). Il existe des points de la suite $x_n$ aussi voisins qu'on le désire de tout point entre 0 et 1. Mais cette suite est dénombrable et de mesure nulle. Elle ne recouvre pas (au sens mathématique du mot) le segment (0,1). Cependant, si l'on considère une suite de segments de longueur $\epsilon$ ayant pour centres les points $x_n$, elle recouvre effectivement le segment (0,1), et cela si petit que soit $\epsilon$.

Le problème du calcul numérique des intégrales par la méthode de Monte-Carlo est donc transformé en un problème arithmétique. Il s'agit maintenant de former des suites équiréparties, et d'étudier l'approximation qu'on doit attendre en remplaçant

$$ \int_C f(P) d\omega \text{ par } \frac{1}{N} \sum_{n=0}^{N-1} f(P_n) $$

V - SUITES EQUIREPARTIES - ETUDE DE LA SUITE An -

Le théorème 2 du paragraphe 4 permet de démontrer facilement que.

THEOREME 4 -

Si A est un nombre irrationnel, la suite An des multiples de $\Lambda$ est équirépartie modulo 1.

Cela veut dire que la suite $x_n$ des parties fractionnaires de An est équirépartie sur (0,1).

En effet,
$$ \frac{1}{N} \sum_{n=0}^{N-1} e^{2i\pi lAn} = \frac{1}{N} \frac{e^{2i\pi lAn} - 1}{e^{2i\pi lAn} - 1} \quad (5.1) $$

Si l est un entier non nul, 1A est irrationnel et $e^{2i\pi lA} - 1$ n'est pas nul. Le numérateur est majoré par 2. Le second membre tend donc vers 0 avec $\frac{1}{N}$.

Pour représenter la suite An modulo 1, il est commode de considérer $2\pi An$ comme l'abscisse curviligne d'un point sur une circonférence de rayon 1. Si A est irrationnel, on démontre facilement que les points ayant pour abscisses curvilignes $2\pi A$, $4\pi A$, $6\pi A$,... sont denses sur le cercle. Ce qui est remarquable et moins élémentaire, c'est qu'ils le soient _uniformément_, au sens indiqué ci-dessus.

<figure><img src="image_3.png" /><figcaption>Figure 2</figcaption></figure>

La somme

$$ \sum_{n=0}^{N-1} e^{2i\pi lAn} $$

se construit d'une façon analogue. C'est une somme de nombres complexes de module 1 qui, dans le plan complexe, sont représentés par des vecteurs unitaires tels que chacun fasse avec le précédent un angle constant égal à $2\pi lA$. Ces vecteurs sont les côtés d'un polygône régulier non fermé, inscrit dans un cercle de rayon $\frac{1}{2 \sin \pi lA}$.

<figure><img src="image_4.png" /><figcaption>Figure 3</figcaption></figure>

On notera que An ne s'annule jamais mod. 1, sauf pour $n = 0$, mais qu'il existe des valeurs de n pour lesquelles An est aussi proche qu'on le désire de 0, mod. 1.
## VI - SUITES DE H. WEYL - CORRELATION DANS LES SUITES ARITHMETIQUES -

De la suite An, on peut déduire par récurrence d'autres suites équiréparties. H. Weyl a démontré le théorème suivant :

### THEOREME 5 -

$\text{Si } \varphi(t) = A_1t^\nu + A_2t^{\nu-1} + \dots + A_\nu \text{ est un polynôme de degré } v \ge 1, \\
\text{dont au moins un coefficient, autre que } A_\nu, \text{ est irrationnel, la suite } \\
\varphi(n) \text{ est équirépartie modulo } 1.$

On en trouvera dans [1] une démonstration utilisant la notion de fonction de corrélation.

Ce théorème met donc à la disposition du calculateur un grand nombre de suites qui peuvent jouer le rôle de suites de nombres aléatoires. Le problème qui se pose est de savoir quelles sont parmi ces suites les plus favorables au calcul. Quelle est l'influence du degré $v$ du polynôme et de la nature des nombres irrationnels ? Quelle est d'autre part l'influence du fait que, dans la pratique, on remplace obligatoirement tout nombre irrationnel par une approximation rationnelle ?

Avant de répondre à ces questions, nous allons examiner un peu plus en détail la nature des suites équiréparties qui sont fournies par les théorèmes de H. Weyl.

Nous aurons besoin de nous servir d'un lemme, qui est dû à H. Weyl, et qui constitue le

### THEOREME 6 -

Pour que la suite de points $P_n$, ayant pour coordonnées $x_n, y_n, z_n, \dots$ soit équirépartie dans l'hypercube $C$, il faut et il suffit que la suite

$$l_1x_n + l_2y_n + l_3z_n + \dots$$

soit équirépartie mod. 1 sur le segment $(0.1)$, quels que soient les en-
tiers $l_1, l_2, l_3, \dots$, non tous nuls.

Nous allons déduire du théorème 6 que

### THEOREME 7 -

$\text{Si } x_n = \varphi(n) \bmod 1 \text{ est une suite satisfaisant aux hypothèses du théorème 5, avec } v \ge 2, \text{ la corrélation entre } x_n \text{ et } x_{n+h} \text{ est nulle pour tout } h \text{ entier } > 1.$
## Démonstration.

Montrons d'abord que la suite $x_n$, $x_{n+h}$ est équirépartie dans le carré C. D'après le théorème 6, il suffit de vérifier que la suite

$$l_1 x_n + l_2 x_{n+h}$$

est équirépartie mod. 1 pour tout couple d'entiers $l_1, l_2$ non simultanément nuls. Celà résulte de ce que $l_1 \varphi(t) + l_2 \varphi(t+h)$ est un polynôme de degré $v$ en général, de degré $v-1$ si $l_1 + l_2 = 0$, satisfaisant aux hypothèses du théorème 5.

D'après le théorème 3, on a alors

$$\lim_{N \to \infty} \frac{1}{N} \sum_{n=0}^{N-1} x_n x_{n+h} = \iint_C xy \, dx \, dy = \frac{1}{4}$$

Comme $x_n$ et $x_{n+h}$ ont pour moyenne $\frac{1}{2}$, la moyenne du produit

$$(x_n - \frac{1}{2})(x_{n+h} - \frac{1}{2})$$

est nulle. Le coefficient de corrélation (envisagé au sens des moyennes arithmétiques) entre $x_n$ et $x_{n+h}$ est nul. Nous dirons que la suite $x_n$ ne présente pas de corrélation interne.

## Remarque.

Ce résultat n'est pas exact dans le cas de la suite An. Désignons par $x_n$ la partie fractionnaire de An, par $y_n$ celle de $A(n + h)$, par $R$ celle de Ah. On a

$$y_n = x_n + R \quad \text{si} \quad x_n + R < 1,$$

$$y_n = x_n + R - 1 \quad \text{si} \quad x_n + R > 1.$$

$y_n$ est donc relié à $x_n$ par une relation indépendante de n. Le point $(x, y)$ se trouve sur l'un des deux segments dessinés sur la figure 5. La suite $(x_n, y_n)$ n'est pas dense dans le carré C.

On a d'ailleurs

$$\begin{align*} \frac{1}{N} \sum_{n=0}^{N-1} x_n y_n &= \frac{1}{N} \sum_{x_n < 1-R} x_n (x_n + R) + \frac{1}{N} \sum_{x_n > 1-R} x_n (x_n + R - 1) \\ &= \frac{1}{N} \sum_{n=0}^{N-1} x_n (x_n + R) - \frac{1}{N} \sum_{x_n > 1-R} x_n. \end{align*}$$
<figure><img src="image_5.png" /><figcaption>Figure 5</figcaption></figure>

La suite $x_n$ étant équirépartie, la première somme est égale à

$$ \int_0^{\frac{1}{R}} x(x+R)dx = \frac{1}{3} + \frac{R}{2} $$

La seconde est égale à l'intégrale de la fonction $f(x)$ définie par

$$ \begin{aligned} f(x) &= 0 \text{ si } x < 1 - R \\ &= x \text{ si } x > 1 - R \end{aligned} $$

Elle vaut

$$ \int_{1-R}^{\frac{1}{R}} xdx = R - \frac{R^2}{2} $$

$x_n y_n$ a donc pour moyenne

$$ \frac{1}{3} - \frac{R}{2} + \frac{R^2}{2} $$

Comme $x_n$ et $y_n$ ont pour moyenne $\frac{1}{2}$, pour moyenne quadratique $\frac{1}{3}$, le coefficient de corrélation entre $x_n$ et $y_n$ a pour valeur :

$$ \frac{\frac{1}{3} - \frac{R}{2} + \frac{R^2}{2} - \frac{1}{4}}{\frac{1}{3}} = \frac{1}{4} - \frac{3}{2} R(1-R) $$

Il est compris entre $-\frac{1}{8}$ et $+\frac{1}{4}$, et ne s'annule que pour des valeurs exceptionnelles de h.

## VII - SUITES A DEUX VALEURS -

En vue des applications numériques, (problèmes de cheminements discontinus), il est important de connaître des suites arithmétiques à deux valeurs, jouant le rôle de suites de nombres aléatoires à deux valeurs également probables. Ces suites ont autant d'applications que les suites équiréparties.

En arithmétique comme en statistique, les suites équiréparties permettent de former des suites à deux valeurs. En statistique, si X
est une variable aléatoire équirépartie sur (0, 1), la valeur aléatoire Y égale à 0 si $0 < X < \frac{1}{2}$, à 1 si $\frac{1}{2} < X < 1$, prend bien deux valeurs équiprobables

En arithmétique, on démontre le théorème suivant :

## THEOREME 8 -

Soit $\varphi(t)$ un polynôme de degré $v \ge 2$ ayant un coefficient irrationnel (autre que le terme constant). La suite des parties entières de $\varphi(n)$ prend modulo 2 les valeurs 0 et 1 avec les fréquences $\frac{1}{2}$ et $-\frac{1}{2}$. Elle n'a pas de corrélation interne.

Désignons en effet par $x_n, y_n$ les parties fractionnaires de $\frac{1}{2}\varphi(n)$ et de $\frac{1}{2}\varphi(n + h)$, où h est un entier positif.

On peut démontrer (cf. [1], p.26) que la suite des points $(x_n, y_n)$ est équirépartie dans le carré C du plan. Par homothétie, on en déduit que la suite $(2x_n, 2y_n)$ est équirépartie dans le carré C' défini par

$$0 < x < 2, \quad 0 < y < 2.$$

Désignons par $\xi_n, \eta_n$ les parties entières de $2x_n, 2y_n$ prises modulo 2. $\xi_n$ est égal à 0 ou à 1 suivant que la partie entière de $2x_n$ est paire ou impaire. Le point $(\xi_n, \eta_n)$ peut occuper quatre positions, qui sont les sommets d'un carré et ont pour coordonnées

$$(0.0), \quad (0.1), \quad (1.0), \quad (1.1)$$

Les droites $x = 1$ et $y = 1$ séparent le carré C' (figure 6) en quatre carrés C, $C_1$, $C_2$, $C_3$.

Lorsque le point $2x_n, 2y_n$ est par exemple dans le carré $C_1$, le point $(\xi_n, \eta_n)$ occupe la position (1.0). Or la "probabilité" (au sens de fréquence arithmétique) pour que le point $(2x_n, 2y_n)$ soit dans $C_1$ est $1/4$. C'est aussi la "probabilité" pour que le point $(\xi_n, \eta_n)$ se trouve en (1.0). Le même raisonnement s'appliquant aux trois autres positions, on voit que les quatre positions du point $(\xi_n, \eta_n)$ ont la même fréquence $1/4$. On voit

<figure><img src="image_6.png" /><figcaption>Figure 6</figcaption></figure>
d'autre part sans difficulté que $\xi_n$, partie entière de $2x_n$, est aussi la partie entière de $\varphi(n)$.

Le produit $\xi_n \eta_n = \xi_n \xi_{n+h}$ prend, avec des fréquences égales, les quatre valeurs 0, 0, 0, 1. Sa moyenne vaut donc $1/4$. Comme $\xi_n$ et $\eta_n$ ont pour moyenne $1/2$, le coefficient de corrélation entre $\xi_n$ et $\xi_{n+h}$ est nul, ce qui démontre le théorème.

## VIII - APPROXIMATIONS DANS LE CALCUL D'UNE INTEGRALE SIMPLE -

Nous nous proposons ici de chercher une limite supérieure de l'erreur commise en remplaçant l'intégrale

$$ \int_0^1 f(x) dx \qquad (8.1) $$

par

$$ \frac{1}{N} \sum_{n=0}^{N-1} f(x_n), \qquad (8.2) $$

où $x_n$ est la partie fractionnaire de An, et où A est un nombre positif donné. Nous savons que l'erreur tend vers 0 lorsque A est irrationnel. Cependant, dans la pratique, on doit naturellement remplacer le nombre irrationnel par une fraction irréductible $\frac{p}{q}$. Nous examinerons d'abord le cas théorique A irrationnel et ensuite le cas pratique $A = \frac{p}{q}$.

Nous aurons besoin dans la démonstration de résultats élémentaires concernant les fractions continues(1). Soit A un nombre positif non entier, $\hat{A}$ sa partie entière. On pose

$$ A = \hat{A} + \frac{1}{a_1} \qquad a_1 > 1 $$

$$ a_1 = \hat{a}_1 + \frac{1}{a_2} \qquad a_2 > 1 $$

$$ a_2 = \hat{a}_2 + \frac{1}{a_3} \qquad a_3 > 1 $$

etc.

On forme par ce procédé une suite $a_1, a_2, a_3, ...$ de nombres

<footnote>(1) Consulter à ce sujet: G. Valiron, Théorie des fonctions (Masson), p. 17.</footnote>
supérieurs à 1 dont les parties entières sont appelées *quotients in-complets*. Si A est rationnel, il arrive au bout d'un nombre fini d'opérations que $a_m$ devienne entier. Si A est irrationnel, $a_m$ n'est jamais entier, et le calcul se poursuit indéfiniment. Si l'on remplace $a_m$ par sa partie entière, on obtient une expression approchée de A, par des opérations rationnelles sur des entiers.

Cette représentation approchée est une fraction irréductible $\frac{P_m}{Q_m}$, appelée *réduite*. On démontre les résultats suivants :

1/ Lorsque $m \to \infty$, $P_m$ et $Q_m$ augmentent indéfiniment.

2/ $\frac{P_m}{Q_m}$ tend vers une limite égale à A.

$$3/ \left| A - \frac{P_m}{Q_m} \right| < \frac{1}{Q_m Q_{m+1}}$$

<figure><img src="image_7.png" /><figcaption>Forme primitive de $f(x)$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Forme modifiée de $f(x)$</figcaption></figure>

Figure 7

4/ La suite $Q_m$ est croissante et $Q_{m+2} > 2Q_m$. Si A est rationnel, la suite $\frac{P_m}{Q_m}$ a un nombre fini de termes. A partir d'un certain rang, $\frac{P_m}{Q_m}$ reste égal à A.

Ceci posé, nous allons d'abord faire l'hypothèse que $f(x)$ est une *fonction continue*, et qu'elle est représentable entre 0 et 1 par une série de Fourier absolument convergente. Remarquons que ce qui est continu, c'est non pas seulement $f(x)$, fonction définie sur $(0,1)$, mais bien la fonction périodique engendrée par $f(x)$. Il en résulte que nécessairement $f(0) = f(1)$. Cette condition n'est pas toujours vérifiée dans la pratique. Il faut alors faire subir à $f(x)$ une modification qui la rende continue sans que le calcul de l'intégrale s'en trouve perturbé. On représente
par exemple $f(x)$ à une échelle des abscisses plus petite, c'est-à-dire
sur un intervalle $[0, a]$, $a < 1$, et l'on prolonge $f(x)$ par une fonction
linéaire, comme l'indique la figure 7.

On considère donc la fonction $f_1(x)$ qui est égale à $f(\frac{x}{a})$ si $0 < x < a$, et qui varie linéairement de $f(1)$ à $f(0)$ lorsque $x$ varie de $a$ à 1. On a

$$ \int_0^1 f_1(x)dx = \int_0^a f(\frac{x}{a})dx + \frac{1-a}{2}[f(1)-f(0)] = a\int_0^1 f(x)dx + \frac{1-a}{2}[f(1)-f(0)] $$

Le calcul numérique de l'intégrale de $f_1(x)$ fournit la valeur de l'intégrale cherchée.

On peut aussi doubler $f(x)$ par symétrie, c'est-à-dire introduire
(figure 8) sur l'intervalle $(0, 1)$, la fonction $g(x)$ égale à

$$ \begin{cases} f(2x) & \text{si } 0 < x < \frac{1}{2} \\ f(2 - 2x) & \text{si } \frac{1}{2} < x < 1 \end{cases} $$

On a

$$ \begin{aligned} \int_0^1 g(x)dx &= 2 \int_0^{\frac{1}{2}} f(2x)dx = \\ &= \int_0^1 f(x)dx \end{aligned} $$

On remarque que, si $f(1) \neq f(0)$,
ce qu'on a supposé, $f(x)$, considérée comme fonction périodique, n'est pas une fonction paire, alors que $g(x)$ est nécessairement paire. Son développement en série de Fourier n'a que des termes en cosinus, dont seuls subsistent ceux dont le rang est pair.

Supposons d'abord que $f(x)$ soit un polynôme trigonométrique ayant un nombre fini de termes :

$$ \begin{gathered} f(x) = \sum_{k=-r}^{r} c_k e^{i\pi kx}, \\ c_{-k} = \bar{c}_k \text{ (conjugué de } c_k \text{).} \end{gathered} $$

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8</figcaption></figure>
Comme $\int_0^1 f(x)dx = c_0$, l'erreur commise dans le calcul de l'intégrale vaut

$$ 
\begin{align*} 
I_N &= \frac{1}{N} \sum_{n=0}^{N-1} \sum_{k \neq 0} c_k e^{2i\pi k \Lambda n} \\ 
&= \frac{1}{N} \sum_{k \neq 0} c_k \sum_{n=0}^{N-1} e^{2i\pi k \Lambda n} \\ 
&= \frac{1}{N} \sum_{k \neq 0} c_k \frac{e^{2i\pi k \Lambda N} - 1}{e^{2i\pi k \Lambda} - 1} 
\end{align*} 
$$

La somme qu'on doit diviser par N comporte un nombre fini de termes bornés. L'erreur tend donc vers 0 comme $\frac{1}{N}$. Si l'on désigne par $c_0$ celui des entiers 1, 2, ..., r pour lequel $e^{2i\pi k\Lambda} - 1$ a le plus petit module, on a

$$ |I_N| < \frac{4}{N} \frac{\sum_{k=1}^{r} |c_k|}{|e^{2i\pi k_0 \Lambda} - 1|} \quad (8.3) $$

## THEOREME 9 -

L'erreur commise en appliquant au calcul de l'intégrale d'un polynôme trigonométrique la méthode de Monte-Carlo, sous la forme arithmétique que justifie le théorème de H. Weyl, avec la suite $x_n = \lfloor n \Lambda \rfloor \mod 1$ (A irrationnel) tend vers 0 comme $\frac{1}{N}$.

Le résultat est donc meilleur qu'avec des nombres aléatoires véritables, car alors l'approximation est en $\frac{1}{\sqrt{N}}$, et non en $\frac{1}{N}$. Mais en pratique le coefficient de $\frac{1}{N}$ dans la formule (8.1) risque d'être grand, à cause du dénominateur. Il faut donc dans les deux cas beaucoup de termes (grande valeur de N) pour obtenir une bonne précision.

Si maintenant la série qui représente $f(x)$ est infinie, on a

$$ f(x) = \sum_{k=0}^{\infty} c_k e^{2i\pi kx}, \quad c_{-k} = \bar{c}_k. \quad (8.4) $$

On va supposer qu'il existe deux nombres positifs M et $\lambda$ tels que

$$ |c_k| < \frac{M}{k^{1+\lambda}} \quad (8.5) $$

La série est alors absolument convergente. On sait que, si $f(x)$ est
continue et si $f'(x)$ existe et est continue, sauf en un nombre fini de points où $f'(x)$ a des sauts, alors $\lambda = 1$. Si $f'$ est partout continue et si $f''$ existe et a les propriétés précédemment requises pour $f'$, $\lambda = 2$, etc. le reste de la série $|c_k|$ est majoré, pour $k \neq 0$, par

$$R_k = 2M \sum_{k=k+1}^{\infty} \frac{1}{k^{1+\lambda}} < \frac{2M}{\lambda k^{\lambda}} \quad (8.6)$$

La différence entre l'intégrale et sa valeur approchée est égale à

<figure><img src="image_11.png" /><figcaption>Figure 9</figcaption></figure>

$$ 
\begin{aligned}
I_N &= \sum_{k \neq 0} c_k \frac{1}{N} \sum_{n=0}^{N-1} e^{2i\pi k An} \\
&= \sum_{k \neq 0} c_k \frac{1}{N} \frac{e^{2i\pi k AN} - 1}{e^{2i\pi k A} - 1} \quad (8.7)
\end{aligned}
$$

Les dénominateurs $e^{2i\pi kA} - 1$ ne sont jamais nuls. Mais nous saisons que, lorsque $k \to \infty$, leur module prend des valeurs aussi voisines de 0 qu'on le désire. Pour majorer $I_N$, séparons les termes de rang supérieur à Q, et les termes de rang inférieur.

Les premiers ont une somme inférieure en module à

$$2 \sum_{k=Q+1}^{\infty} |c_k| \quad (8.8)$$

car

$$\frac{1}{N} \sum_{n=0}^{N-1} e^{2i\pi k n}$$

est en module inférieur à 1.

Sur ces termes, on commet donc une erreur inférieure à

$$R_0 = \frac{2M}{\lambda Q^{\lambda}} \quad (8.9)$$

Désignons maintenant par $\eta$ le minimum de $|e^{2i\pi kA} - 1|$ lorsque $k$ prend les valeurs 1, 2, ..., Q. $\eta$ est une fonction de Q. Pour évaluer $\eta$, utilisons la représentation approchée de A par des réduites $\frac{P_m}{Q_m}$. On peut toujours mettre A sous la forme

$$A = \frac{P_m}{Q_m} + \frac{\vartheta_m}{Q_m Q_{m+1}}, \quad |1\vartheta_m| < 1. \quad (8.10)$$
Si l'on néglige le terme en $\vartheta_m$, la suite des nombres complexes

$$e^{\frac{2 i \pi k \frac{P_m}{Q_m}}{Q_m}}$$

devient périodique. Si k parcourt l'ensemble des valeurs

$1, 2, \dots, Q_{m-1},$

la plus petite valeur de $\left| e^{\frac{2 i \pi k \frac{P_m}{Q_m}}{Q_m}} - 1 \right|$ est égale à $\left| e^{\frac{2 i \pi}{Q_m}} - 1 \right|$ (8.11)

On sera sûr de ne pas dépasser cette valeur si l'on choisit pour Q un nombre inférieur ou égal à $Q_{m-1}$. Choisissons $Q=Q_{m-1}$. Comparons alors la suite des points

$$e^{\frac{2 i \pi k \frac{P_m}{Q_m}}{Q_m}} \quad \text{et} \quad e^{2 i \pi k A}$$

sur le cercle trigonométrique. On a

$$kA = \frac{kP_m}{Q_m} + \frac{k\vartheta_m}{Q_m Q_{m+1}}.$$

Or $\left| \frac{k \vartheta_m}{Q_m} \right|$ est toujours inférieur à $\frac{k}{Q_m} < \frac{Q_{m-1}}{Q_m}$.

Sur le cercle trigonométrique, les arcs $\frac{kP_m}{Q_m}$ et $kA$ (au facteur $2\pi$ près) diffèrent de moins de $\frac{Q_{m-1}}{Q_m Q_{m+1}}$. Lorsque $m \to \infty$, $Q_m$ augmente indéfiniment, au moins aussi vite que $2^{\frac{2}{2}}$, et $\frac{Q_{m-1}}{Q_{m+1}}$ reste inférieur à $\frac{1}{2}$.

Donc $\frac{Q_{m-1}}{Q_m} < \frac{Q_{m+1}}{Q_m}$ tend vers 0 et les points

$$e^{\frac{2 i \pi k \frac{P_m}{Q_m}}{Q_m}}, e^{2 i \pi k A}$$

sont très voisins. L'argument du terme de la suite $e^{2 i \pi k A}$ ($k=1, 2, \dots, Q_{m-1}$) qui est le plus voisin de 1 est au moins égal à

$$2\pi \left( \frac{1}{Q_m} - \frac{Q_{m-1}}{Q_m Q_{m+1}} \right) = \frac{2\pi}{Q_m} \left( 1 - \frac{Q_{m-1}}{Q_{m+1}} \right)$$

Finalement, le module de $e^{2 i \pi k A} - 1$ est au moins égal à
$$ \left| e^{\frac{2i\pi}{Q_m}\left(1 - \frac{Q_{m-1}}{Q_{m+1}}\right)} - 1 \right| = 2 \sin \frac{\pi}{Q_m} \left( 1 - \frac{Q_{m-1}}{Q_{m+1}} \right) > \frac{4}{Q_m} \left( 1 - \frac{Q_{m-1}}{Q_{m+1}} \right) $$

Comme $Q_{m+1} > 2Q_{m-1}$, cette quantité est elle-même supérieure à

$$ \frac{4}{Q_m} \left(1 - \frac{1}{2}\right) = \frac{2}{Q_m} \qquad (8.12) $$

C'est cette valeur que nous allons attribuer à la quantité $\eta$.

Au total, l'erreur commise sur l'intégrale est majorée par

$$ \frac{2M}{\lambda Q_{m-1}^{\lambda}} + \frac{4}{N\eta} \sum_{k=1}^{Q_{m-1}} |c_k| $$

Soit S la somme de la série $|c_1| + ... + |c_k| + ...$ D'une façon générale, S est inférieure à $c_1 + \frac{M}{\lambda} < \frac{M(\lambda + 1)}{\lambda}$, mais cette majoration n'est pas très utile, et peut en général être remplacée par une majoration plus précise. Si par exemple $\lambda = 1$, on a $|c_k| < \frac{M}{k^2}$, et $S < M \frac{\pi^2}{6}$. Quoi qu'il en soit, on obtient

$$ |I_N| < \frac{2M}{\lambda Q_{m-1}^{\lambda}} + \frac{4S}{\eta N} $$

c'est-à-dire

$$ |I_N| < \alpha_m + \frac{\beta_m}{N} \qquad (8.13) $$

avec

$$ \alpha_m = \frac{2M}{\lambda} \frac{1}{Q_{m-1}^{\lambda}}, \quad \beta_m = 2S Q_m. $$

Nous allons interpréter la formule (8.13).

$\lambda$ était positif, la suite $\alpha_m$ décroît et tend vers 0. Pour une valeur fixée de m, l'erreur $|I_N|$ décroît linéairement en fonction de $\frac{1}{N}$ et tend, lorsque $N \to \infty$, vers une limite $\alpha_m$ d'autant plus petite que m est plus grand.

Etudions ce qui se passe lorsqu'on remplace m par m + 1.

Dans le plan des $(u, N)$, on remplace l'hyperbole $H_m$ qui a pour équation
$$u = \alpha_m + \frac{\beta_m}{N}$$

par l'hyperbole $H_{m+1}$ qui a pour équation

$$u = \alpha_{m+1} + \frac{\beta_{m+1}}{N}$$

Or $\beta_{m+1} = 2SQ_{m+1}$ est supérieur à $\beta_m$ et augmente indéfiniment avec m, pendant que $\alpha_m$ décroît et tend vers 0. Les hyperboles $H_m$ et $H_{m+1}$ se coupent en un point qui a pour abscisse le nombre positif.

$$N_m = \frac{\beta_{m+1} - \beta_m}{\alpha_m - \alpha_{m+1}}$$

La disposition relative de $H_m$ et $H_{m+1}$ est celle qu'indique la figure 10 :

<figure><img src="image_12.png" /><figcaption>Figure 10</figcaption></figure>
Pour limiter l'erreur $|I_N|$, on procède donc de la façon suivante : on choisit le nombre A. On calcule les $Q_m$ jusqu'à un rang $m_0$ tel que $\alpha_{m_0}$ soit suffisamment petit. On majore $|I_N|$ par

$$ \alpha_{m_0} + \frac{\beta_{m_0}}{N} $$

Pour améliorer la précision, on construit, à partir d'un certain rang $m_1 < m_0$, les hyperboles $H_{m_1}, H_{m_1+1}, \dots, H_{m_0}$. On conserve de cha-cune d'entre elles l'arc qui se trouve en-dessous de l'hyperbole précé-dente et au-dessus de l'hyperbole suivante. La courbe qui représente en fonction de N la limite supérieure de l'erreur commise est donc consti-tuée par une suite continue et décroissante d'arcs d'hyperboles, tendant vers 0 lorsque $N \to \infty$, et qu'il est facile de construire si les $Q_m$ ont été calculés. Si l'on choisit pour A le nombre $\pi$, qui n'est peut-être pas le nombre irrationnel donnant la convergence la plus rapide, on a

$$ \begin{array}{lclcl} Q_0 & = & 1 & Q_1 & = 7 & Q_2 & = 106 & Q_3 & = 113 \\ Q_4 & = & 33102 & Q_5 & = 33215 & Q_6 & = 66317\dots \end{array} $$

Il n'est pas sûr que cette méthode donne la meilleure évaluation possible de l'erreur commise. Elle est probablement assez grossière, car elle consiste à remplacer tous les dénominateurs $e^{2i\pi kA}$ - 1 par le plus petit d'entre eux, alors que certains ne sont pas très petits. Elle a l'avantage sur les méthodes dont il sera question plus loin de ne partir d'aucune hypothèse sur le nombre irrationnel A.

L'influence de la nature de l'irrationnel A peut s'interpréter de la façon suivante. Il y a intérêt à rester le moins longtemps possible sur une hyperbole $H_m$ fixée, et à passer à une hyperbole plus rappro-chée de l'axe des N. Il y a donc intérêt à ce que la suite $\beta_m$ soit régulièrement croissante, et cela de telle sorte qu'il faille changer le plus fréquemment possible d'hyperbole. Il semble que cette circonstance se produise lorsque A est un nombre algébrique, et plus particulièrement un nombre quadratique.

Si par exemple

$$ A = \frac{1 + \sqrt{5}}{2}, $$

on a

$$ A^2 = A + 1, \quad A = 1 + \frac{1}{A}. $$

Les quotients incomplets sont tous égaux à 1, et $Q_m$ est solution de l'é-quation aux différences finies
$$Q_m = Q_{m-1} + Q_{m-2},$$

avec

$$Q_0 = 1, \quad Q_1 = 1,$$

On trouve

$$Q_m = \frac{1}{\sqrt{5}} \left[ \left( \frac{1+\sqrt{5}}{2} \right)^{m+1} - \left( \frac{1-\sqrt{5}}{2} \right)^{m+1} \right], \quad (8.14)$$

$$Q_2 = 2, \quad Q_3 = 3, \quad Q_4 = \frac{135}{32} = 4, 21, \dots,$$

Lorsque $m \to \infty$, $(\frac{1-\sqrt{5}}{2})^{m+1}$ tend vers 0. On a donc approximativement

$$Q_m = \frac{1}{\sqrt{5}} A^{m+1}. \quad (8.15)$$

Il existe deux constantes $C$ et $C'$ telles que, si $m$ est suffisamment grand, on ait

$$|I_N| < \frac{C}{A^{m\lambda}} + C' \frac{A^n}{N}. \quad (8.16)$$

On vérifie que $H_m$ et $H_{m+1}$ se coupent en un point d'abscisse positive égale à

$$\frac{C'}{C} \frac{A - 1}{A^\lambda - 1} A^{\lambda(m+1)},$$

dont l'abscisse augmente avec $m$. La disposition est celle de la figure 11.

## IX - APPROXIMATION DANS LE CAS RATIONNEL -

Si l'on remplace $A$ par une fraction irréductible $\frac{P}{q}$, la suite des réduites $\frac{P_m}{Q_m}$ a un nombre fini de termes. A partir d'un certain rang, les termes ne varient plus. Tout ce qui a été dit au paragraphe 8 reste valable, sauf que la suite des hyperboles $H_m$ est limitée. Il y a une dernière asymptote qu'il n'est pas possible de dépasser. Cela ne prouve d'ailleurs pas que l'erreur commise ait une borne inférieure infranchissable. Cela indique seulement que, par la méthode employée, on ne peut plus réduire indéfiniment cette erreur en donnant à $N$ des valeurs de plus en plus grandes. Ce dernier résultat est cependant exact, et nous allons maintenant expliquer pourquoi.

L'erreur commise est égale à
<figure><img src="image_13.png" /><figcaption>Figure 11</figcaption></figure>

$$I_N = \sum_{k \neq 0} c_k \frac{1}{N} \sum_{n=0}^{N-1} e^{2i\pi k \frac{p}{q} n} \quad (9.1)$$

Si $\frac{p}{q}$ est une bonne approximation de A, $q$ est un nombre très grand, et N, bien qu'atteignant des valeurs très grandes, reste inférieur à $q$, et même petit devant $q$. Si le nombre $k \frac{p}{q}$ est entier, cela provient de $k$, et non de $n$. Distinguons donc deux cas :

1. $k$ est multiple de $q$.

$e^{2i\pi k \frac{p}{q} n}$ est alors égal à 1. La somme $\frac{1}{N} \sum_{n=0}^{N-1}$ vaut 1. Le terme correspondant de $I_N$ est exactement égal à $c_k$. Dans le calcul de $I_N$, on trouve donc la somme
$$c_q + c_{2q} + c_{3q} + \dots \qquad (9.2)$$

qui est indépendante de N et irréductible.

2/ *k* n'est pas multiple de *q*.

Le coefficient de $c_k$ est égal à

$$\frac{1}{N} \frac{e^{\frac{2i\pi k}{q} \frac{\omega}{N}} - 1}{e^{\frac{2i\pi k}{q}} - 1} \qquad (9.3)$$

Il est majoré par

$$\frac{1}{N} \frac{2}{\left| e^{\frac{2i\pi k}{q}} - 1 \right|} \qquad (9.4)$$

Il tend donc vers 0 comme $\frac{1}{N}$. On retrouve les résultats généraux de la discussion du paragraphe précédent. Mais on obtient ici une expression de la limite de $I_N$ lorsque $N \to \infty$. On est sûr que la valeur absolue de l'erreur commise est au moins égale à

$$c_q + c_{2q} + c_{3q} + \dots \qquad (9.5)$$

Si *q* est grand, cette quantité est naturellement très petite.

## X - RESULTATS DE RICHTMYER [4] -

Richtmyer semble avoir été le premier à signaler, en 1952, l'intérêt éventuel de l'utilisation de la suite An dans les calculs d'intégrales par la méthode de Monte-Carlo. Il a pu obtenir des majorations intéressantes de l'erreur $I_N$ lorsque *A* est un nombre algébrique. Voici quelques-uns de ses résultats (qu'il a réussi à généraliser aux intégrales multiples).

Rappelons d'abord deux théorèmes concernant les nombres algébriques.

Le théorème de Liouville est relatif à l'approximation des nombres algébriques par les nombres rationnels. Supposons que *A* soit racine d'une équation algébrique à coefficients entiers de degré *r* (et ne soit pas racine d'une équation de degré moindre). Soit $\frac{\rho}{q}$ une fraction ir-
réductible quelconque telle que $|A - \frac{p}{q}| < 1, q > 0$. Il existe alors une constante $c$, indépendante de $p$ et $q$, telle que

$$|qA - p| > \frac{c}{q^{r-1}}, \quad (10.1)$$

sauf peut-être pour un nombre fini de valeurs de $p$ et $q$ (l'inégalité est donc sûrement valable dès que $q$ est suffisamment grand).

Le polynôme dont $A$ est racine se met en effet sous la forme

$$g(x) = H(x - A)(x - A_1) \dots (x - A_r), \quad \text{H entier}$$

Donc

$$g(\frac{p}{q}) = H(\frac{p}{q} - A) (\frac{p}{q} - A_1) \dots (\frac{p}{q} - A_r)$$

Or

$$|\frac{p}{q} - A_i| = |\frac{p}{q} - A + A - A_i| < |\frac{p}{q} - A| + |A - A_i| < 1 + |A - A_i|$$

Les termes $\frac{p}{q} - A_i$ sont donc majorés par des quantités indépendantes de $\frac{p}{q}$. Il existe un nombre $c$, dépendant seulement du polynôme $g$, tel que

$$|g(\frac{p}{q})| < \frac{1}{c} |\frac{p}{q} - A|$$

D'autre part,

$$g(\frac{p}{q}) = A \frac{p^r}{q^r} + \dots = \frac{\text{entier}}{q^r}.$$

L'entier du numérateur a un module au moins égal à 1, sauf si $\frac{p}{q}$ est racine du polynôme $g$. Ce cas excepté, on a

$$|g(\frac{p}{q})| > \frac{1}{q^r}$$

Finalement

$$\frac{1}{c} |\frac{p}{q} - A| > \frac{1}{q^r}$$

D'où l'inégalité cherchée.
Dans le cas des nombres quadratiques, r est égal à 2, et le théorème de Liouville s'écrit

$$|qA - p| > \frac{c}{q}. \qquad (10.2)$$

Sauf dans ce dernier cas, on peut beaucoup améliorer le théorème de Liouville. D'après le théorème de Roth, l'inégalité

$$|qA - p| > \frac{1}{q^{\delta+1}} \qquad (10.3)$$

est valable quel que soit le nombre $\delta > 0$, aussi petit qu'on le désire, sauf peut-être pour un nombre fini de valeurs de p et q. On trouvera la démonstration de ce théorème dans le livre de Cassels [3].

Ceci posé, la méthode de Richtmyer consiste à partir de l'expression

$$I_N = \sum_{k \neq 0} c_k \frac{1}{N} \frac{e^{2i\pi kA} - 1}{e^{2i\pi kA} - 1}, \qquad (10.4)$$

et à étudier d'une façon précise le comportement du dénominateur. On a

$$|e^{2i\pi kA} - 1| = 2|\sin \pi kA|.$$

On remarque ensuite que :

$$\text{si } 0 < x < \frac{1}{2}, \qquad |\sin \pi x| > 2x,$$

$$\text{si } \frac{1}{2} < x < 1, \qquad |\sin \pi x| > 2(1 - x),$$

et que, d'une façon générale,

$$|\sin \pi x| > 2 |x - 1|, \qquad (10.5)$$

où l est l'entier le plus voisin de x.

On a donc

$$|e^{2i\pi kA} - 1| > 4|kA - 1|.$$

D'après le théorème de Roth, on a, quel que soit l, l'inégalité
$$|kA - 1| > \frac{1}{k^{\delta+1}}, \qquad k > 0,$$

(sauf peut-être pour un nombre fini de valeurs de k).

Il en résulte que, à partir d'une certaine valeur de $|k|,$

$$ \left| \frac{e^{2i\pi kAN} - 1}{e^{2i\pi kA} - 1} \right| < \frac{|k|^{\delta+1}}{2}, \qquad (10.6) $$

et que la série dont $I_N$ est la somme est majorée par la série ayant pour terme général

$$ \frac{1}{N} \frac{|c_k|}{2} |k|^{\delta+1}. \qquad (10.7) $$

Supposons pour fixer les idées que l'inégalité de Roth soit valable sans restrictions. Dans le cas contraire, il faudrait isoler et calculer à part les termes, en nombre fini, qui ne sont pas justiciables de l'inégalité. On a alors

$$ |I_N| < \frac{1}{N} \sum_{k=1}^{\infty} |c_k| k^{\delta+1}, \qquad (10.8) $$

la somme étant étendue aux valeurs positives de k.

Cette inégalité fournit une majoration de l'erreur $I_N$ lorsque la série du second membre est convergente. Or par hypothèse,

$$ |c_k| < \frac{M}{|k|^{\lambda+1}}, $$

d'où

$$ |k|^{\delta+1} |c_k| < \frac{M}{k^{\lambda-\delta}}. \qquad (10.9) $$

Il y a convergence si $\lambda > 1 + \delta$. Comm\o $\delta$ peut être choisi arbitrairement petit, il suffit que $\lambda > 1$. Le résultat de Richtmyer s'énonce ainsi de la façon suivante :

## THEOREME 10 -

Si $f(x)$ admet sur l'intervalle $(0,1)$ un développement en série de Fourier dont les coefficients sont majorés par ceux de la série $\frac{M}{k^{1+\lambda}}$ où
$\lambda$ est un nombre supérieur à 1, l'erreur commise sur le calcul de l'in-
tégrale $\int_{0}^{1} f(x)dx$ en employant la suite arithmétique An (mod 1), où A
est un nombre algébrique, tend vers 0 comme $\frac{1}{N}$. Ce résultat est en
particulier valable si $f(x)$, considérée comme une fonction périodique
de période 1, admet une dérivée seconde continue, et une dérivée troi-
sième qui soit continue sauf peut-être en un nombre fini de sauts dans
chaque période.

Lorsque A est un nombre quadratique, on peut remplacer $\delta$ par 0, et la démonstration de ce théorème ne fait pas appel au théorème de Roth, mais seulement au théorème élémentaire de Liouville.

XI - QUELQUES RESULTATS SUR LE CALCUL ARITHMETIQUE DES INTEGRALES A L'AIDE DES POLYNOMES DE DEGRE SUPERIEUR AU PREMIER -

Les méthodes utilisées pour étudier la suite An tombent en défaut dès qu'on remplace An par un polynôme de degré au moins égal à 2. On connaît assez peu de résultats précis sur le comportement de la différence

$$
\int_{0}^{1} f(x) dx - \frac{1}{N} \sum_{n=0}^{N-1} f(x_n)
$$

lorsque $x_n$ est de la forme $\varphi(n)$ mod 1, où $\varphi$ est un polynôme ayant un
coefficient irrationnel. Cette différence tend vers 0, et nous n'étu-
dierons pas plus en détail comment elle tend vers 0.

Nous nous contenterons d'examiner ce qui se passe pour

$\varphi(n) = An^2$

et lorsque A est remplacé par une fraction irréductible $\frac{p}{q}$.

Il est commode de supposer que $f(x)$ est développable en série de Fourier. Tout revient donc à l'étude du cas particulier où

$f(x) = e^{2 
ament}$

et nous nous limiterons même à l'étude des sommes

$$
\frac{1}{q} \sum_{n=0}^{q-1} e^{2i\pi k \frac{p}{q} n^2} \quad (11.1)
$$

pour $N = q$.
Si k est multiple de q, ces sommes sont égales à 1, comme les sommes $\frac{1}{N} \sum_{n=0}^{N-1} e^{\frac{2i\pi k}{q}n}$ précédemment étudiées.

Si k n'est pas multiple de q, ce n'est pas une restriction de choisir k = 1. Cela revient à faire jouer à la lettre p le rôle qui était primitivement dévolu à kp. Considérons donc la moyenne

$$S = \frac{1}{q} \sum_{n=0}^{q-1} e^{2i\pi \frac{p}{q} n^2} \quad (11.2)$$

Introduisons la fonction g(t), périodique, de période q, qui est égale à

$$e^{\frac{2i\pi}{q} n^2} \quad \text{sin } n < t < n + 1.$$

C'est une fonction en escalier, développable en une série de Fourier. g(t) est de carré sommable, et la formule de Parseval lui est applicable. S est égale à la valeur moyenne de g(t) sur une période. Si donc on écrit

on a

$$g(t) = \sum_{j=-\infty}^{\infty} a_j e^{\frac{2i\pi}{q} j t}, \quad (11.3)$$

$$S = a_0.$$

On sait d'ailleurs que, si

$$\gamma(h) = \frac{1}{q} \int_{0}^{q} g(t) g(t+h) dt \quad (11.4)$$

désigne la "fonction de corrélation" de g(t), on peut représenter $\gamma(h)$ par une série de Fourier ayant comme coefficients les carrés des modèles de ceux du développement de g(t) :

$$\gamma(h) = \sum_{j=-\infty}^{\infty} |a_j|^2 e^{\frac{2i\pi}{q} jh} \quad (11.5)$$

Il en résulte que $|a_0|^2$ est égal à la moyenne de $\gamma(h)$ sur une période : le carré du module de la somme S est égal à la moyenne de la fonction $\gamma(h)$ sur l'intervalle $0 < h < q$.

On démontre d'autre part [1] que la fonction $\gamma(h)$ est continue et que, dans chaque intervalle $(n, n+1)$, elle varie linéairement. Ses
valeurs pour h entier suffisent donc à la déterminer. Nous allons les calculer.

Si h est *entier*, on a

$$ \gamma(h) = \frac{1}{q} \sum_{n=0}^{q-1} e^{\frac{2i\pi}{q} [(n+h)^2 - n^2]} = e^{\frac{2i\pi h^2}{q}} \frac{1}{q} \sum_{n=0}^{q-1} e^{\frac{4i\pi}{q} nh} \quad (11.6) $$

L'expression de $\gamma(h)$ est plus simple que celle de S. Il n'y figure plus qu'une progression géométrique. On trouve naturellement que $\gamma(h) = 1$ pour $h = 0$ et $h = q$. $\gamma(h)$ est nul pour presque toutes les autres valeurs de h. Il ne peut y avoir exception que si $\frac{2h}{q}$ est entier.

Si q est impair, cette circonstance n'a jamais lieu, et $\gamma(h) = 0$ pour tout $h \neq 0$ mod q.

Si $q = 2q'$, on a

$$ \gamma(q') = e^{i\pi q'} \frac{1}{2q'} \sum_{n=0}^{2q'-1} e^{2i\pi n} = e^{i\pi q'}. \quad (11.7) $$

$\gamma(q') = 1$ ou $-1$ suivant que $q'$ est pair ou impair.

Dans le premier cas ($q$ impair), $\gamma(h)$ a pour moyenne $\frac{1}{q}$.

Dans le second cas ($q = 0 \mod 4$), $\gamma(h)$ a pour moyenne $\frac{2}{q}$.

Dans le troisième cas ($q = 2 \mod 4$), $\gamma(h)$ a pour moyenne 0.

Il en résulte que, suivant les cas, $|S|$ a pour valeur

$$ 0, \quad \frac{1}{\sqrt{q}}, \quad \frac{\sqrt{2}}{\sqrt{q}}. $$

## THEOREME 11 -

Le module de la somme $S = \frac{1}{q} \sum_{n=0}^{q-1} e^{2i\pi \frac{p}{q} n^2}$ est égal à

$$ \frac{1}{\sqrt{q}} \quad \text{si q est impair} $$

$$ \frac{\sqrt{2}}{\sqrt{q}} \quad \text{si q = 0 mod 4} $$

$$ 0 \quad \text{si q = 2 mod 4.} $$
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Figure 12</figcaption></figure>
On remarque que la somme

$$ \frac{1}{q} \sum_{n=0}^{q-1} e^{2i\pi \frac{p}{q}n} $$

étudiée dans les paragraphes précédents, est toujours nulle. Ceci montre que, lorsqu'on se limite aux sommes $\frac{1}{q} \sum_{n=0}^{q-1}$, l'approximation d'une intégrale à l'aide de suites $\frac{p}{q} n^2$ est moins bonne qu'à l'aide des suites $\frac{p}{q} n$. Dans le premier cas, la différence entre l'intégrale et sa valeur approchée contient des termes qui sont petits avec $\frac{1}{q}$, mais ne sont pas nuls, alors que ces termes disparaissent identiquement dans le second.

Remarque.

Les sommes $\sum_{n=0}^{q-1} e^{2i\pi \frac{p}{q}n^2}$ contiennent comme cas particulier, pour $p = 1$, les sommes de Gauss

$$ \frac{1}{q} \sum_{n=0}^{q-1} e^{2i\pi \frac{p}{q}n^2} $$

En intégrant la fonction $\frac{e^{2i\pi \frac{p}{q}z^2}}{e^{2i\pi z} - 1}$ sur un contour convenable du plan complexe, on démontre que les sommes de Gauss valent

$$ \frac{1+i}{\sqrt{q}}, \quad \frac{1}{\sqrt{q}}, \quad 0, \quad \frac{i}{\sqrt{q}} $$

suivant que $q$ prend, mod 4, les valeurs

$$ 0, \quad 1, \quad 2, \quad 3 $$

La méthode employée dans ce paragraphe ne fournit que les modèles des sommes de Gauss, qui suffisent d'ailleurs pour discuter le problème posé.

XII - CONCLUSIONS -

De l'étude qui précède, on ne peut pas tirer de conclusions ab solues et définitives. Sous réserve de vérifications plus poussées, il semble qu'on puisse adopter les règles de conduite que voici :

Pour calculer une intégrale par la méthode de Monte-Carlo, il y a intérêt à employer des suites arithmétiques plutôt que des suites de véritables nombres aléatoires.
Comme suites arithmétiques uniformément denses, on connaît celles qui résultent des théorèmes de H. Weyl et quelques suites de structure voisine. Les suites de polynômes $\varphi(n)$ mod 1, que nous avons étudiées, ne sont pas les seules qu'on sache construire, mais ce sont très probablement les plus maniables. Elles se prêtent bien au calcul électronique.

Parmi les polynômes $\varphi(n)$, ceux qui donnent la meilleure approximation sont les polynômes du premier degré $A_n$, où $A$ est un nombre irrationnel. L'approximation est d'autant meilleure que la fonction à intégrer est plus régulière. On sait évaluer, assez grossièrement, une limite supérieure de l'erreur commise en remplaçant

$$ \int_{0}^{1} f(x) dx \quad \text{par} \quad \frac{1}{N} \sum_{n=0}^{N-1} f(x_n), $$

où $x_n = A_n \mod 1$, dès que $f(x)$ admet un développement en série de Fourier dont les coefficients $c_k$ sont majorés par $\frac{M}{k^{1+\lambda}}$, $\lambda > 0$.

Il semble qu'il y ait intérêt à choisir pour $A$ non pas un nombre transcendant, mais un nombre algébrique, et même un nombre quadratique. (1)

Si $A$ est quadratique, et si $\lambda > 1$, on sait que l'erreur commise est inférieure à $\frac{H}{N}$, où $H$ est un nombre que l'on sait calculer.

La méthode arithmétique s'applique au calcul des intégrales multiples. On peut par exemple évaluer approximativement l'intégrale double

$$ \iint_C f(x, y) dx dy $$

étendue au carré $0 < x < 1$, $0 < y < 1$, en la remplaçant par la somme

$$ \frac{1}{N} \sum_{n=0}^{N-1} f(x_n, y_n), $$

où la suite $x_n, y_n$ est équirépartie dans le carré $C$. On sait qu'il en est ainsi si la suite

$$ l_1 x_n + l_2 y_n $$

<footnote>(1) Par une étude plus détaillée de la structure du nombre irrationnel $A$, on peut sensiblement améliorer ces conclusions. Voir à ce sujet un prochain article de M. Bertrandias, où seront en outre développés plusieurs exemples numériques.</footnote>
est équirépartie mod 1 pour tous les couples d'entiers $l_1, l_2$ non simul-
tanément nuls.

Dans le cas de suites $x_n = A_n, y_n = B_n$, cette condition est réa-
lisée si les nombres irrationnels A et B ne sont pas liés par une rela-
tion à coefficients entiers de la forme

$$
l_1 A + l_2 B = l_3
$$

Exemple:

$$
A = \sqrt{2},
$$

$$
B = \sqrt{3}.
$$

La suite $x_n = A_n, y_n = B_n(n+1) \pmod 1$ ne satisfait pas à cette
condition.

Si φ(n) est un polynôme de degré 2 au moins, ayant un coefficient
irrationnel (autre que le terme constant), la suite

$x_n = \phi(n),$

$y_n = \phi(n + 1),$

mod 1

est équirépartie dans le carré C. Elle peut servir au calcul des inté-
grales doubles.

La suite An semble la plus favorable au calcul des intégrales.
Mais elle présente une corrélation interne. Si l'on a besoin de résou-
dre par la méthode de Monte-Carlo des problèmes dans lesquels il est
indispensable que la suite de nombres aléatoires ne présente pas de
corrélation interne, il faut faire appel à des polynômes φ(n) de degré
2 au moins, par exemple An² mod 1. La précision des résultats est
alors difficile à prévoir. On ne sait pas calculer d'avance le nombre N
de termes qu'il faut calculer pour que l'erreur commise soit inférieure
à une valeur imposée.

La suite des parties entières mod 2 de φ(n) est, du point de vue
arithmétique, équivalente à une suite de nombres aléatoires prenant
les valeurs 0 et 1, sans corrélation interne.
BIBLIOGRAPHIE

[1] J.BASS - Suites uniformément denses, moyennes trigonométriques, fonctions pseudo-aléatoires. Bulletin de la Soc. Math. de France, 87, 1959, p. 1 à 64.

[2] J.BASS et J.GUILLOUD - Méthode de Monte-Carlo et suites uniformément denses - Chiffres, 1, 1958, p. 149 à 155.

[3] J.W.S. CASSELS - An introduction to diophantine approximation, Cambridge University Press, 1957.

[4] R.D. RICHTMYER - The evaluation of definite integrals, and a quasi Monte-Carlo method based on the properties of algebraic numbers.

Los Alamos scientific laboratory - University of California, 1952.

[5] H. WEYL - Ueber die Gleichverteilung von Zahlen modulo Eins
Math Annalen, 77, 1916, p. 313-352.