/* compute pi using Monte Carlo method */
#include <math.h>
#include "mpi.h"
#include "mpe.h"
#define CHUNKSIZE      1000
/* We'd like a value that gives the maximum value returned by the function
   random, but no such value is *portable*.  RAND_MAX is available on many 
   systems but is not always the correct value for random (it isn't for 
   Solaris).  The value ((unsigned(1)<<31)-1) is common but not guaranteed */
#define INT_MAX 1000000000

/* message tags */
#define REQUEST  1
#define REPLY    2
int main( int argc, char *argv[] )
{
    int iter;
    int in, out, i, iters, max, ix, iy, ranks[1], done, temp;
    double x, y, Pi, error, epsilon;
    int numprocs, myid, server, totalin, totalout, workerid;
    int rands[CHUNKSIZE], request;
    MPI_Comm world, workers;
    MPI_Group world_group, worker_group;
    MPI_Status status;

    MPI_Init(&argc,&argv);
    world  = MPI_COMM_WORLD;
    MPI_Comm_size(world,&numprocs);
    MPI_Comm_rank(world,&myid);
    server = numprocs-1;	/* last proc is server */
    if (myid == 0)
        sscanf( argv[1], "%lf", &epsilon );
    MPI_Bcast( &epsilon, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD );
    MPI_Comm_group( world, &world_group );
    ranks[0] = server;
    MPI_Group_excl( world_group, 1, ranks, &worker_group );
    MPI_Comm_create( world, worker_group, &workers );
    MPI_Group_free(&worker_group);
    if (myid == server) {	/* I am the rand server */
	do {
	    MPI_Recv(&request, 1, MPI_INT, MPI_ANY_SOURCE, REQUEST,
		     world, &status);
	    if (request) {
		for (i = 0; i < CHUNKSIZE; ) {
		        rands[i] = random();
			if (rands[i] <= INT_MAX) i++;
		}
		MPI_Send(rands, CHUNKSIZE, MPI_INT,
                         status.MPI_SOURCE, REPLY, world);
	    }
	}
	while( request>0 );
    }
    else {			/* I am a worker process */
        request = 1;
	done = in = out = 0;
	max  = INT_MAX;         /* max int, for normalization */
        MPI_Send( &request, 1, MPI_INT, server, REQUEST, world );
        MPI_Comm_rank( workers, &workerid );
	iter = 0;
	while (!done) {
	    iter++;
	    request = 1;
	    MPI_Recv( rands, CHUNKSIZE, MPI_INT, server, REPLY,
		     world, &status );
	    for (i=0; i<CHUNKSIZE-1; ) {
	        x = (((double) rands[i++])/max) * 2 - 1;
		y = (((double) rands[i++])/max) * 2 - 1;
		if (x*x + y*y < 1.0)
		    in++;
		else
		    out++;
	    }
	    MPI_Allreduce(&in, &totalin, 1, MPI_INT, MPI_SUM,
			  workers);
	    MPI_Allreduce(&out, &totalout, 1, MPI_INT, MPI_SUM,
			  workers);
	    Pi = (4.0*totalin)/(totalin + totalout);
	    error = fabs( Pi-3.141592653589793238462643);
	    done = (error < epsilon || (totalin+totalout) > 1000000);
	    request = (done) ? 0 : 1;
	    if (myid == 0) {
		printf( "\rpi = %23.20f", Pi );
		MPI_Send( &request, 1, MPI_INT, server, REQUEST,
			 world );
	    }
	    else {
		if (request)
		    MPI_Send(&request, 1, MPI_INT, server, REQUEST,
			     world);
	    }
	}
	MPI_Comm_free(&workers);
    }

    if (myid == 0) {
        printf( "\npoints: %d\nin: %d, out: %d, <ret> to exit\n",
	       totalin+totalout, totalin, totalout );
	getchar();
    }
    MPI_Finalize();
}
