用于蒙特卡罗积分的线程安全随机数生成

Question

我试着写一些非常快速计算随机数的东西,可以应用于多个线程.我目前的代码是:

/* 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);
}

这可以工作,但从观察资源管理器我知道它不使用所有线程.rand()是否适用于多线程代码？如果不是,那还有一个好的选择吗？非常感谢.插口

Answer 1

是rand()线程安全的？也许,也许不是:

rand()函数不需要是可重入的.不需要可重入的函数不需要是线程安全的."

一个测试和良好的学习练习将替换调用rand(),例如,一个固定的整数,看看会发生什么.

The way I think of pseudo-random number generators is as a black box which take an integer as input and return an integer as output. For any given input the output is always the same, but there is no pattern in the sequence of numbers and the sequence is uniformly distributed over the range of possible outputs. (This model isn't entirely accurate, but it'll do.) The way you use this black box is to choose a staring number (the seed) use the output value in your application and as the input for the next call to the random number generator. There are two common approaches to designing an API:

1. 两个函数,一个用于设置初始种子(例如srand(seed)),另一个用于从序列中检索下一个值(例如rand()).PRNG的状态存储在内部的全局变量中.生成一个新的随机数要么不是线程安全的(很难说,但是输出流不可重复),要么在多线程代码中会很慢(最终会在状态值周围进行一些序列化).
2. PRNG状态向应用程序编程人员公开的接口.在这里,您通常有三个函数:init_prng(seed)返回PRNG状态的一些不透明表示get_prng(state),它返回一个随机数并更改状态变量,并且destroy_peng(state)只清除已分配的内存等等.具有此类API的PRNG应该都是线程安全的并且没有锁定并行运行(因为您负责管理(现在是线程本地)状态变量.

我通常用Fortran编写并使用Ladd实现的Mersenne Twister PRNG(该链接值得一读).在C中有许多合适的PRNG,它们将状态暴露给你的控制.PRNG看起来很好并且使用它(在并行区域和私有状态变量中初始化和销毁​​调用)应该给你一个不错的加速.

最后,如果您一次性要求完整的随机数序列(例如,编译器可以对PRNG内部进行矢量化),通常情况下PRNG可以表现得更好.因为这个库经常会有类似get_prng_array(state)函数的东西,它会让你返回一个充满随机数的数组,好像你放入get_prng一个填充数组元素的循环 - 它们只是更快地完成它.这将是第二次优化(并且需要在并行for循环中添加for循环.显然,您不希望用完每个线程的堆栈空间!