#define  A  1.0    // Integral starting position
#define  B  2.0    // Integral ending position
#define  W   .001  // Integral step size

#include <stdio.h>
#include <mpi.h>
#define  round(v) (float)(int)(v / W) * W

double f(double x) {
    return x * x + 5 * x + 3;
}

int main(int argc, char* argv[]) {
    int i;
    int id, p;
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &p);
    MPI_Comm_rank(MPI_COMM_WORLD, &id);

    // Calculates job boundaries and prints them in order
    double a, b, da, db, wpp;
    wpp   = (B - A) / p;
    da    = wpp * id;
    db    = wpp * (id + 1);
    a     = A + round(da);
    b     = A + round(db);
    for (i = 0; i < p; i++) {
        if (i == id) printf("Process #%d: f(%f) to f(%f)\n", id, a, b);
        MPI_Barrier(MPI_COMM_WORLD);
    }

    // Processes job
    double sum = 0, x;
    for (x = a; x < b; x += W)
        sum += f(x) * W;

    // Reduces result and displays
    double grand_sum = 0;
    MPI_Reduce(&sum, &grand_sum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
    if (id == 0) printf("Result of integral: %f\n", grand_sum);

    MPI_Finalize();
}