/*	compute pi by integrating f(x) = 4/(1 + x**2)     
 *
 *	each process:
 *		- receives the # of intervals used in the approximation
 *		- calculates the areas of it's rectangles
 *		- synchronizes for a global summation
 *	process 0 prints the result and the time it took
 */

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

/*
 *	the integration function
 */
double
f(a)

double			a;

{
	return(4.0 / (1.0 + a*a));
}


main(argc, argv)

int			argc;
char			*argv[];

{
	int		myid, numprocs;
	int		n, i;
	int		namelen;
	char		processor_name[MPI_MAX_PROCESSOR_NAME];
	double		startwtime, endwtime;
	double		PI25DT = 3.141592653589793238462643;
	double		mypi, pi;
	double		h, sum, x;

	MPI_Init(&argc, &argv);

	MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
	MPI_Comm_rank(MPI_COMM_WORLD, &myid);

	MPI_Get_processor_name(processor_name, &namelen);
	printf("Process %d on %s\n", myid, processor_name);
/*
 * Assumes that only process 0 originally knows the # intervals.
 * It broadcasts it to the others and times the application.
 */
	if (myid == 0) {
		n = 100;
		startwtime = MPI_Wtime();
	}

	MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);
/*
 * Each process computes its contribution.
 */
	h = 1.0 / (double) n;
	sum = 0.0;

	for (i = myid + 1; i <= n; i += numprocs) {
		x = h * ((double) i - 0.5);
		sum += f(x);
	}

	mypi = h * sum;
/*
 * The final result is available at process 0.
 */
	MPI_Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

	if (myid == 0) {
		endwtime = MPI_Wtime();

		printf("pi is approximately %.16f, Error is %.16f\n",
			pi, fabs(pi - PI25DT));
		printf("wall clock time = %f\n", endwtime - startwtime);
	}

	MPI_Finalize();
	return 0;
}