Using OpenMP in for loop results in incorrect output

Question

I am writing a program to integrate a function using both trapezium and midpoint rules. I wanted to use parallelism to speed up execution. In the code below I assumed since each integration is independent of the next this would be ideal for parallel processing, but when OpenMP is used the output is incorrect, albeit much faster.

I tried using "#pragma omp parallel private" to for the variables used but that resulted in correct output, at single thread speeds. The same thing occurred when using "#pragma omp critical".

Why is this occurring?

Thanks in advance for any advice/help.

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

#define lowerbound -3
#define upperbound -1

// The function to be integrated
long double f(long double x) {
    return ((6 * pow(x,-2) + 5 * pow(x,-7) - 1) / (2 * sin(x) - pow(tan(x),2)) + 2);
}

// Trapezium integration
long double trapezium(unsigned long n) {
    // let h = trapezium widths
    long double h = ((upperbound - lowerbound) / (n * 1.0f));

    // add f0 and fn
    long double output = f(lowerbound) + f(upperbound);

    // add f1 to fn-1
    #pragma omp parallel for
    for (unsigned long i=1; i < n; i++) {
        // h multiplied by i gives the x coordinate to find the f(x) of
        long double val = f(lowerbound + h * i) * 2;
        output += val;
    }

    // multiply by half h as per trapezium rule
    return ((h/2) * output);
}

// Midpoint integration
long double midpoint(unsigned long n) {

    // let h = width
    long double h = ((upperbound - lowerbound) / (n * 1.0f));

    long double output = 0;

    // add middle rectangles
    #pragma omp parallel for
    for (unsigned long i=0; i < n; i++) {
        // adding h2 due to inital offset of centres of rectangles, multiplying by h to get area
        long double val = f(lowerbound + h/2 + h * i) * h;
        output += val;
    }

    return output;
}



int main() {
    omp_set_num_threads(8);

    // Print answers

    for (int i=1; i < 30; i++) {
        printf("%lu", (unsigned long)pow(2,i));
        printf("    ");
        printf("%20.20Lf", trapezium(pow(2,i)));
        printf("    ");
        printf("%20.20Lf\n", midpoint(pow(2,i)));
    }
}