#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#include <mpi.h>

const bool debug = false;
const double a = 0.0;
const double b = 1.0;

/* Function declarations */
void readinput(int argc, char* argv[], int iam, double* n_p);
double trapezoidal_rule(double left_endpt, double right_endpt, int trap_count,
   double base_len);
double f(double x);

int main(int argc, char** argv) {
   int iam, np, local_n;
   double n, h, local_a, local_b;
   double local_int, total_int;
   double start, finish, loc_elapsed, elapsed;

   /* Star the MPI environment and discover who we are */
   MPI_Init(NULL,NULL);
   MPI_Comm_rank(MPI_COMM_WORLD, &iam);
   MPI_Comm_size(MPI_COMM_WORLD, &np);

   readinput(argc, argv, iam, &n); // Read the user input from the command line

   /*Note: h and local_n are the same for all processes (load balance!)*/
   h = (b-a)/n;     /* basis of each trapezoid */
   local_n = n/np;  /* number of trapezoids per process */

   /* Length of each process' interval of integration = local_n*h. */
   local_a = a + iam*local_n*h;
   local_b = local_a + local_n*h;
   if (debug) printf("Process %d a = %f b = %f\n",iam,local_a,local_b);
   /* Start a timer */
   MPI_Barrier(MPI_COMM_WORLD);
   start = MPI_Wtime();
   /* We compute the local integral on each process using the local endpoints*/
   local_int = trapezoidal_rule(local_a, local_b, local_n, h);
   if (debug) printf("Process %d local integral = %f\n",iam,local_int);
   finish = MPI_Wtime();
   loc_elapsed = finish-start;
   /* Do a max reduce to get the time */
   MPI_Reduce(&loc_elapsed, &elapsed, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);

   /* Do a sum reduce of the integrals computed on each process */
   MPI_Reduce(&local_int, &total_int, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);


   if (iam == 0) {
      printf("Employing n = %.0f trapezoids\n", n);
      printf("Trapezoid per process n = %d trapezoids\n", local_n);
      printf("Number of processes: %d\n",np);
      printf("of the integral of f(x) from %.0f to %.0f = %1.16f\n", a, b, total_int);
      printf("Total elapsed time = %f milliseconds \n", elapsed * 1000);
   }

   /* Shut down MPI */
   MPI_Finalize();

   return 0;
}

void readinput(int argc, char* argv[], int iam, double* n_p){
/*READINPUT gets the user input: the number of trapezoids
 * Input args:   iam:  process rank in MPI_COMM_WORLD
 *               np:  number of processes in MPI_COMM_WORLD
 * Output args:  n_p:  pointer to number of trapezoids
 */
	if (iam == 0) {
		if (argc!= 2){
		    fprintf(stderr, "usage: srun %s <number of trapezoids> \n", argv[0]);
            fflush(stderr);
            *n_p = -1;
		} else {
			*n_p = atoi(argv[1]);
		}
	}
	/* Broadcasts value of n to each process */
	MPI_Bcast(n_p, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	
	// negative n ends the program
    if (*n_p <= 0) {
        MPI_Finalize();
        exit(-1); // we return an error!
    }
}

double trapezoidal_rule(double left_endpt, double right_endpt, int trap_count, double base_len) {
/* TRAPEZOIDAL_RULE the serial function applying the trapezoidal rule
 * Input:   left_endpt
 *          right_endpt
 *          trap_count
 *          base_len
 * Output:  Approximated value of the integral   
 */
   double estimate, x;
   int i;

   estimate = (f(left_endpt) + f(right_endpt))/2.0;
   for (i = 1; i <= trap_count-1; i++) {
      x = left_endpt + i*base_len;
      estimate += f(x);
   }
   estimate = estimate*base_len;

   return estimate;
}


double f(double x) {
/*F function to be integrated */
   return (4.0*sqrt(1.0 - x*x));
}