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

#ifndef EXACT
#define EXACT false
#endif

#define NUM_THREADS 6

#if EXACT
#include "ratio.h"
#define VALUETYPE rational
#define EXPTYPE unsigned int
#else
#include "double.h"
#define VALUETYPE double
#define EXPTYPE double
#endif

/*length of the support of f*/
static const int N=16;

/*order of the derivative to consider. Should not be larger than (D-N)/2 because then the support of f^{(K)} reaches outside of our domain.*/
static const int K=12;

/*length of the domain, minimum N+2*K to make sure the support of f^{(k)} belongs to the domain*/
static const int D=N+2*K;


/*given function df[0] on domain [0,D-1], compute derivatives f' until f^{(K)} and store them in df[1] to df[K]*/
void differentiate(VALUETYPE* f, VALUETYPE* df){
	/*Set zeroth derivative to be f.*/
	for(int i=0; i<D; i++){
		df[i] = f[i];
	}
	for(int k=1; k<=K; k++){
		/*Compute kth derivative of f from (k-1)th.*/
		/*Only compute derivatives at arguments i < D-k, because for larger i we would need data from outside the domain of f.*/
		for(int i=0; i<D-k; i++){
			df[i] = difference(df[i+1], df[i]);
		}
		/*Make them zero for arguments where their dependence reaches outside the domain.*/
		for(int i=D-k; i<D; i++){
			df[i] = convert_int(0);
		}
	}

	/*
	printf("df    ");
	for(int i=0;i<D;i++) printf("%+6.1f ",to_double(df[i]));
	printf("\n");
	*/
}

/*given function f on domain [0,D-1] compute pth root of integral of |f|^p*/
VALUETYPE integratep(VALUETYPE* f, EXPTYPE p){
	VALUETYPE integralp = convert_int(0);
	for(int i=0;i<D;i++){
		VALUETYPE padd = power(absolute(f[i]),p);
		integralp = sum(integralp,padd);
	}
	return integralp;
}

void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf){
	/*Sf[i][j] will be the integral of f on [min(i,j),max(i,j)]*/
	//VALUETYPE Sf[D][D];
	/*Af[i][j] will be the average of f on [min(i,j),max(i,j)]*/
	//VALUETYPE Af[D][D];
	/*Apparently may become too big for stack or something.*/
	VALUETYPE* Sf[D];
	VALUETYPE* Af[D];
	for(int i=0; i<D;i++){
		Sf[i] = malloc(D*sizeof(VALUETYPE));
		Af[i] = malloc(D*sizeof(VALUETYPE));
	}

	for(int i=0;i<D;i++) {
		Sf[i][i] = f[i];
		Af[i][i] = f[i];
	}
	
	for(int n=1;n<D;n++){
		/*Recursively compute all integrals and averages over intervals of increasing length*/
		for(int i=0;i+n<D;i++){
			Sf[i][i+n] = sum(Sf[i][i+n-1], f[i+n]);
			Sf[i+n][i] = Sf[i][i+n];
	
			Af[i][i+n] = ratio(Sf[i][i+n],convert_int(n+1));
			Af[i+n][i] = Af[i][i+n];
		}
	}
	
	/*Compute maximal function by picking the largest average*/
	for(int i=0;i<D;i++){
		Mf[i] = Af[i][i];
		for(int j=0;j<D;j++){
			if(is_greater(Af[i][j], Mf[i])) Mf[i] = Af[i][j];
		}
	}

	/*Print computed functions and averages*/
	/*
	printf("Mf    ");
	for(int i=0;i<D;i++) printf("%+0.1f ",to_double(Mf[i]));
	printf("\n");
	for(int i=0;i<D;i++){
		for(int j=0;j<D;j++) printf("%0.1f ",to_double(Af[i][j]));
		printf("\n");
	}
	*/

	for(int i=0; i<D; i++){
		free(Sf[i]);
		free(Af[i]);
	}
}

void compute(EXPTYPE p, int t){
	/*allocate memory for f*/
	VALUETYPE f[D];
	/*Initiate f to be zero everywhere.*/
	for(int i=0;i<D;i++) f[i] = convert_int(0);
	//for(int i=0;i<N;i++) f[N+i] = rand() %2;
	/*In the middle of the domain set f to the values encoded in bit string t*/
	for(int i=0;i<N;i++){
		/*Since we care about the Kth derivative, which in i depends on f on [i,i+K], shift f to the right by K so that the support of f^{(K)} stays within the domain.*/
		f[(D-(N+2*K))/2+K+i] = convert_int((t >> i) & 1);
		//if(i%3==0) f[2*N+i+K/2] = 1;
	}

	double max = 1.0;
	unsigned int max_operations_per_variable = 1;

	VALUETYPE Mf[D];
	compute_maximalfunction(f, Mf);

	max_operations_per_variable += (D-1)*(1+max_operations_per_variable); //sum of nominators and denominators
	max_operations_per_variable += 2*(1+max_operations_per_variable)+1; //division of nominator and denominator

	/*Allocate memory for derivatives.*/
	VALUETYPE df[D];
	VALUETYPE dMf[D];

	/*Compute Kth derivative of f and Mf*/
	differentiate(f,df);
	differentiate(Mf,dMf);

	max = max*(1<<K); //times 2^K
	for(int i=1; i<= K; i++) max_operations_per_variable = max_operations_per_variable*2+1;

	
	/*Print derivatives*/
	/*
	for(int k=0;k<=K;k++){
		printf("f  %d: ",k);
		for(int i=0;i<D;i++) printf("%+0.1f ",df[k][i]);
		printf("\n");
		printf("Mf %d: ",k);
		for(int i=0;i<D;i++) printf("%+0.1f ",dMf[k][i]);
		printf("\n");
	}
	*/

	/*Compute L^p norm of derivatives*/
	VALUETYPE intdfp = integratep(df,p);
	VALUETYPE intdMfp = integratep(dMf,p);

	max = max*pow(D,to_double(ratio(convert_int(1),convert_int(p))));
	max_operations_per_variable = D*(max_operations_per_variable+1) + D-1 + 1;
	//printf("%d: %f / %f = %f\n",k,intdMfp[k],intdfp[k],intdMfp[k]/intdfp[k]);

	/*Compute ||Mf^{(k)}||_p/||f^{(k)}||_p.*/

	double divisor = to_double(intdfp);

	VALUETYPE r = ratio(intdMfp, intdfp);

	max_operations_per_variable = 2*max_operations_per_variable+1;

	//printf("%.3d: %.3f  \n",t,r);
	/*Print f and ||Mf^{(k)}||_p/||f^{(k)}||_p if the latter is close to 1/2.*/
	//if(to_double(r)>.4997)
	if(to_double(r)>.65)
	{
		printf("f: ");
		for(int i=0;i<D;i++) printf("%1.0f ",to_double(f[i]));
		printf("\n");
		printf("%.4f\n",to_double(r));
		if(!EXACT) printf("max error: %u * %.1f * 2^{-51} / %f = %.12f\n", max_operations_per_variable, max, divisor, max_operations_per_variable*max/divisor*pow(2.0,-51));
	}
}

struct Args {EXPTYPE* p; int* t; };

void* perform_work(void* arguments){
	struct Args *args = arguments;
	EXPTYPE p = *(args->p);
	int* s = args->t;

	while(*s <= (1 << N)-1){
		int t = *s;
		*s = t+1;
		compute(p,t);
	}
	return NULL;
}

int main() {
	/*exponent p of the L^p norm to consider*/
	EXPTYPE p = 1;

	/*Iterate over all strings of 0s and 1s with length N. Those will represent f.*/
	int t = 1;

	//for(int t=1; t<=(1 << N)-1; t++){
		//compute(N,K,D,p,t);
	//}

	struct Args args = {&p, &t};

	pthread_t threads[NUM_THREADS];
	int result_code;

	for (int i = 0; i < NUM_THREADS; i++) {
		printf("In main: Creating thread %d.\n", i);
		result_code = pthread_create(&threads[i], NULL, perform_work, &args);
		assert(!result_code);
	}

	// wait for each thread to complete
	for (int i = 0; i < NUM_THREADS; i++) {
		result_code = pthread_join(threads[i], NULL);
		assert(!result_code);
		printf("In main: Thread %d has ended.\n", i);
	}

	return 0;
}