#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//for multithreading
#include <pthread.h>

#define DOUBLEMODE 0
#define DOUBLEERRORMODE 1
#define RATIOMODE 2

#ifndef MODE
#define MODE DOUBLEMODE
#endif

#define NUM_THREADS 6

#if MODE == DOUBLEERRORMODE
#include "double-error.h"
#define VALUETYPE double_error
#define EXPTYPE double_error
#elif MODE == RATIOMODE
#include "ratio.h"
#define VALUETYPE rational
#define EXPTYPE unsigned int
#else
#include "double.h"
#define VALUETYPE double
#define EXPTYPE double
#endif

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

/*given function df[0] on domain [0,M-1], compute derivatives f' until f^{(K)} and store f^{(K)} in df*/
void differentiate(VALUETYPE* f, VALUETYPE* df, int D, int K){
	VALUETYPE df0[D];
	/*Set zeroth derivative to be f.*/
	for(int i=0; i<D; i++){
		df0[i] = f[i];
		df[i] = f[i];
	}
	for(int k=1; k<=K; k++){
		/*Compute kth derivative of f from (k-1)th.*/
		for(int i=0; i<D; i++) df[i] = difference(df0[(i+1)%D], df0[i]);
		for(int i=0; i<D; i++) df0[i] = df[i];
	}

	/*
	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, int D){
	VALUETYPE integralp = int_to_valuetype(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, int D){
	/*Sf[i][j] will be the integral of f on [min(i,j),max(i,j)]*/
	//VALUETYPE Sf[N][N];
	/*Af[i][j] will be the average of f on [min(i,j),max(i,j)]*/
	//VALUETYPE Af[N][N];
	/*Apparently may become too big for stack or something so have to use malloc instead.*/
	VALUETYPE* Sf[N];
	VALUETYPE* Af[N];
	for(int i=0; i<D;i++){
		Sf[i] = malloc(N*sizeof(VALUETYPE));
		Af[i] = malloc(N*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<D; i++){
			Sf[i][(i+n)%D] = sum(Sf[i][(i+n-1)%D], f[(i+n)%D]);
			Af[i][(i+n)%D] = ratio(Sf[i][(i+n)%D],int_to_valuetype(n+1));
		}
	}
	
	/*Compute maximal function by picking the largest average*/
	for(int i=0; i<D; i++){
		Mf[i] = Af[i][i];
		for(int j=1; j<=D; j++){
			Mf[i] = maximum(Af[i][(i+j)%D], Mf[i]);
			Mf[i] = maximum(Af[(i-j+D)%D][i], Mf[i]);
		}
	}

	/*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 D, VALUETYPE* f){
	VALUETYPE Mf[N];
	/*This is the only O(D^2) operation in here so makes a lot of sense to only compute once and avoid repeating it.*/
	compute_maximalfunction(f,Mf,D);

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

	for(int i=0; i<D; i++){
		df[1][i] = f[i];
		df[0][i] = f[i];
		dMf[0][i] = Mf[i];
		dMf[1][i] = Mf[i];
	}

	VALUETYPE tmp[N];

	for(int k=1; k<=16; k++){
		/*Compute kth derivative of f and Mf*/
		differentiate(df[(k+1)%2],df[k%2],D,1);
		differentiate(dMf[(k+1)%2],dMf[k%2],D,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[k%2],p,D);
		VALUETYPE intdMfp = integratep(dMf[k%2],p,D);

		//printf("%d: %f / %f = %f\n",k,intdMfp[k],intdfp[k],intdMfp[k]/intdfp[k]);

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

		VALUETYPE r = ratio(intdMfp, intdfp);
		double t = to_double(r);
		if(
		//(k==1 && t>.9)
		//||
		//(k==2 && t>=.5)
		//||
		(k==3 && t>.53)
		//||
		//(k==4 && t>=.5)
		||
		(k==5 && t>=.58)
		||
		(k==6 && t>=.58)
		||
		(k==7 && t>=.69)
		||
		(k==8 && t>=.83)
		||
		(k==9 && t>=.8699)
		||
		(k==10 && t>=.919)
		||
		(k==11 && t>=.97)
		||
		(k==12 && t>=.97)
		||
		(k==13 && t>=.98)
		||
		(k==14 && t>=.98)
		||
		(k==15 && t>=.9817)
		||
		(k==16 && t>=.9817)
		){
			printf("f: ");
			for(int i=0; i<D; i++) printf("%1.0f ",to_double(f[i]));
			printf("\n");
			char s[128];
			root_to_string(s,r,p);
			printf("ratio for %dth derivative: %s\n",k,s);
		}
	}
}

/*Writes into f the values of the function indexed by i. Returns the size of the support of f, or -1 if there is no function with that index.*/
int generate_function(VALUETYPE* f, int i){
	int s=0;
	int d=1;
	while(d <= N){
		/*number of strings of length d that are not all 0s or all 1s*/
		int powd = (1<<d)-2;
		if(i-s >= powd){
			s += powd;
			d++;
		}
		else {
			int t = i-s+1;
			/*Set f to the values encoded in bit string t which is a value between 1 and powd = (1<<d)-2.*/
			for(int n=0; n<d; n++) f[n] = int_to_valuetype((t >> n) & 1);
			return d;
		}
	}
	return -1;

}

/*Goes through all functions indexed by those numbers which equal the thread number module NUM_THREADS.*/
void* compute_chunk(void* arguments){
	int i = *((int*)arguments);
	EXPTYPE p = int_to_exptype(1);

	VALUETYPE f[N];
	int d = generate_function(f,i);
	while(d >= 0){
		compute(p,d,f);
		i += NUM_THREADS;
		d = generate_function(f,i);
	}	
}

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

	pthread_t threads[NUM_THREADS];
	int args[NUM_THREADS];
	int result_code;

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

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

	return 0;
}