Centre Blaise Pascal » Bench4GPU

//

// Estimation of Pi using Monte Carlo exploration process

// gcc -std=c99 -O3 -o Pi Pi.c -lm 

//

#include <math.h>

#include <stdio.h>

#include <stdlib.h>

#include <omp.h>

// Marsaglia RNG very simple implementation

#define znew  ((z=36969*(z&65535)+(z>>16))<<16)

#define wnew  ((w=18000*(w&65535)+(w>>16))&65535)

#define MWC   (znew+wnew)

#define SHR3  (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))

#define CONG  (jcong=69069*jcong+1234567)

#define KISS  ((MWC^CONG)+SHR3)

#define MWCfp MWC * 2.328306435454494e-10f

#define KISSfp KISS * 2.328306435454494e-10f

#define ITERATIONS 1000000000

#define PROCESS 8

#ifdef LONG

#define LENGTH unsigned long

#else

#define LENGTH unsigned int

#endif

LENGTH MainLoopGlobal(LENGTH iterations,unsigned int seed_w,unsigned int seed_z)

{

   unsigned int z=seed_z;

   unsigned int w=seed_w;

   LENGTH total=0;

   for (LENGTH i=0;i<iterations;i++) {

      float x=MWCfp ;

      float y=MWCfp ;

      // Matching test

      int inside=((x*x+y*y) < 1.0f) ? 1:0;

      total+=inside;

   }

   return(total);

}

int main(int argc, char *argv[]) {

  unsigned int seed_w=10,seed_z=10;

  LENGTH iterations=ITERATIONS;

  LENGTH inside[1024],insides=0;

  if (argc > 1) {

    iterations=(LENGTH)atol(argv[1]);

  }

  else {

    printf("\n\tPi : Estimate Pi with Monte Carlo exploration\n\n\t\t#1 : number of iterations (default 1 billion)\n\n");

  }

#pragma omp parallel for

  for (int i=0 ; i<PROCESS; i++) {

    inside[i]=MainLoopGlobal(iterations/PROCESS,seed_w,seed_z);

    printf("%lu\n",inside[i]);

  }

  for (int i=0 ; i<PROCESS; i++) {

    insides+=inside[i];

  }

  float pi=4.*(float)insides/(float)iterations;

  printf("\tPi=%f with error %f and %lu iterations\n\n",pi,

         fabs(pi-4*atan(1))/pi,(unsigned long)iterations);

}