rand_r() versus rand() in approximation of pi

Question

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

Answer 1

You're not updating seed each time through the loop.

rand_r() updates the local variable seed in the randomNumber() function, but that doesn't affect the variable in main(). The result is that you're generating the same random number every time you call it.

You need to pass a pointer to main's variable to 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;
}