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

/* The work in f(i) is proportional to i. */
double f(long i) {
    long first = i * (i + 1) / 2;
    long last = first + i;
    double result = 0.0;

    for (long j = first; j <= last; j++) {
        result += sin(j);
    }
    return result;
}

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

    double *finished = calloc(thread_count, sizeof(double));
    double sum = 0.0;

    double start = omp_get_wtime();
#   pragma omp parallel num_threads(thread_count) \
        default(none) shared(n, finished, start) reduction(+: sum)
    {
#       pragma omp for schedule(runtime) nowait
        for (long i = 0; i <= n; i++) {
            sum += f(i);
        }
        finished[omp_get_thread_num()] = omp_get_wtime() - start;
    }
    double elapsed = omp_get_wtime() - start;

    const char *schedule = getenv("OMP_SCHEDULE");
    printf("OMP_SCHEDULE: %s\n", schedule != NULL ? schedule : "(not set)");
    printf("Sum: %.6f\n", sum);
    for (int t = 0; t < thread_count; t++) {
        printf("Thread %d finished at %.3f seconds\n", t, finished[t]);
    }
    printf("Elapsed time: %.3f seconds\n", elapsed);
    free(finished);
    return 0;
}