rand_r() 与 rand() 在 pi 的近似值

rand_r() versus rand() in approximation of pi

我正在尝试通过使用蒙特卡洛方法,通过在单位框中随机采样点,并计算落在由框包围的单位圆内的比率来近似 pi 的家庭作业。见下文。我被要求使用多线程并行执行此操作,但我决定在并行化之前先让事情正常工作。因此,我的老师明确要求我使用 rand_r() 来保证线程安全。我知道存在更好的伪随机数生成器。但是,我似乎无法把事情做好,而且我认为我以错误的方式播种 rand_r() 。我曾尝试使用计算机时间播种,但我得到的值是错误的(大约 2.8)。如果我改用一些随机数和种子 rand() 和 srand(),我可以很容易地近似 pi。谁能告诉我我在这里缺少什么?我编写的代码如下所示:

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>

double randomNumber( unsigned int seed){
  /*_Thread_local*/
  double maxRand      =   (double)RAND_MAX;           // Maximum random number, cast to double
  double randNum      =   (double)rand_r( &seed );    // Generate pseudo-random number from seed, cast to double

  return 2 * randNum/maxRand - 1;                     // Recast number between -1 and 1
}

int main( void ){

  unsigned int seed         =   time(NULL);

  int numOfPts              =   (int)1e8  ;
  int ptsInCircle           =   0         ;
  double unitCircleRadius   =   1.0       ;

  double xpos   =   0;
  double ypos   =   0;

  for ( int iteration = 0; iteration < numOfPts; iteration++ ){
    xpos = randomNumber(seed);
    ypos = randomNumber(seed);

    if ( sqrt( pow(xpos, 2) + pow(ypos, 2) ) <= unitCircleRadius ){
      ptsInCircle++;
    }

  }

  double myPiApprox = 4.0*((double)ptsInCircle)/((double)numOfPts);

  printf("My approximation of pi = %g\n", myPiApprox);

  return 0;
}

你是:

  1. 在每次迭代时重新播种随机数生成器而不是让它做它的事情,并且

  2. 您不仅为每次迭代重新播种,而且为 xy 使用相同的种子,因此您仅在对角线上创建元素(其中 x == y)

    编辑:在您使用 randomNumber(iteration) 的原始问题中,它会在对角线上生成排序随机数,在您的编辑中,您将其更改为 randomNumber(seed),它始终是完全相同的值.

更改randomNumber的定义:

double randomNumber(unsigned int* seedp) {
  // ...
  double randNum      =   (double)rand_r(seedp);
  // ...
}

和你的循环:

for ( int iteration = 0; iteration < numOfPts; iteration++ )
  xpos = randomNumber(&seed);
  ypos = randomNumber(&seed);
  // ...
}

您不会在循环中每次都更新 seed

rand_r() 更新 randomNumber() 函数中的局部变量 seed,但这不会影响 main() 中的变量。结果是您每次调用它时都会生成相同的随机数。

您需要将指向 main 的变量的指针传递给 randomNumber()

double randomNumber( unsigned int *seed){
  /*_Thread_local*/
  double maxRand      =   (double)RAND_MAX;           // Maximum random number, cast to double
  double randNum      =   (double)rand_r(seed);       // Generate pseudo-random number from seed, cast to double

  return 2 * randNum/maxRand - 1;                     // Recast number between -1 and 1
}

int main( void ){
  unsigned int seed         =   time(NULL);
  int numOfPts              =   (int)1e8  ;
  int ptsInCircle           =   0         ;
  double unitCircleRadius   =   1.0       ;
  double xpos   =   0;
  double ypos   =   0;

  for ( int iteration = 0; iteration < numOfPts; iteration++ ){
    xpos = randomNumber(&seed);
    ypos = randomNumber(&seed);
    if ( sqrt( pow(xpos, 2) + pow(ypos, 2) ) <= unitCircleRadius ){
      ptsInCircle++;
    }
  }

  double myPiApprox = 4.0*((double)ptsInCircle)/((double)numOfPts);
  printf("My approximation of pi = %g\n", myPiApprox);

  return 0;
}