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

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()线程安全吗？也许是，也许不是：

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

一个测试和很好的学习练习是用一个固定的整数替换对rand()的调用，看看会发生什么。

我认为伪随机数生成器就像一个黑盒，它接受一个整数作为输入，并返回一个整数作为输出。对于任何给定的输入，输出总是相同的，但是数字序列中没有模式，并且序列均匀分布在可能的输出范围内。(这个模型并不完全准确，但它可以做到。)使用这个黑盒的方法是选择一个起始数(种子)，在应用程序中使用输出值，并将其作为下一次调用随机数生成器的输入。有两种常见的API设计方法：

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

我通常用Fortran编写，并使用Mersenne Twister PRNG的Ladd's实现(该链接值得一读)。在C中有很多合适的PRNG，它们将状态暴露给你的控制。PRNG看起来不错，使用它(通过并行区域内的初始化和销毁调用以及私有状态变量)应该会给您带来不错的加速。

最后，通常的情况是，如果你一次性请求一个完整的随机数序列，可以让PRNG执行得更好(例如，编译器可以向量化PRNG的内部结构)。正因为如此，库通常有一些类似get_prng_array(state)的函数，它们会给你返回一个充满随机数的数组，就像你把get_prng放在一个循环中，填充数组元素一样--它们只是做得更快。这将是第二次优化(并且需要在并行for循环中添加一个for循环。显然，您不希望在这样做时耗尽每个线程的堆栈空间！