// Parallel calculation of integral

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <windows.h>

const int MAX_THREADS = 16;

// Thread body function
DWORD WINAPI computeIntegral(void* arg);
double simpson(double (*f)(double), double a, double b, int n);

// Argument block
class ArgBlock {
public:
    double (*f)(double);    // Input
    double a;
    double b;
    int n;
    double integral;        // Output
};

double f(double);
double antider_f(double);

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

double antider_f(double x) {
    return (-cos(x));
}

int main() {
    HANDLE hThread[MAX_THREADS];
    DWORD threadID[MAX_THREADS];
    ArgBlock threadArg[MAX_THREADS];
    
    double a, b; int n;
    int numThreads;
    while (true) {
        printf("Number of threads: ");
        if (scanf("%d", &numThreads) < 1 || numThreads <= 0)
            break;
        if (numThreads > MAX_THREADS)
            numThreads = MAX_THREADS;
        
        printf("a, b, n:\n");
        if (scanf("%lf%lf%d", &a, &b, &n) < 3 || n <= 0)
            break;
        
        DWORD startTime = 
            GetTickCount();  // Time in milliseconds since the start of system
        
        int k = n/numThreads;
        if (k < 1)
            k = 1;
        double dx = (b - a)/numThreads;
        for (int i = 0; i < numThreads; ++i) {
            // Define arguments of i-th thread
            threadArg[i].a = a + i*dx;
            threadArg[i].b = threadArg[i].a + dx;
            threadArg[i].n = k;
            threadArg[i].f = &f;
            
            // Create i-th thread
            hThread[i] = CreateThread(
                NULL,           // lpThreadAttributes
                0,              // dwStackSize - default
                &computeIntegral,  // thread starting function
                (void *) &(threadArg[i]), // parameter of thread starting function
                0,              // dwCreationFlags
                &(threadID[i])
            );
        }   

        WaitForMultipleObjects(
            (DWORD) numThreads,
            hThread,
            TRUE,
            INFINITE
        );
            
        double s = 0.;
        for (int i = 0; i < numThreads; ++i) {
            s += threadArg[i].integral;
        }

        DWORD endTime = 
            GetTickCount();  // Time in milliseconds since the start of system

        DWORD dt = endTime - startTime;            
            
        double s0 = antider_f(b) - antider_f(a);
        printf("integral = %.12f\n", s);
        printf("must be:   %.12f\n", s0);
        printf("Computation time: %d ms\n", dt);
    
        for (int i = 0; i < numThreads; ++i) {
            CloseHandle(hThread[i]);
        }
    }        
            
    return 0;
}

DWORD WINAPI computeIntegral(void* threadArgs) {
    ArgBlock* args = (ArgBlock*) threadArgs;
    
    double s = simpson(args->f, args->a, args->b, args->n);
    args->integral = s;
    return 0;
}

double simpson(double (*f)(double), double a, double b, int n) {
    assert(n > 0);
    double h = (b - a)/n;
    double s1 = (*f)(a) + (*f)(b);
    double s2 = 0.;
    double x = a + h;
    for (int i = 0; i < n - 1; ++i) {
        s2 += (*f)(x);
        x += h;
    }
    double s4 = 0.;
    x = a + h/2.;    
    for (int i = 0; i < n; ++i) {
        s4 += (*f)(x);
        x += h;
    }
    return (s1 + 2.*s2 + 4.*s4)*h/6.;
}