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

#include <sys/time.h>
#include <unistd.h>

#define START 0
#define STOP 1
#define MIN_VALUE 0
#define MAX_VALUE 1

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

double my_random(double start, double stop);
int getmilis();

double monte_carlo(int steps);
double rectangles(int steps);
double trapezoids(int steps);

#define TEST(func) do { \
        t1 = getmilis(); \
        i = iterations; \
        while(i--) result = func(n); \
        t2 = getmilis(); \
        printf("integral (" #func "): %lf, time: %d ms\n", result, t2 - t1); \
    while(0)

int main(int argc, char** argv) {
    int n, i;
    int t1, t2;
    int iterations;
    double result;

    if(argc < 3) {
        fprintf(stderr, "usage: %s [number of shots] [iterations]\n", argv[0]);

        return EXIT_FAILURE;
    }

    srand(time(NULL));

    n = atoi(argv[1]);
    iterations = atoi(argv[2]);

    TEST(monte_carlo);
    TEST(rectangles);
    TEST(trapezoids);

    return EXIT_SUCCESS;
}

double trapezoids(int steps) {
    double a;
    double area = 0;
    double size = (double)(STOP - START) / steps;

    a = 0;
    while(a + size < STOP) {
        area += (size * (func(a) + func(a + size))) / 2;

        a += size;
    }

    return area;
}

double rectangles(int steps) {
    double a;
    double area = 0;
    double size = (double)(STOP - START) / steps;

    a = 0;
    while(a + size < STOP) {
        area += size * func((a + a + size) / 2);

        a += size;
    }

    return area;
}

double monte_carlo(int steps) {
    int count;
    double values = 0;
    double c = 0, t, y, v;

    /* using kahan summation to minimize errors. */
    for(count = 0; count < steps; count++) {
        double x = my_random(START, STOP);
        v = func(x);
        y = v - c;
        t = values + y;
        c = (t - values) - y;
        values = t;
    }

    return (((double) (STOP - START) / steps) * values);
}

double my_random(double start, double stop) {
    return (((double)rand() / RAND_MAX) * (stop - start)) + start;
}

int getmilis() {
    struct timeval t;

    gettimeofday(&t, NULL);
    
    return ((t.tv_sec * 1000) + (t.tv_usec / 1000)) + 0.5;
}