#include <stdlib.h>
#include <stdio.h>
#include <math.h>
//for multithreading
#include <pthread.h>
#include <stdbool.h>
#include <string.h>
//getLine
#include "misc.h"

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

#ifndef MODE
#define MODE DOUBLEMODE
#endif

#define NUM_THREADS 7

#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

#define STRING_SIZE 65536

/*maximum length of the support of f*/
#define N 32

/*maximal order of derivative*/
#define K 32

/*maximal number of exponents*/
#define P 5

/*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 l=1; l<=k; l++){
		/*Compute lth derivative of f from (l-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 ",valuetype_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++){
		if(exptype_is_infinite(p)) integralp = maximum(integralp,absolute(f[i]));
		else integralp = sum(integralp,power(absolute(f[i]),p));
	}
	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((2*D+1)*sizeof(VALUETYPE));
		Af[i] = malloc((2*D+1)*sizeof(VALUETYPE));
	}

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

	for(int l=2; l < 2*D; l++){
		/*Recursively compute all integrals and averages over intervals of increasing length*/
		for(int i=0; i<D; i++){
			Sf[i][l] = sum(Sf[i][l-1], f[(i+l-1)%D]);
			Af[i][l] = ratio(Sf[i][l],int_to_valuetype(l));
		}
	}

	Af[0][2*D] = Af[0][D];

	int Pbest[N];
	int Lbest[N];
	for(int i=0; i<D; i++){
		Pbest[i]=0;
		Lbest[i]=2*D;
	}

	for(int l=2*D-1; l>0; l--) {
		for(int i=0; i<D; i++){
			VALUETYPE A = Af[i][l];
			int n=i;
			while(n<i+l && n<i+D){
				if(is_greater_certainly(Af[Pbest[n%D]][Lbest[n%D]],Af[i][l])) n += (Pbest[n%D]-n-D)%D + Lbest[n%D];
				else{
					A = maximum(A,Af[Pbest[n%D]][Lbest[n%D]]);
					int k = n+1;
					while((Pbest[k%D] == Pbest[n%D]) && (Lbest[k%D] == Lbest[n%D]) && k<i+l && k<i+D){
						Pbest[k%D] = i;
						Lbest[k%D] = l;
						k++;
					}
					Pbest[n%D] = i;
					Lbest[n%D] = l;
					n=k;
				}
			}
			Af[i][l] = A;
		}
	}
	
	/*Compute maximal function by picking the largest average*/
	for(int i=0; i<D; i++) Mf[i] = Af[Pbest[i]][Lbest[i]];

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

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

/*Generates all characteristic functions of length up N.*/
int generate_each_charf(VALUETYPE* f, int i){
	int s=0;
	int d=0;
	while(d+2 <= N){
		/*number of strings of length d that are not all 0s or all 1s*/
		int powd = 1<<d;
		if(i-s >= powd){
			s += powd;
			d++;
		}
		else {
			int t = i-s;
			/*Set f to the values encoded in bit string t which is a value between 1 and powd = (1<<d)-2.*/
			f[0] = int_to_valuetype(1);
			f[1] = int_to_valuetype(0);
			for(int n=0; n<d; n++) f[2+n] = int_to_valuetype((t >> n) & 1);
			return d+2;
		}
	}
	return -1;
}

int generate_triangle(VALUETYPE* f, int i){
	int j = i+1;
	if(j <= N/2){
		for(int n=0; n < N; n++) f[n] = int_to_valuetype(0);
		for(int n=0; n < j; n++) f[n] = int_to_valuetype(j-n);
		for(int n=1; n < j; n++) f[N-n] = int_to_valuetype(j-n);
		return N;
	}
	else return -1;
}

int generate_random(VALUETYPE* f, int i){
	f[0] = int_to_valuetype(1);
	f[1] = int_to_valuetype(0);
	for(int n=2; n < N; n++) f[n] = int_to_valuetype(rand() % 2);
	return N;
}

/*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){
	return generate_each_charf(f,i);
	//return generate_random(f,i);
	//return generate_triangle(f,i);
}

#define FORMAT_TEXT 0
#define FORMAT_LATEX 1
void format_result(char* s, int index, int k, EXPTYPE p, VALUETYPE r, int format){
	VALUETYPE f[N];
	int d = generate_function(f,index);
	if(format == FORMAT_TEXT){
		strcpy(s,"f: ");
		int l = 3;
		for(int i=0; i<d; i++){
			char v[8];
			valuetype_to_string(v,f[i]);
			l += sprintf(s+l,"%s ",v);
		}
		for(int i=d; i<N; i++) l += sprintf(s+l,"  ");
		char e[16];
		char rts[128];
		exptype_to_string(e,p);
		root_to_string(rts,r,p);
		l += sprintf(s+l,"|f^(%2d)|_%s: %s",k,e,rts);
		/*Compute and print also Mf.*/
		/*Compute and print also df and dMf.*/
		/*
		VALUETYPE Mf[N];
		compute_maximalfunction(f, Mf, d);
		l += sprintf(s+l,"\nM: ");
		for(int i=0; i<d; i++){
			char v[8];
			valuetype_to_string(v,Mf[i]);
			l += sprintf(s+l,"%s ",v);
		}
		VALUETYPE df[2][N];
		VALUETYPE dMf[2][N];
		for(int i=0; i<=1; i++){
			for(int n=0; n<d; n++){
				df[i][n] = f[n];
				dMf[i][n] = Mf[n];
			}
		}
		for(int m=0; m+1<=k; m++){
			differentiate(df[(m+1)%2],df[m%2],d,1);
			differentiate(dMf[(m+1)%2],dMf[m%2],d,1);
			l += sprintf(s+l,"\nf%d: ",m+1);
			for(int i=0; i<d; i++){
				char v[8];
				valuetype_to_string(v,df[(m+1)%2][i]);
				l += sprintf(s+l,"%s ",v);
			}
			l += sprintf(s+l,"\nM%d: ",m+1);
			for(int i=0; i<d; i++){
				char v[8];
				valuetype_to_string(v,dMf[(m+1)%2][i]);
				l += sprintf(s+l,"%s ",v);
			}
		}
		*/
	}
	else if(format == FORMAT_LATEX){
		int l = 0;
		char v[128];
		root_to_latex(v,r,p);
		l += sprintf(s+l,"$%d$ & $%s$",k,v);
		for(int i=0; i<d; i++){
			valuetype_to_latex(v,f[i]);
			l += sprintf(s+l,"& $%s$ ",v);
		}
		for(int i=d; i<N; i++) l += sprintf(s+l,"& ");
	}
}

bool over_threshold_charf(double t, int k){
	return (k==1 && t>=1) || (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) || (k==17 && t>=.987) || (k==18 && t>=.991) || (k==19 && t>=.994) || (k==20 && t>=1.001) || (k==21 && t>=1.009) || (k==22 && t>=1.009) || (k==23 && t>=1.003) || (k==24 && t>=1.174) ;
}

int compute(int index, EXPTYPE exponents[P], VALUETYPE (*records_ratio)[K+1][P], int (*records_index)[K+1][P]){
	VALUETYPE f[N];
	int D = generate_function(f,index);
	/*Immediately abort if index is out of bounds.*/
	if(D <= 0) return D;

	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<=1; i++){
		for(int n=0; n<D; n++){
			df[i][n] = f[n];
			dMf[i][n] = Mf[n];
		}

	}

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

			if(
			//over_threshold_charf(valuetype_to_double(r),k) &&
			is_greater_possibly(r,(*records_ratio)[k][p])){
				/*extra check for printing only because in error mode for some reason floats randomly seem to increase by tiny amounts*/
				if(is_greater_certainly(r,(*records_ratio)[k][p])){
					(*records_index)[k][p] = index;
					char s[1024];
					format_result(s, index, k, exponents[p], r, FORMAT_TEXT);
					printf("%s\n",s);
				}
				(*records_ratio)[k][p] = maximum((*records_ratio)[k][p],r);
			}
		}
	}
	return D;
}

void format_results(char* text, EXPTYPE exponents[P], VALUETYPE ratios[K+1][P], int indeces[K+1][P], int format){
	char beginning[1024];
	char end[1024];
	char beginningp[P][1024];
	char endp[P][1024];
	char newline[4];
	char e[8];
	
	if(format == FORMAT_TEXT){
		strcpy(beginning,"Current records:\n");
		strcpy(end,"");

		for(int p=0; p<P; p++){
			exptype_to_string(e,exponents[p]);
			sprintf(beginningp[p],"exponent %s:",e);
			strcpy(endp[p],"");
		}

		strcpy(newline,"");
	}
	else if(format == FORMAT_LATEX){
		strcpy(beginning,"");
		strcpy(end,"");

		char h[1024];
		strcpy(h,"{cc");
		for(int n=0; n<N; n++) strcat(h,"c");
		strcat(h,"}");
		char f[1024];
		strcpy(f,"$k$ & $r$ & $f$ ");
		for(int n=1; n<N; n++) strcat(f,"& ");
		strcat(f,"\\\\\n\\hline");
		for(int p=0; p<P; p++){
			exptype_to_latex(e,exponents[p]);
			sprintf(beginningp[p],"\\begin{table}\n\\begin{tabular}%s\n%s",h,f);
			sprintf(endp[p],"\\end{tabular}\n\\caption{$p=%s$}\n\\end{table}\n",e);
		}
		strcpy(newline,"\\\\");
	}
	
	char s[STRING_SIZE];
	int l = sprintf(text, "%s\n", beginning);
	for(int p=0; p<P; p++){
		l += sprintf(text+l, "%s\n", beginningp[p]);
		for(int k=0; k<=K; k++){
			format_result(s, indeces[k][p], k, exponents[p], ratios[k][p], format);
			l += sprintf(text+l,"%s%s\n",s,newline);
		}
		if(format == FORMAT_LATEX){
			/*Replace last \\\n by \n */
			l -= 4;
			l += sprintf(text+l,"\n");
		}
		l += sprintf(text+l, "%s\n", endp[p]);
	}
	l += sprintf(text+l, "%s\n", end);
}

void print_records(EXPTYPE exponents[P], VALUETYPE ratios[K+1][P], int indeces[K+1][P]){
	char s[STRING_SIZE];
	format_results(s, exponents, ratios, indeces, FORMAT_TEXT);
	printf("%s\n",s);
}

void print_latex_records_to_file(EXPTYPE exponents[P], VALUETYPE ratios[K+1][P], int indeces[K+1][P], char* filename){
	char s[STRING_SIZE];
	format_results(s, exponents, ratios, indeces, FORMAT_LATEX);

	FILE *fptr;
	fptr = fopen(filename, "w");
	fprintf(fptr, "%s\n",s);
	fclose(fptr);
}

typedef struct { EXPTYPE exponents[P]; VALUETYPE (*records_ratio)[K+1][P]; int (*records_index)[K+1][P]; int num_thread; int* domain_current; int* cont; } Args;

/*Goes through all functions indexed by those numbers which equal the thread number module NUM_THREADS.*/
void* compute_chunk(void* arguments){
	Args* args = arguments;
	int i = args -> num_thread;
	int* domain_current = args -> domain_current;
	int d = 0;
	while(d >= 0){
		d = compute(i, args -> exponents, args -> records_ratio, args -> records_index);
		if(d > *domain_current){
			*domain_current = d;
			printf("Start considering length: %d\n",d);
		}
		i += NUM_THREADS;
	}
	(*(args -> cont))++;
	if(*(args -> cont) >= NUM_THREADS) printf("Calculation finished. Press any button to stop.\n");
	return NULL;
}

int main() {
	EXPTYPE exponents[P];
	exponents[0] = int_to_exptype(1);
	exponents[1] = int_to_exptype(2);
	exponents[2] = int_to_exptype(4);
	exponents[3] = int_to_exptype(8);
	exponents[4] = infinity_to_exptype();

	pthread_t threads[NUM_THREADS];
	Args args[NUM_THREADS];
	VALUETYPE records_ratio[K+1][P];
	int records_index[K+1][P];
	for(int k=0; k<=K; k++) for(int p=0; p<P; p++) records_ratio[k][p] = int_to_valuetype(0);
	int domain_current = 0;
	int cont = 0;

	for (int i = 0; i < NUM_THREADS; i++) {
		printf("In main: Creating thread %d.\n", i);
		memcpy(args[i].exponents,exponents,sizeof exponents);
		args[i].records_ratio = &records_ratio;
		args[i].records_index = &records_index;
		args[i].num_thread = i;
		args[i].domain_current = &domain_current;
		args[i].cont = &cont;
		pthread_create(&threads[i], NULL, compute_chunk, args+i);
	}

	size_t sz = 1024;
	char prmpt[] = "To exit enter q, quit or exit.\nTo print current records to terminal enter r.\nTo print current records to a file enter p.\nTo do so with tex formatting enter t.\n";
	char buff[sz];
	do {
		getLine(prmpt, buff, sz);
		if( 0 == strcmp(buff,"q") || 0 == strcmp(buff,"quit") || 0 == strcmp(buff,"exit") ){
			print_records(exponents,records_ratio,records_index);
			exit(EXIT_SUCCESS);
		}
		else if( 0 == strcmp(buff,"r")) print_records(exponents,records_ratio,records_index);
		else if( 0 == strcmp(buff,"t")){
			getLine("Enter file name:\n", buff, sz);
			print_latex_records_to_file(exponents,records_ratio,records_index,buff);
		}
	} while(cont < NUM_THREADS);

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

	print_records(exponents,records_ratio,records_index);

	return 0;
}