#include <stdio.h>
#include <mpi.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define MIN_EXP 5
#define MAX_EXP 6  // running time is about 1 minute for MAX_EXP=9

int main(int argc,char **argv)
{
  int myid, numprocs, i;
  long long int count, n, j;
  double PI25DT = 3.141592653589793238462643;
  double x, y, mypi, mysigma, pi, sigma, error;

  MPI_Init(&argc,&argv);
  MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
  MPI_Comm_rank(MPI_COMM_WORLD,&myid);

  srand((unsigned)time(NULL)+myid);

  n = 1;
  for (i=1;i<MIN_EXP;i++) n = n*10;
  for (i=MIN_EXP;i<=MAX_EXP;i++)
    {
      n = n*10;
      count = 0;
      for (j=0;j<n;j++)
        {
          x = rand()/((double)RAND_MAX+1);
          y = rand()/((double)RAND_MAX+1);
          if ((x*x+y*y) <= 1.0) count++;
        }
      mypi = (double)(count*4)/(double)n;
      mysigma = pow((mypi-PI25DT),2);
      MPI_Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
      MPI_Reduce(&mysigma, &sigma, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
      if(myid==0)
        {
          pi = pi/(double)numprocs;
          sigma = sqrt(sigma/(double)numprocs);
          error = pi - PI25DT;
          printf("For %d trials with n = %lld:\n",numprocs,n);
          printf("pi= %.12e, error=%.12e, sigma=%.12e\n",pi,error,sigma);
        }
    }
  MPI_Finalize();
  return 0;
}