lundi 1 mars 2021

rand_r() versus rand() in approximation of pi

I am trying to approximate pi for homework by using a Monte-Carlo method, by randomly sampling points in a unit box, and counting the ratio that falls inside a unit circle enclosed by the box. See below. I am asked to do this in parallel using multithreading but decided to get things working first before parallelizing things. My teacher has hence explicitly asked me to use rand_r() for thread safety. I am aware that better pseudo-random number generators exist. However, I cannot seem to get things right, and I figure I am seeding rand_r() the wrong way. I have tried to seed using the computer time, but the value I get is wrong (around 2.8). If I use some random number and seed rand() with srand() instead, I am able to approximate pi pretty easily. Can anyone enlighten me as to what I am missing here? The code I have written looks like this:

#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;
}



Aucun commentaire:

Enregistrer un commentaire