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

// kernel
double integrate(int n, int myid, int numprocs) {
  double h, x, sum;
  int i;

  sum = 0.0;
  h   = 1.0 / (double) n;
  for (i = myid + 1; i <= n; i += numprocs) {
    x = h * ((double)i - 0.5);
    sum += 4.0 / (1.0 + x*x);
  }
  return sum * h;
}

int main(int argc, char *argv[]) {
  int n, myid, numprocs;
  double PI25DT = 3.141592653589793238462643;
  double mypi, pi;
  double startwtime, endwtime;
  int  namelen;
  char processor_name[MPI_MAX_PROCESSOR_NAME];
  
  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("MPI Task %2d on %20s\n", myid, processor_name);
  sleep(1); // wait for everyone to print
  while (1) {
    if (myid == 0) {
      printf("Enter the number of intervals: (0 quits) ");fflush(stdout);
      scanf("%d",&n);

      startwtime = MPI_Wtime();
    }
    MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);
    if (n == 0)
      break;
    mypi = integrate(n, myid, numprocs);
    
    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 = %.6f\n",
	     endwtime-startwtime);	       
    }
  }
  MPI_Finalize();

  return 0;
}