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

#include <mpi.h>

#include <sys/time.h>

/* internal variables and functions for timer */
static double get_time(void)
{
  struct timeval tv;
  struct timezone tz;

  gettimeofday(&tv,&tz);
  double systime = tv.tv_sec;
  systime += ((double) tv.tv_usec)/1000000.0;
  return systime;
}

static double func(double x)
{
    return 4.0/(1.0+x*x);
}

int main( int argc, char *argv[] )
{
    const double pi_known = 3.141592653589793238462643;

    /* 
     * mypi      - the area under the curve calculated by a single process
     * sumf      - intermediate sum of slices
     * x         - x value at the center of a slice
     * tstart/tstop - for timing measurements
     * argstr    - command line argument string
     * n         - number of rectangless in which to divide the curve
     * i         - loop counter on a given node
     */

    double mypi,h,sumf,x,tstart,tstop;
    int n, i, ncpu, me;

    n = MPI_Init(&argc,&argv);
    if (n != MPI_SUCCESS) {
        puts("problem setting up the MPI system");
        return n;
    }

    MPI_Comm_size(MPI_COMM_WORLD,&ncpu);
    MPI_Comm_rank(MPI_COMM_WORLD,&me); 
      
    if (argc < 2) {
        if (me == 0 )
            fprintf(stderr, "Usage: %s <intervals> "
                    "(with intervals > number processes)\n", argv[0]);

        MPI_Finalize();
        return 1;
    }

    n=atoi(argv[1]);
    if (n < ncpu) {
        if (me == 0)
            fprintf(stderr, "Usage: %s <intervals> "
                    "(with intervals > number procs / %d)\n",argv[0],ncpu);
        MPI_Finalize();
        return 1;
    }

    if (me == 0) printf("Integration with %d intervals\n", n);

    /* get time stamp at starting point */
    tstart = get_time();

    /* compute part */
    h    = 1.0 / (double) n;
    sumf = 0.0;
    for (i = me; i < n; i += ncpu) {
        x = h * ((double)i + 0.5);
        sumf += func(x);
    }
    sumf *= h;
    MPI_Reduce(&sumf,&mypi,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);

    /* get time stamp at completion */
    tstop = get_time();

    if (me == 0) {
        printf("pi estimate= %.16f, Error is %.16f\n",
               mypi, fabs(mypi - pi_known));

        printf("time used: %g\n", tstop-tstart);
    }

    MPI_Finalize();
    return 0;
}