Génération de nombres aléatoires thread-safe pour l'intégration de Monte-Carlo

Question

J'essaie d'écrire quelque chose qui calcule très rapidement des nombres aléatoires et peut être appliqué sur plusieurs threads. Mon code actuel est:

/* Approximating PI using a Monte-Carlo method. */ 

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <time.h> 
#include <omp.h> 
#define N 1000000000 /* As lareg as possible for increased accuracy */ 

double random_function(void); 

int main(void) 
{ 
    int i = 0; 
    double X, Y; 
    double count_inside_temp = 0.0, count_inside = 0.0; 
    unsigned int th_id = omp_get_thread_num(); 
    #pragma omp parallel private(i, X, Y) firstprivate(count_inside_temp) 
    { 
     srand(th_id); 
     #pragma omp for schedule(static) 
     for (i = 0; i <= N; i++) { 
     X = 2.0 * random_function() - 1.0; 
     Y = 2.0 * random_function() - 1.0; 
     if ((X * X) + (Y * Y) < 1.0) { 
     count_inside_temp += 1.0; 
    } 
    } 
    #pragma omp atomic 
     count_inside += count_inside_temp; 
    } 
    printf("Approximation to PI is = %.10lf\n", (count_inside * 4.0)/ N); 
    return 0; 
} 

double random_function(void) 
{ 
    return ((double) rand()/(double) RAND_MAX); 
} 

Cela fonctionne, mais d'observer un gestionnaire de ressources Je sais que ne pas utiliser tous les fils. Est-ce que rand() fonctionne pour du code multithread? Et sinon, y a-t-il une bonne alternative? Merci beaucoup. Jack

Answer 1

Le rand() est-il sécuritaire? Peut-être, peut-être pas:

The rand() function need not be reentrant. A function that is not required to be reentrant is not required to be thread-safe."

Un test et un bon exercice d'apprentissage serait de remplacer l'appel à rand() avec, disons, un entier fixe et voir ce qui se passe. La façon dont je pense aux générateurs de nombres pseudo-aléatoires est comme une boîte noire qui prend un entier en entrée et renvoie un entier en sortie. Pour toute entrée donnée, la sortie est toujours la même, mais il n'y a pas de motif dans la séquence de nombres et la séquence est uniformément répartie sur la plage des sorties possibles. (Ce modèle n'est pas tout à fait exact, mais il le fera.) La façon dont vous utilisez cette boîte noire est de choisir un nombre fixe (la graine) utiliser la valeur de sortie dans votre application et comme entrée pour le prochain appel à la générateur de nombres aléatoires. Il existe deux approches communes pour la conception d'une API:

1. Deux fonctions, une pour définir la graine initiale (par exemple srand(seed)) et un pour récupérer la valeur suivante de la séquence (par exemple rand()). L'état du PRNG est stocké en interne dans une sorte de variable globale. Générer un nouveau nombre aléatoire ne sera pas thread-safe (difficile à dire, mais le flux de sortie ne sera pas reproductible) ou sera lent dans le code multithreded (vous vous retrouvez avec une sérialisation autour de la valeur d'état).
2. Une interface où l'état PRNG est exposé au programmeur d'application. Ici vous avez généralement trois fonctions: init_prng(seed), qui renvoie une représentation opaque de l'état PRNG, get_prng(state), qui renvoie un nombre aléatoire et modifie la variable d'état, et destroy_peng(state), qui nettoie juste la mémoire allouée et ainsi de suite. PRNGs avec ce type d'API doivent tous être thread-safe et fonctionnent en parallèle sans verrouillage (parce que vous êtes en charge de la gestion du (fil maintenant locale) variable d'état.

J'écris généralement Fortran et utiliser Ladd's mise en œuvre du PRNG de Mersenne Twister (ce lien vaut la peine d'être lu) Il y a beaucoup de PRNG appropriés en C qui exposent l'état à votre contrôle PRNG ça a l'air bien et ça (avec l'initialisation et la destruction des appels dans la région parallèle et les variables privées) Enfin, il est souvent possible que les PRNG soient plus performants si vous demandez une séquence entière de nombres aléatoires en une seule fois (e). .g. le compilateur peut vectoriser les internes PRNG). A cause de cela, les bibliothèques ont souvent quelque chose comme get_prng_array(state) fonctions qui vous renvoie un tableau plein de nombres aléatoires comme si vous mettez get_prng dans une boucle remplissant les éléments du tableau - ils le font simplement plus rapidement. Ceci serait une deuxième optimisation (et aurait besoin d'une boucle for ajoutée à l'intérieur de la boucle parallèle.) Evidemment, vous ne voulez pas manquer d'espace de pile par threads!