# Une méthode pour choisir les paramètres d'un générateur de nombres aléatoires congruentiel

Yoshiharu Kurita

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

## To cite this version:

Yoshiharu Kurita. Une méthode pour choisir les paramètres d'un générateur de nombres aléatoires congruentiel. IRIA-RR-199, IRIA. 1976, pp.13. hal-04831699

HAL Id: hal-04831699

https://inria.hal.science/hal-04831699v1

Submitted on 11 Dec 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
IRIA
laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B.P. 105 78 150 - Le Chesnay
France
Tél. 954 90 20

laboratoire de recherche
en informatique
et automatique

UNE MÉTHODE POUR
CHOISIR LES PARAMÈTRES
D'UN GÉNÉRATEUR
DE NOMBRES ALÉATOIRES
CONGRUENTIEL

Yoshiharu KURITA

Rapport de Recherche N° 199
Novembre 1976
# UNE MÉTHODE POUR CHOISIR LES PARAMÈTRES D'UN GÉNÉRATEUR DE NOMBRES ALÉATOIRES CONGRUENTIEL

Yoshiharu Kurita

## Résumé :

Nous étudions dans ce papier un générateur de nombres aléatoires où la génération des nombres est obtenue par une méthode congruentielle. Nous indiquons comment il faut choisir les paramètres du générateur pour lui assurer une période maximale d'une part et de bonnes qualités statistiques d'autre part.

## Abstract

In this paper we are studying a random number generator where the sequence of the random numbers is obtained from a congruential method. We point out the influence of the different parameters and we indicate how we have to choose them to ensure a maximal period and good statistical properties.
UNE METHODE POUR CHOISIR LES PARAMETRES
UN GENERATEUR DE NOMBRES ALEATOIRES CONGRUENTIEL

Yoshiharu KURITA

National de Recherche en Métrologie (Tokyo)*
et

- LABORIA**

INTRODUCTION

Une des caractéristiques des ordinateurs contemporains est leur comportement déterministe. Aussi lorsqu'on génère par ordinateur une suite de nombres dans le but de représenter les réalisations d'une variable aléatoire, parle-t-on de nombres pseudo-aléatoires pour indiquer le déterminisme des algorithmes qui ont servi à la génération des nombres. Par ailleurs, s'il est possible de définir le caractère aléatoire d'une suite infinie, il n'existe pas de définition généralement acceptée pour une suite finie, ce qui est inévitablement le cas d'une suite générée par ordinateur. Nous voyons donc le paradoxe qui recouvre le mot pseudo. Nous dirons qu'une suite de nombres est pseudo-aléatoire, si "c'est une suite finie dans laquelle il est difficile de constater des régularités évidentes."

De nombreuses méthodes ont été proposées pour générer des nombres pseudo-aléatoires. Les algorithmes les plus fréquemment utilisés ont été jusqu'à maintenant ceux qui sont basés sur les méthodes congruentielles. Pourtant une technique récente fondée sur les registres à décalage bouclés génèrent des suites de TAUSWORTHE [6] tend à se répandre. Le défaut majeur de cette méthode est de ne avoir été suffisamment testée statistiquement du fait de sa nouveauté, ce qui n'est pas le cas des générateurs congruentiels.

Dans ce papier nous présentons la méthodologie qu'il importe de suivre pour choisir les paramètres d'un générateur congruentiel et nous justifions les choix que nous proposons par des résultats numériques.

<footnote>* 4 - 10, 1= chôme, Kaga Itabashi - TOKYO, 173 - JAPON</footnote>

<footnote>** B. P. 5 - 78150 LE CHESNAY - FRANCE</footnote>
## LA METHODE_CONGRUENTIELLE_MULTIPLICATIVE.

Parmi les différentes méthodes congruentielles de génération de nombres aléatoires, nous nous bornerons à l'étude du générateur congruentiel multiplicatif proposé. 1951 par LEHAER. Le deuxième nombre aléatoire fourni par ce générateur est donné par la récurrence:

$$x_{n+1} = x_n \mod m$$

Dans la formule (1) nous remarquons qu'il nous faut choisir $x_0$, $m$ et $a$. Le choix de $x_0$ n'est pas essentiel, en pratique la valeur de $x_0$ sera déterminée par des tests statistiques effectués dans les conditions des expériences par l'utilisateur. Par contre on essaiera d'avoir une période pour la suite de nombres générés aussi grande que possible. Cette contrainte impose à m d'avoir certaines propriétés.

### Recherche d'une période maximale

D'après (1) il est bien évident que la période de la suite des nombres générés ne peut pas dépasser $m$, on choisira donc $m$ aussi grand que possible. En général $m$ sera la plus grand nombre représentable dans le registre de l'ordinateur utilisé, mais ce n'est pas pour autant que la période $T$ du générateur sera grande.

De (1) on tire aisément : $x_n = a^n x_0 \mod m$, et la période $T$ est alors :

$$T = \min(\lambda \mid a^{\lambda} = 1 \mod m).$$

Par définition on dira que $\lambda$ est racine primitive de $m$ si et seulement $T = m-1$.

Ainsi, si l'on choisit $m$ tel que ce nombre admette des racines primitives et parmi ces racines primitives, on sera assuré d'avoir une période maximale ($T = m - 1$) et grande si $m$ est grand.

Pour satisfaire les conditions d'existence de racines primitives, $m$ devra avoir les formes suivantes :

$$ (2) \qquad m = p^s \text{ ou } 2p^s $$

où $p$ est un nombre premier différent de 2 et $s = 1,2,...$. Pour pouvoir effectuer sur un ordinateur, l'opération de congruence modulo $m$ sans division [4], on choisit :

$$ (3) \qquad m = \dots $$
## 2-2 Coefficient d'autocorrélation d'ordre 1.

Il semble normal d'imposer que le coefficient d'autocorrélation d'ordre 1 sur une période soit petit, ce qui indique qu'il y aura peu de liaison, dans la suite générée, entre un nombre et son successeur. L'avantage de ce critère est de pouvoir être calculé assez aisément quand la période du générateur est maximale.

Le coefficient d'autocorrélation d'ordre 1 est défini:

$$= \frac{\sum_{n=1}^{m-1} x_n (ax_n \mod m) - \left(\sum_{n=1}^{m-1} x_n\right)^2}{\sum_{n=1}^{m-1} x_n^2 - \left(\sum_{n=1}^{m-1} x_n\right)^2}$$

Si la suite {$x_n$} des nombres générés a une période maximale (m-1), on peut alors calculer C par une somme de Dedekind généralisée [2].

Nous donnons sur la figure 1 les valeurs de C obtenues pour $m = 2^{31}$, $a = q, q+1, q+2, ..., q+50$ avec $q = 2^9, 2^{10}, \sqrt{m}, ..., m-\sqrt{m}$.

<figure><img src="image_3.png" /><figcaption>Fig. 1 : coefficient d'autocorrélation C calculé sur la période complète.</figcaption></figure>
un voit sur la figure que si $a/m \ll 1$, alors $C \approx 1/a$.

Dans le cas contraire, lorsque $a$ est voisin de $m$, $C$ peut devenir très grand, mais même dans ce cas, il est clair que l'on peut trouver un $a$ pour lequel $C$ est petit.

Il faut noter que le coefficient $C$ est calculé sur toute la période du générateur et n'apporte donc aucune garantie sur les qualités statistiques d'une sous-suite qui en est extraite. Or, dans la pratique, on n'utilise habituellement que des sous-suites. Ce critère ne semble donc pas suffisant pour s'assurer de la qualité d'un générateur de nombres aléatoires.

## 2-3 Analyse spectrale des nombres générés

Nous proposons maintenant un test plus sévère connu sous les termes de "structure d'hyperplan" ou de "régularité cristalline". Le but de ces tests est d'étudier les répétitions dans la suite des nombres générés des sous-suites de longueur $r$. $r = 2, \dots, m-1$. On peut envisager deux approches, l'une géométrique [3], l'autre analytique par une méthode d'analyse spectrale [1].

### 2-3.1 Représentation géométrique

Soit $\{x_n\}$ la suite définie par

$$ \underbrace{x}_{n} \equiv \frac{x_{n-1}}{\mod m} $$

où $a$ et $m$ sont choisis de telle sorte que tous les points de $\{x_n\}$ appartient-à un seul cycle, et $x_i$, $i = 1, 2, \dots$ les sous-suites de $\{x_n\}$ de longueur $r$ définies par

$$ x_n \equiv x, \quad r $$
L'approche géométrique consiste à considérer les $n_i$ comme des points de coordonnées $(x_i, x_{i+1}, ..., x_{i+r-1})$ dans un espace à r dimensions. En raison de leur construction, ces points forment un réseau d'hyperplans parallèles et équidistants, comme représenté sur la Figure 2 pour r = 2, et on en déduit que le générateur sera d'autant plus régulier que le nombre minimal $v_r$ d'hyperplans obtenu sera plus grand. On montre que ce nombre est donné par

$$v_r = \min (s_1^2 + s_2^2 + \dots + s_r^2)^{1/2}$$

où les $s_i$, i = 1, 2, ..., r, satisfont à la congruence

$$s_1 + a s_2 + \dots + a^{r-1} s_r = 0 \mod m.$$

<figure><img src="image_4.png" /><figcaption>Figure 2 : structures d'hyperplan pour r=2</figcaption></figure>
### 2.3.2 Représentation par analyse spectrale

Cette approche a été développée en [1]. Considérons une suite $\{x_n\}$ telle que tous ses points appartiennent à un seul cycle, et soit $F(t_1, t_2, ..., t_r)$ la densité limite du nombre d'apparitions des sous-suites $\{t_1, t_2, ..., t_r\}$ parmi toutes les sous-suites d'ordre r. Nous pouvons écrire

$$ \lim_{N \to \infty} \sum_{k=0,1,...,N-1} \delta_{x_k, t_1} \delta_{x_{k+1}, t_2} \dots \delta_{x_{k+r-1}, t_r} $$

$$ \delta_{i,j} = \begin{cases} 1, & \text{si } j < i \\ 0, & \text{sinon} \end{cases} $$

suite $\{x_n\}$ est uniforme, nous devons avoir

$$ F(t_1, t_2, ..., t_r) = 1/m^r, $$

et dans ce cas nous dirons que la suite est uniforme d'ordre r. Afin d'analyser cette propriété calculons la transformée de Fourier de $F(t_1, t_2, ..., t_r)$.
La transformée de Fourier de $s$ écrit

$$ (4) \quad f(s_1, s_2, ..., s_r) = \delta[\Phi(a)/m] $$

$$ \text{où } \delta(z) = 1 \text{ si } z \in \mathbb{N} \quad , \quad 0 \text{ sinon,} \\ \text{et } \Phi(a) = s_1 + as_2 + \dots + a^{r-1}s_r. $$

Si $\{x_n\}$ est uniforme d'ordre r, alors

$$ f(s_1, s_2, ..., s_r) = 1 \text{ si } s_1 \equiv s_1 \equiv s_2 \equiv \dots \equiv s_r \equiv 0 \pmod{m}, \\ 0 \text{ sinon.} $$

Une mesure de l'irrégularité de $\{x_n\}$ sera donc obtenue en examinant les vecteurs $(s_1, s_2, ..., s_r) \neq 0$ pour lesquels $f(s_1, s_2, ..., s_r) \neq 0$. Chaque nombre tel $f/m^r$ est l'amplitude de l'harmonique de longueur d'onde $L_r$ dans la décomposition de Fourier de $F(t_1, t_2, ..., t_r)$, où

$$ L_r = 1/(s_1^2 + s_2^2 + \dots + s_r^2)^{\frac{1}{2}} $$

Cet harmonique est source de perturbations dans la suite si sa longueur d'onde est élevée, et donc on cherchera à minimiser le nombre d'onde $v_r = 1/L_r$. Pour une suite multiplicative $x_n = a x_{n-1} \pmod{m}$, $f \neq 0$ implique, d'après (4), $\varphi(a) = 0 \pmod{m}$, et on est donc ramené à calculer

$$ \varphi(a) = \min \left[ s_1^2 + s_2^2 + \dots + s_r^2 \right]^{\frac{1}{2}} $$

avec

$$ a s_2 + \dots + a^{r-1} s_r = 0 \pmod{m} $$
Sous ces deux approches, ce test se fera donc en calculant les nombres $v_r$ suivant un algorithme présenté en [2], et le choix de $a$ doit être tel que $v_r$ soit le plus grand possible. Cependant, $v_r$ ne peut être arbitrairement grand, en raison du Théorème d'Hermite sur le minimum d'une forme quadratique définie positive qui entraîne :

$r \quad \gamma_r \stackrel{1/r}{\sim} \quad \text{où } \gamma_r \text{ est une constante [7].}$

La conséquence importante de cette limitation est qu'il est impossible d'avoir un générateur de nombres pseudo-aléatoires congruentiel arbitrairement bon par un choix convenable des paramètres.

Nous présentons dans le paragraphe suivant des exemples numériques de calcul de $v_r$ tel que $v_r$ soit maximisé.

### III. RESULTATS NUMÉRIQUES

A partir des résultats sur le coefficient d'autocorrélation présentés dans le paragraphe 2.2, nous pouvons choisir, pour $m = 2^{31} - 1$, $a$ dans l'intervalle suivant :

$$2 \ 100 \ 000 \ 000 \le a \le 2 \ 100 \ 000 \ 000 + 30 \ 206$$

Le nombre de racines primitives situées dans cet intervalle est 7 440, et chacune a été examinée par "passage au crible" de la façon suivante :

1) prendre la première racine primitive

2) **pour** $r = 2$ **par pas de** $1$ **jusqu'à** $7$
   **début** si $v_r < \bar{v}_r$ **alors** **aller à** 4)
   **sinon**
   **fin**

3) calculer $v_8$ et adopter cette racine primitive

4) prendre la racine primitive suivante, **aller à** 2)

Les $\bar{v}_r$, $r=1,..,7$ représentent la finesse du crible. Ils ont été choisis entre 60 % et 75 % de la borne supérieure $\gamma_r \stackrel{1/r}{\sim} 30 000$.

$\bar{v}_3 = 1 086$, $\bar{v}_4 = 192$, $\bar{v}_5 = 60$, $\bar{v}_6 = 30$, $\bar{v}_7 = 19$, 11 racines primitives ont passé le crible sur les 7 440 possibles. Elles sont rapportées dans la table 1 avec les valeurs de $v_r$, $r=1,..,7$ et du coefficient d'autocorrélation $C$.
<table><thead><tr><th>a</th><th>v<sub>2</sub></th><th>v<sub>3</sub></th><th>v<sub>4</sub></th><th>v<sub>5</sub></th><th>v<sub>6</sub></th><th>v<sub>7</sub></th><th>C</th></tr></thead><tbody><tr><td>2100004246</td><td>35497.3</td><td>1919.9</td><td>212.3</td><td>64.0</td><td>30.9</td><td>19.5</td><td>-2.8</td></tr><tr><td>21000055341</td><td>43487.0</td><td>1201.5</td><td>209.7</td><td>65.2</td><td>31.9</td><td>19.3</td><td>1.1</td></tr><tr><td>2100008090</td><td>41396.6</td><td>1223.1</td><td>197.8</td><td>60.8</td><td>31.7</td><td>19.6</td><td>-1.3</td></tr><tr><td>2100009038</td><td>39752.6</td><td>1100.5</td><td>194.1</td><td>68.9</td><td>34.3</td><td>19.5</td><td>-5.1</td></tr><tr><td>2100009048</td><td>30328.3</td><td>1155.3</td><td>206.9</td><td>65.0</td><td>31.0</td><td>19.5</td><td>-0.5</td></tr><tr><td>2100009052</td><td>34471.2</td><td>1334.2</td><td>195.3</td><td>69.0</td><td>33.1</td><td>19.9</td><td>-2.2</td></tr><tr><td>2100014651</td><td>31168.0</td><td>1279.3</td><td>209.2</td><td>64.0</td><td>35.2</td><td>20.3</td><td>-4.7</td></tr><tr><td>2100017018</td><td>35193.2</td><td>1154.9</td><td>223.4</td><td>70.3</td><td>34.0</td><td>19.3</td><td>8.3</td></tr><tr><td>2100017008</td><td>44643.0</td><td>1196.5</td><td>194.2</td><td>68.5</td><td>31.8</td><td>21.2</td><td>-1.4</td></tr><tr><td>2100019187</td><td>41742.2</td><td>1177.4</td><td>203.6</td><td>74.2</td><td>33.2</td><td>19.2</td><td>-0.8</td></tr><tr><td>2100023033</td><td>31118.3</td><td>1164.5</td><td>209.7</td><td>65.2</td><td>33.3</td><td>20.5</td><td>-3.8</td></tr></tbody><tfoot><tr><td>borne sup.</td><td>(44746.6)</td><td>(1448.2)</td><td>(256.0)</td><td>(40.5)</td><td>(46.4)</td><td>(28.0)</td><td></td></tr></tfoot></table>

Table 1 : racines primitives candidates au rôle de a
pour la suite $x_n = a \cdot x_{n-1} \pmod{2^{31} - 1}$ maximisant
, r=1,...,7.

La Table 2 présente des résultats semblables avec la condition supplémentaire
que a soit de la forme

$$a = 2^\alpha + 2^\beta + 1$$

afin d'effectuer rapidement la multiplication $a \cdot x_n$.

<table><thead><tr><th>m</th><th>a</th><th>v<sub>2</sub></th><th>v<sub>3</sub></th><th>v<sub>4</sub></th><th>v<sub>5</sub></th><th>C</th></tr></thead><tbody><tr><td rowspan="3">m=2<sup>47-2<sup>7</sup>+1</sup></td><td>2<sup>36-2<sup>12</sup></sup></td><td>8.4 x 10<sup>6</sup></td><td>4.7 x 10<sup>4</sup></td><td>2.3 x 10<sup>3</sup></td><td>5.6 x 10<sup>2</sup></td><td>-4.4 x 10<sup>-11</sup></td></tr><tr><td>2<sup>44-2<sup>28</sup>+1</sup></td><td>8.9 x 10<sup>6</sup></td><td>4.9 x 10<sup>4</sup></td><td>2.1 x 10<sup>3</sup></td><td>5.5 x 10<sup>2</sup></td><td>-5.8 x 10<sup>-11</sup></td></tr><tr><td>borne sup</td><td>1.2 x 10<sup>7</sup></td><td>5.8 x 10<sup>4</sup></td><td>4.1 x 10<sup>3</sup></td><td>8.3 x 10<sup>2</sup></td><td></td></tr><tr><td rowspan="3">m=2<sup>35-2<sup>5</sup>+1</sup></td><td>2<sup>27-2<sup>18</sup>-1</sup></td><td>1.7 x 10<sup>5</sup></td><td>2.1 x 10<sup>3</sup></td><td>3.1 x 10<sup>2</sup></td><td>1.1 x 10<sup>2</sup></td><td>2.6 x 10<sup>-10</sup></td></tr><tr><td>2<sup>30-2<sup>18</sup>-1</sup></td><td>1.3 x 10<sup>5</sup></td><td>2.9 x 10<sup>3</sup></td><td>3.4 x 10<sup>2</sup></td><td>1.1 x 10<sup>2</sup></td><td>-5.4 x 10<sup>-9</sup></td></tr><tr><td>borne sup</td><td>1.9 x 10<sup>9</sup></td><td>3.6 x 10<sup>3</sup></td><td>5.1 x 10<sup>2</sup></td><td>1.5 x 10<sup>2</sup></td><td></td></tr><tr><td rowspan="3">m=2<sup>31-1</sup></td><td>2<sup>22-2<sup>14</sup>+4</sup></td><td>3.1 x 10<sup>4</sup></td><td>1.2 x 10<sup>3</sup></td><td>1.8 x 10<sup>2</sup></td><td>6.6 x 10<sup>2</sup></td><td>1.3 x 10<sup>-7</sup></td></tr><tr><td>2<sup>26-2<sup>16</sup>+2</sup></td><td>3.4 x 10<sup>4</sup></td><td>1.2 x 10<sup>3</sup></td><td>1.8 x 10<sup>2</sup></td><td>5.1 x 10<sup>2</sup></td><td>-5.4 x 10<sup>-9</sup></td></tr><tr><td>borne sup.</td><td>4.9 x 10<sup>4</sup></td><td>1.4 x 10<sup>3</sup></td><td>2.6 x 10<sup>2</sup></td><td>9.0 x 10<sup>2</sup></td><td></td></tr></tbody></table>

Table 2 : racines primitives candidates au rôle de a avant une forme spéciale
pour $x_n = a \cdot x_{n-1} \pmod m$
-9-
IV. CONCLUSION

Nous avons présenté dans cet article une méthode pour le choix du paramètre multiplicatif d'un générateur de nombres aléatoires congruentiels. La méthode est fondée sur deux tests : le premier test permet de déterminer un intervalle de valeurs du paramètre qui assure une faible autocorrélation d'ordre 1 ; les valeurs acceptables du paramètre dans cet intervalle sont ensuite passées dans un crible, et seules sont retenues celles qui engendrent une répartition aussi uniforme que possible des sous-suites d'ordre 2,3. Les tests sont illustrés par des exemples numériques pour différentes valeurs de paramètres des générateurs.

Il resterait à vérifier que les qualités des générateurs ainsi construits se traduisent par une amélioration des résultats obtenus lorsqu'ils sont utilisés par exemple dans une expérience de simulation. Une étude portant sur la simulation d'une file d'attente M/1 M/1 est en cours au LABORIA, et les premiers résultats indiquent une amélioration sensible. D'autres expériences seraient nécessaires et il s'agit là d'un champ de recherche encore largement inexploré.
Références

COVEYOU, R.R. ; Mac PHERSON, R.D., "Fourier analysis of Uniform
random numbers" J. A. C. M. Vol. 14 (1967)

KNUTH,D.E., "The art of Computer Programming", Vol. 2 - Seminu-
merical Algorithms - Adison Wesley (1969)

[3] MARSAGLIA, G., "Random number fall mainly on the planes"
Proc. N. A. S. Vol. 61 (1968)

[4] PAYNE, W.H., "Coding the Lehmer pseudo-random number generator"
C.A.C.M. Vol. 12 (1969)

[5] SIERPINSKI, W., "Elementary theory of numbers" Porska Academia
Nauk (1964)

[6] TAUSWORTHE, R.C. ; "Random number generated by linear recurrence
modulo two" Math. Comp. Vol. 19 (1965)

[7] CASSELS, J.W.S., "An introduction to the Geometry of Numbers"-
Springer, Berlin (1959)