Nombres aléatoires et biais dû au modulo

La génération de nombres aléatoires est un problème étonnamment complexe. Heureusement, sauf besoin spécifique, le développeur y est rarement confronté directement : la bibliothèque C fournit des générateurs pseudo-aléatoires performants et prêts à l'emploi. Les systèmes basés sur FreeBSD proposent notamment random() et arc4random(), qui sont de très bons remplacements à l'historique et défectueuse fonction rand().

En revanche, un problème auquel le développeur est souvent confronté est de générer des nombres pseudo-aléatoires compris entre 0 et N à partir d'un générateur qui produit des nombres compris entre 0 et RAND_MAX. Une technique souvent rencontrée consiste à appliquer un simple modulo N. Par exemple, pour générer un nombre pseudo-aléatoire entre 0 et 4 inclus :

int x = random() % 5;

Le problème de cet algorithme naïf (et extrêmement répandu…) est qu'il introduit un biais : les nombres 0, 1 et 2 ont plus de chance de sortir que les nombres 3 et 4. Cette perte d'uniformité s'explique par le fait que RAND_MAX + 1, qui vaut typiquement une puissance de 2 sur la plupart des plateformes, n'est pas un multiple de 5. Plus généralement, à partir du moment où (RAND_MAX + 1) % N != 0, alors il existe un biais.

Pour comprendre ce qui se passe, imaginons que la fonction random() retourne un nombre compris entre 0 et 7 inclus (ce n'est pas une hypothèse absurde, les générateurs pseudo-aléatoires retournent souvent par construction des nombres compris en 0 et 2ⁿ-1, or 7 vaut justement 2³-1). Quels nombres donneront ces valeurs après un modulo 5 ?

0 donnera 0.
1 donnera 1.
2 donnera 2.
3 donnera 3.
4 donnera 4.
5 donnera 0.
6 donnera 1.
7 donnera 2.

Les nombres entre 0 et 2 auront donc deux chances sur huit de sortir, tandis que les nombres 3 et 4 n'auront qu'une chance sur huit. Le générateur n'est plus uniforme.

Dans cet exemple trivial, pour éviter ce problème de biais dû au modulo, la solution consiste à ignorer les sorties du générateur supérieures à 4. Si jamais on tombe sur une valeur trop grande, on boucle jusqu'à tomber sur une valeur qui convienne ; seulement ensuite, on applique l'opérateur modulo. Bien sûr, en théorie, boucler jusqu'à obtenir un nombre inférieur à la bonne limite peut prendre un temps imprévisible, voire infini ! En pratique, il devrait être tout à fait exceptionnel d'avoir besoin de plus de deux ou trois itérations.

Dans le cas général, il faut ignorer les sorties supérieures à RAND_MAX - (RAND_MAX + 1) % N. Une implémentation possible est donnée ci-dessous :

int x, limit;

limit = RAND_MAX - (RAND_MAX + 1) % N;
do
{
  x = random();
}
while (x > limit);
x %= N;

Attention, ce code naïf ne tient pas compte d'un éventuel débordement de capacité lors du calcul RAND_MAX + 1 ; il faut donc l'adapter selon la plateforme, en choisissant un type d'entier suffisamment large en fonction de la valeur de RAND_MAX.

Notez enfin que si vous utilisez le générateur pseudo-aléatoire arc4random(), ce dernier fournit une fonction arc4random_uniform() qui implémente directement le code nécessaire pour éviter ce biais dû au modulo.