#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

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

int main(int argc, char *argv[]) {
    int thread_count = 4;
    if (argc > 1) {
        thread_count = (int)strtol(argv[1], NULL, 10);
    }

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

    double start = omp_get_wtime();
#   pragma omp parallel num_threads(thread_count)
    {
        int my_rank = omp_get_thread_num();
        int team_size = omp_get_num_threads();
        double my_area = 0.0;

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

#       pragma omp critical
        total_area += my_area;
    }
    double elapsed = omp_get_wtime() - start;

    printf("Trapezoids: %d\n", n);
    printf("Threads: %d\n", thread_count);
    printf("Estimated area: %.12f\n", total_area);
    printf("Exact area:     %.12f\n", 1.0 / 3.0);
    printf("Elapsed time: %.3f seconds\n", elapsed);
    return 0;
}