#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <signal.h>
#include "mpi.h"

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

double mid_point_rule(int start,int intervals, double delta_x){
    double area = 0.0; 
    int i;
    for (i = start; i <= start+intervals; i++) {
        area += delta_x * f((double) (delta_x*i+delta_x/2.));
    }
    return area  ;
}

double rectangular_rule(int start,int intervals, double delta_x){
    double area = 0.0; 
    int i;
    for (i = start; i <= start+intervals; i++) {
        area += delta_x * f((double) delta_x*i);
    }
    return area  ;
}

double simpsons_rule(int start,int intervals, double delta_x){
    double area = f(delta_x*start)-f(delta_x*(start+intervals)); 
    int i;
    for (i = 1; i <= (intervals/2); i++) {
        area += 4.0*f(delta_x*(start+ (2*i-1)))+2.0*f(delta_x*(start+ 2*i));
    }
    return delta_x*area/3.0 ;
}

double (*rule)(int start,int intervals, double delta_x);

int main(int argc,char *argv[]){
    int done = 0, n, myid, numprocs, mpi_error,namelen;
    int buffer[2],local_start,local_intervals;
    double PI_to_25dp = 3.141592653589793238462643;
    double mypi, pi, delta_x;
    double startwtime, endwtime,beginwtime;
    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("Process number %d on processor %s \n", myid, processor_name, numprocs);
    fflush(stdout);
    n = 0;
    while (!done) {
        if (myid == 0) {
            beginwtime = MPI_Wtime();
            printf("\nWhich rule do you want to use?\n 0= done \n 1=rectangular rule \n 2=mid point rule\n 3=simpsons rule (choose intervals/process(or)s to be even)\n ");
            fflush(stdout);
            scanf("%d",&buffer[0]);
            if ((buffer[0]<1)||(3<buffer[0])) {
                done=1;
                buffer[0]=0;
            }
            if (!done) {
                printf("\nPlease input the number of intervals(0 quits).\n (Note: In general, the error is minimized if the number of intervals is a multiple of the number of processors)\n");
                fflush(stdout); //for running on Windows the out buffer needs to be flushed)
                scanf("%d",&buffer[1]);
            }
            startwtime = MPI_Wtime();
            printf("wall clock begin = %f\n",beginwtime); 
            fflush(stdout);
        }//end setup by the master process..0
        MPI_Bcast(buffer, 2, MPI_INT, 0, MPI_COMM_WORLD);
        if (buffer[0] == 0)
            done = 1;
        else {
            n=buffer[1];
            switch (buffer[0]) {
            case 0:
                done=1;
                break;
            case 1:
                rule=&rectangular_rule;
                break;
            case 2:
                rule=&mid_point_rule;
                break;
            case 3:
                rule=&simpsons_rule;
                break;
            }
            /*******set up this processes part of the calculation********/
            delta_x   = 1.0 / (double) n;
            local_start=(n/numprocs)*(myid);
            local_intervals=n/numprocs;
            /******* The next statments are just for diagnostics ********/
            printf("process=%d local start = %d local intervals = %d\n",myid,local_start,local_intervals); 
            fflush(stdout);
            /**********************************************************/
            mypi=( *rule)(local_start,local_intervals,delta_x);
            /*  The next statment---either sends mypi back to 0 or if 0 does the calculation*/
            MPI_Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

            if (myid == 0) {
                printf("pi is approximately %.16f, Relative Error is %16.8e\n",
                       pi, (double)100 * (pi - PI_to_25dp)/PI_to_25dp);
                fflush(stdout);
                endwtime = MPI_Wtime();
                printf("wall clock time = %f\n", endwtime-startwtime); 
                fflush(stdout);
            }
        }
    }
    mpi_error=MPI_Finalize();
    if (MPI_SUCCESS!=mpi_error)
        return mpi_error;
    else
        return 0;
}