// This program is under development
#include "stdio.h"
#include "math.h"
#include "assert.h"
#include "omp.h"

#define EPS 0.001
#define NUM_THREADS 4

typedef double f_t(double x);
typedef struct {
	double so,se;
	double xo,xe;
} calc_t;

#define T get_thread_num()

double sim(f_t * f, double a, double b, int n) {
	double par = f(a)+f(b);
	calc_t s[NUM_THREADS];  // s means "state"
	double h=(b-a)/n;
	double in = 2*NUM_THREADS*h;
	assert (n % 8 == 0);
	int k = n >> 1, i=0;

	omp_set_num_threads(NUM_THREADS);

	int j;
	for (j=0;j<NUM_THREADS;j++) {
		s[j].se=s[j].so=0.0;
	}

	#pragma omp parallel
        {
		s[T].xe=a+h*(T+1);
		s[T].xo=s[T].xe+h;

		for (i=0;i<k;i+=2*NUM_THREADS) {
			s[T].se+=f(s[T].xe);
			s[T].so+=f(s[T].xo);
			s[T].xe=s[T].xe+in;
			s[T].xo=s[T].xe+h;
		}
	}
	double so=0.0,se=0.0;
  	for (j=0;j<NUM_THREADS;j++) {
		so+=s[j].so;
		se+=s[j].se;
	}
	return (b-a)/(3*n)*(par+4*so+2*se);
}

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

void test_sim_spec() {
	double r = sim(sin_, 0, M_PI, 8*80);
	printf("test_sim_spec..");
	assert (fabs(r)<EPS);
	printf("DONE\n");
}
void test_sim_non_null() {
	double r = sim(sin_, 0, M_PI/2, 8*80);
	printf("test_sim_non_null..");
	printf(" r=%f \n", r);
	assert (fabs(r)>1.0);
	printf("DONE\n");
}


int main() {
	test_sim_spec();
	test_sim_non_null();
}