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

#define LIMIT_A 1
#define LIMIT_B 3
#define NUM_N 100000


double f(double x) {
    return x * sin(x);
}


double integrate(double a, double b, int n) {
    double integral;
    int i;
    double h;
    
    h = (b - a) / ((double) n);
    integral = 0;
    for (i=1; i<=n-1; i++) 
	integral += f(a + ((double) i)*h);
    integral *= h;
    integral += (h/2) * (f(a) + f(b));

    return integral;
}

double exactIntegral(double a, double b) {
    return (sin(b) - sin(a) - (b * cos(b) - a * cos(a)));
}

int main(int argc, double* argv[]) {
    double serial_integral;
    double exact_integral;
    int proc_num, proc_rank;
    double local_a, local_b, local_integral;
    int local_n;
    double parallel_integral;
    double buffer;
    MPI_Status status;
    int i;

    MPI_Init( &argc, (char***) &argv );
    MPI_Comm_size( MPI_COMM_WORLD, &proc_num );
    MPI_Comm_rank( MPI_COMM_WORLD, &proc_rank );
    
    if (proc_rank == 0) {
	printf("---------------------------------------------------------------------\n");
	printf("Serial numerical integration:\n\n");
	printf("Integrating f(x) = x * sin(x) from %e to %e\n",(double) LIMIT_A,(double) LIMIT_B);
	serial_integral = integrate(LIMIT_A,LIMIT_B,NUM_N);
	exact_integral = exactIntegral(LIMIT_A,LIMIT_B);
	printf("Numerical integration value using %d trapezoids: %e\n",NUM_N,serial_integral);
	printf("Exact integration value: %e\n",exact_integral);
	printf("Absolute / relative error: %e / %e\n",fabs(serial_integral - exact_integral),
	       fabs(serial_integral - exact_integral) / fabs(exact_integral));
	printf("---------------------------------------------------------------------\n\n");
	printf("---------------------------------------------------------------------\n");
	printf("Parallel numerical integration:\n\n");
    }
    MPI_Barrier(MPI_COMM_WORLD);
    local_a = (double) LIMIT_A + 
	((((double) LIMIT_B - (double) LIMIT_A) / ((double) proc_num)) * ((double) proc_rank));
    local_b = (double) LIMIT_A + 
	((((double) LIMIT_B - (double) LIMIT_A) / ((double) proc_num)) * ((double) (proc_rank + 1)));
    local_n = NUM_N / proc_num;
    local_integral = integrate(local_a,local_b,local_n);
    printf("|-- Process number %d out of %d processes:\n|--- Integration from %e to %e using %d trapezoids\n"
	   "|--- Local intergration value: %e\n",
	   proc_rank,proc_num,local_a,local_b,local_n,local_integral);
    if (proc_rank == 0) {
	parallel_integral = local_integral;
	for (i=1;i<proc_num;i++) {
	    MPI_Recv(&buffer, 1, MPI_DOUBLE, i, 1, MPI_COMM_WORLD, &status);
	    parallel_integral += buffer;
	}
	printf("Parallel numerical integration value: %e\n",parallel_integral);
	printf("Exact integration value: %e\n",exact_integral);
	printf("Absolute / relative error: %e / %e\n",fabs(parallel_integral - exact_integral),
	       fabs(parallel_integral - exact_integral) / fabs(exact_integral));
	printf("---------------------------------------------------------------------\n\n");
    } else {
	MPI_Send(&local_integral,1,MPI_DOUBLE, 0, 1, MPI_COMM_WORLD);
    }


    fflush(stdout);
    fflush(stderr);
    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Finalize();

    return 0;
}