#include <stdio.h>
#include <mpi.h>

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

int main(void) {
    int my_rank, comm_sz;
    MPI_Init(NULL, NULL);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
    MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

    const double a = 0.0;
    const double b = 1.0;
    const int n = 500000000;
    const double h = (b - a) / n;
    double local_area = 0.0;

    MPI_Barrier(MPI_COMM_WORLD);
    double start = MPI_Wtime();

    for (int i = my_rank; i < n; i += comm_sz) {
        double left = a + i * h;
        double right = a + (i + 1) * h;
        local_area += h * (f(left) + f(right)) / 2.0;
    }

    double total_area = 0.0;
    MPI_Reduce(&local_area, &total_area, 1, MPI_DOUBLE, MPI_SUM,
               0, MPI_COMM_WORLD);

    double elapsed = MPI_Wtime() - start;
    double max_elapsed = 0.0;
    MPI_Reduce(&elapsed, &max_elapsed, 1, MPI_DOUBLE, MPI_MAX,
               0, MPI_COMM_WORLD);

    if (my_rank == 0) {
        printf("Computation and reduction: %.6f seconds\n", max_elapsed);
        printf("Trapezoids: %d\n", n);
        printf("Processes: %d\n", comm_sz);
        printf("Estimated area: %.12f\n", total_area);
        printf("Exact area:     %.12f\n", 1.0 / 3.0);
    }
    MPI_Barrier(MPI_COMM_WORLD); /* PicoCluster shutdown workaround. */
    MPI_Finalize();
    return 0;
}