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


int main(int argc, char *argv[]) 
{
  int n, n0, nn, np, me, ier, i;
  double pi, pi0, x, y;
  MPI_Init(&argc,&argv);
  if ( argc != 2 ) {
    printf("usage: mpirun -np X a.out N_points \n");
    MPI_Finalize();
    exit(1);
  }
  n0 = atoi(argv[1]);
  MPI_Comm_size(MPI_COMM_WORLD,&np);
  MPI_Comm_rank(MPI_COMM_WORLD,&me);
  nn = n0/np;
  n = nn*np;
  printf("me=%d n=%d nn=%d\n",me,n,nn);
  srand(time(NULL)+me);
  pi = 0.0;
  for ( i=0; i<nn; i++ ) {
    x = (double)rand()/RAND_MAX;
    y = (double)rand()/RAND_MAX;
    if ( sqrt( x*x + y*y ) <= 1 )
      pi = pi + 1;
  }
  MPI_Reduce(&pi,&pi0,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
  if ( me==0 ) 
    printf("Monte Carlo approximation to pi= %lf\n",4*pi0/n);
  MPI_Finalize();
}