#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(double start, double delta_x){
    
    return delta_x*f(start+(delta_x/2.))  ;
}

double rectangular_rule(double start, double delta_x){
    
    return delta_x*f(start)    ;
}

double simpsons_rule(double start, double delta_x){
   
    return (delta_x/6.0)*(f(start)+4*f(start+(delta_x/2.))+f(start+delta_x)) ;
}

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

int main(int argc,char *argv[]){
    int done = 0, n, myid, numprocs, mpi_error,namelen,i_counter;
    int buffer[2],local_intervals;
    double local_start;
    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) {
            printf("\nWhich rule do you want to use?\n 0= done \n 1=rectangular rule \n 2=mid point rule\n 3=simpsons rule \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");
                fflush(stdout); //for running on Windows the out buffer needs to be flushed)
                scanf("%d",&buffer[1]);
            }
            startwtime = MPI_Wtime();        
        }//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=(1.0/numprocs)*(myid);
            local_intervals=n/numprocs;
            /******* The next statments are just for diagnostics ********/
            //printf("process=%d local start = %.16f local intervals = %d\n",myid,local_start,local_intervals); 
                  //fflush(stdout);  
            /**********************************************************/
        mypi=0.0;
        i_counter=0;
        for(i_counter=0;i_counter<local_intervals;i_counter++)
                    mypi+=( *rule)(local_start+i_counter*delta_x,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("\npi 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;
}
/*----------------------------------------------------------------*/
/*  Brought To You By CToHTML 
http://www.cs.washington.edu/homes/zahorjan/homepage/Tools/index.htm                     */
/*----------------------------------------------------------------*/