Faire de l'aléatoire convenablement

Bonjour à tous !

Dans le cadre d'un programme (que je dois absolument coder en c), j'aimerai reproduire le comportement de la fonction random.randrange(), qui permet de retourner un nombre $x$ tel que $ x \in [\text{min}; \text{max}]$. Ça veux dire aussi que $x \in \mathbb{R}$.

Une solution comme suis n'est donc pas envisageable :

1
2
3

double randrange(double min, double max) {
    return min + rand() % (max-min);
}

Déjà parce que le modulo donne un nombre entier, ensuite parce qu'il est dit un peu partout qu'il ne faut pas faire comme ça (j'ai pas trop pigé pourquoi). Du coup, sur les conseils de cette page, j'ai écris ceci :

1
2
3
4
5
6
7
8
9

double randrange(double M, double N) {
    if (M > N) {
        double T = M;
        M = N;
        N = T;
    }

    return M + rand()/(RAND_MAX + 1.) * (N-M);
}

Mais je ne suis pas sur qu'elle soie meilleure (déjà, pourquoi + 1. ?!?). En plus, pour l'application que j'écris, j'ai besoin de "vrai" aléatoire ou de ce qui s'en rapproche de plus (donc pas une distribution gaussienne), et la plupart des messages sur internet parlent peu ou prou de "casser l'aléatoire" quand on écrit une fonction du genre.

J'ai aussi lu à propos de drand48(), mais je suis pas sous OpenBSD.

Est ce que ma fonction est correcte ? Est ce qu'il y a des manières plus intelligentes de faire ? D'avance merci

Tu peux aussi regarder le code de randrange() si tu veux avoir le même comportement.

Comme je sais que la fonction random() de python est une fonction uniforme, j'ai donc été voir comment ils faisaient et voici ce que me dit le code de cpython :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

/* generates a random number on [0,0xffffffff]-interval */

static PY_UINT32_T
genrand_int32(RandomObject *self)
{
    PY_UINT32_T y;
    static const PY_UINT32_T mag01[2] = {0x0U, MATRIX_A};
    /* mag01[x] = x * MATRIX_A  for x=0,1 */
    PY_UINT32_T *mt;
    mt = self->state;
    if (self->index >= N) { /* generate N words at one time */
        int kk;

        for (kk=0;kk<N-M;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1U];
        }
        for (;kk<N-1;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1U];
        }
        y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
        mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1U];

        self->index = 0;
    }

    y = mt[self->index++];
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680U;
    y ^= (y << 15) & 0xefc60000U;
    y ^= (y >> 18);
    return y;
}

j'espère que ça t'intéresse.

Si le but est de ne pas se prendre trop le chou et que tu es sous un système Unix, lire le fichier /dev/random devrait suffire.

J'ai écris ça un peu vite ce matin, je vois

Déjà, ce que je craignais est arrivé :

C'est bel et bien un algo de type Monte Carlo que j'essaye d'écrire. Je me doutais bien qu'il y avait un truc du genre qui allait sortir, mais les papiers qui décrivent Monte Carlo partent du principe que tu a un Bac+18 en statistique et que ce genre de chose est obvious. Merci, donc

@Kje et @artragis : merci. Sauf que je me suis rendu compte que randrange() reste dans le domaine des entiers, or je veux un réel Après, le code de génération d'entier de cpython m'as l'air de correspondre à ce que Gabbro mentionnait comme étant "Mersenne Twister", donc je peux éventuellement leur emprunter … Mais faut que je sois sur que ma fonction randrange() du dessus me détruise pas la qualité de l'aléatoire généré.

@Renault: de fait, et à priori, je resterai sur de l'UNIX tout du long. Néanmoins, on ne sais jamais (question de rapidité, /dev/urandom est pas "lent" vu qu'il y a une lecture à faire ?).

Je ne comprens pas très bien ton message : a priori, ce que tu demande correspond à random.uniform() de python et non à randrange.

Pour la recoder, comme pour beaucoup d'autres fonctions aléatoires, il vaut mieux avoir d'abord l'équivalent de random.random, qui elle est relativement simple à coder en C (mais note que la granularité va dépendre de RAND_MAX, et la qualité va dépendre de l'implémentation de rand()) :

1
2
3

double random() {
    return (double) rand() / RAND_MAX;
}

Ensuite, tu peux définir uniform comme suit (correspond à la définition de python) :

1
2
3

double uniform(double a, double b) {
    return a + (b-a) * random();
}

EDIT:

@gabbro: Il me semble que la qualité ne se dégrade pas, mais l'aléatoire devient en quelque sorte cyclique sur de très longues chaines. Typiquement, ça ne dérange en général pas vraiment pour du Monte-Carlo.

@yoch: de fait. Néanmoins, c'est à peut près ce que je fais en une fonction, non?

Oui, quasiment, je n'avais pas vu. Il te manque toutefois le cast en double (très important, et vaut mieux le faire explicitement, je ne saurais dire si le tien est transtypé au bon moment dans ton code), et je ne comprend pas bien pourquoi tu divise par RAND_MAX + 1 ? EDIT: je viens de voir que le + 1 vient de la page de comp.lang.c que tu met en lien, or c'est censé être des experts là bas. Du coup, as-tu regardé cette page ?

Je ne connais pas bien les différents algorithmes de génération de nombres aléatoires, mais /dev/{u}random est un pseudo-fichier. Il te permet d'accéder au kernel globalement. Donc il n'y a pas de lecture disque-dur à proprement parler. Plus d'infos avec man urandom.

@QuentinC : Merci pour ces explications.

Toutefois, ce que tu avance pour le +1 ne me satisfait pas. En effet, RAND_MAX fait partie de l'intervalle (dans ton cas par exemple il vaudra 15). Et on ne va donc jamais dépasser la borne (mais l'atteindre si).

Petit test :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

>>> RAND_MAX = 15
>>> for i in range(RAND_MAX + 1):
    print(i, '=>', 1 + (i / RAND_MAX) * (10-1))


0 => 1.0
1 => 1.6
2 => 2.2
3 => 2.8
4 => 3.4
5 => 4.0
6 => 4.6
7 => 5.2
8 => 5.8
9 => 6.3999999999999995
10 => 7.0
11 => 7.6
12 => 8.2
13 => 8.8
14 => 9.4
15 => 10.0

Je retiens donc plutôt l'explication que l'on veut éviter d'avoir la borne supérieure dans les résultats possible, ce qui n'est pas le cas pour random.uniform.

En ce qui concerne /dev/urandom, il me semble qu'il est très déconseillé de l'utiliser de façon massive, et je suppose que c'est lent de toutes façons.

1 2 3

>> RAND_MAX = 15 >> for i in range(RAND_MAX + 1): print(i, '=>', 1 + (i / RAND_MAX) * (10-1))

Tu ne démontres pas la même chose que moi. Soit tu fixes RAND_MAX=16 et tu boucles jusqu'à RAND_MAX et non RAND_MAX+1, soit tu divises par RAND_MAX+1 et non RAND_MAX. IL faut choisir. C'est évident que si on prend x=16 dans le calcul 1+(x/16)*10, on tombe pile sur 11.

Si RAND_MAX vaut usuellement une puissance de 2 comme 65536 ou 32768, c'est justement parce qu'elle n'est jamais atteinte par l'algorithme de congruence linéaire. Pour rappel l'algorithme de base en gros résumé c'est nombre_suivant = (nombre_actuel * A + B) % RAND_MAX avec A et B des constantes définies par l'implémentation. Sauf erreur, A et B sont censés être proche de sqrt(RAND_MAX) et premiers entre eux et avec RAND_MAX si on veut un résultat optimal: Je n'ai pas le code de la libc mais on dit partout que ces constantes A et B n'ont justement pas été bien choisis.

Merci à tous, c'est très intéressant. Si je résume, ma fonction pour avoir un range est "correcte", mais le générateur de nombre aléatoire est à ch***

Je vais probablement emprunter le code de Mersenne/Twister (merci fromvega)

@Renault:

Yep, mais la page wikipédia précise que sous Linux, l'appel à /dev/random est bloquant. À priori, j'ai besoin de générer une quantité importante de nombre aléatoire à la seconde (bon, pas des millions, mais une centaine dans le meilleur des cas), et j'ai peur que la sortie ne suive pas¹. D'autant qu'à terme, c'est pour faire tourner sur un super-calculateur et plus sur ma machine, donc je ne peut pas prévoir le temps de réaction. Mais c'est une alternative.

@Yoch: en fait, je me suis trompé sur un point. Je viens de tester rapidement (j'aurais dû le faire avant) et chez moi avec GCC 4.7 windows, RAND_MAX = 32767, et non pas 32768 comme je l'avais cru. Donc ta démonstration est correcte mais ja'vais quand même vu juste, le +1 sert bien à s'assurer qu'on ne tombera jamais sur la borne supérieure; sauf qu'il s'avère indispensable du coup. ET hop on a l'explication définitive sur pourquoi il est là.

Tu fais erreur sur ce point. Comme déjà dit, RAND_MAX fait bien partie de l'intervalle, et sa valeur est usuellement de la forme $2^n -1$ (l'implémentation étant alors nombre_suivant = (nombre_actuel * A + B) & RAND_MAX).

EDIT: bon bah grilled, tant pis je laisse

C'est forcément en C que tu dois écrire ton code ? En C++ Mersenne-Twister est standard depuis C++11, en Fortran tu as random_number avec un vrai générateur de nombrs alératoires dans gfortran et ifort. Et avec rust, c'est dans une lib sympa !

Si tu n'as pas le choix, j'ai collecté quelques codes Monte-Carlo et j'en ai au moins un en C, avec une implémentation de générateur de nombres aléatoires. je dois pouvoir t'en passer une partie, à vérifier avec la licence.

C'est forcément en C que tu dois écrire ton code ? En C++ Mersenne-Twister est standard depuis C++11, en Fortran tu as random_number avec un vrai générateur de nombrs alératoires dans gfortran et ifort. Et avec rust, c'est dans une lib sympa !

Je sais pas coder en C++ (je fais du C with class, donc autant faire du C), j'aime pas le Fortran (oui, je sais, tout les codes scientifiques de l'univers sont écrit en Fortran¹, je lis le Fortran et j'aime pas ça) et je connais pas le Rust. Et je savais très bien que quelqu'un allais encore me sortir le C++11 et son (ses?) super-générateur-de-nombre-aléatoire-de-la-mort-qui-tue-disponible-en-25-milliards-de-possibilités-que-tu-sais-pas-quoi-choisir-et-les-défenseurs-de-C++-me-tapent-sur-le-système-avec-ça.²

C'est le genre de commentaires qui m'ennuient profondément. Je fais du C si je veux, et c'est pas la question. Sinon, j'aurais pas posé la question et j'aurais fait du Python (pas pour rien que je cherche l'équivalent C de la fonction Python correspondante).

Edit @Kass'Peuk: Je tape parce que on m'as encore sorti le classique "change de langage, le mien est mieux, t'as vu". Le reste est de la mauvaise foi tout à fait assumée.