#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

/*Determines if computations are done with floating point type, floating point with error bounds or fractions.
Floating point without error bounds are the default.*/
#ifndef MODE
#define MODE DOUBLEMODE
#endif

/*number of threads for parallel execution.*/
#define NUM_THREADS 7

/*Make sure the appropriate files are included for the type of computation and sets the data types accordingly.*/
#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

typedef unsigned long index_t;

/*maximum length of the circle*/
#define N 36

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

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

/*Given function f on domain [0,D-1], compute derivative of order k and store result 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.*/
/*For p=∞ instead compute the maximum, i.e. ||f||_∞.*/
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 [i,i+j).*/
	VALUETYPE* Sf[N];
	/*Af[i][j] will be the average of f on [i,i+j).*/
	VALUETYPE* Af[N];
	/*Would use C arrays, but apparently they may become too big for stack or something on my machine so have to use malloc instead.*/
	//VALUETYPE Sf[N][2*N+1];
	//VALUETYPE Af[N][2*N+1];

	for(int i=0; i<D; i++){
		Sf[i] = malloc((2*D+1)*sizeof(VALUETYPE));
		Af[i] = malloc((2*D+1)*sizeof(VALUETYPE));
	}

	/*Recursively compute all integrals and averages over intervals of increasing length*/
	/*We have to consider intervals of length up to twice the length of the circle, i.e. 2*D-1.*/
	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++){
		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));
		}
	}
	/*Because it simplifies code later on we also assign the average on twice the circle, starting in 0.*/
	Af[0][2*D] = Af[0][D];

	/*Finding the maximal average actually looks like the most costly computation, in the whole algorithm, being of order D^2.
	Hence we put effort into making it efficient.
	The strategy is to go trough all possible intervals, starting with those of largest length.
	For each [i,i+l) that we encounter we check for each n∈[i,i+l) if the average on [i,i+l) beats the current best average that we have so far computed for intervals containing n.
	The starting point of the current best interval is stored in Pbest[n] and the length in Lbest[n].
	*/
	int Pbest[N];
	int Lbest[N];
	/*Start with the largest possible interval, [0,2*D), which contains each point n∈[0,D).*/
	for(int i=0; i<D; i++){
		Pbest[i]=0;
		Lbest[i]=2*D;
	}
	/*Go through all lenghts l <= 2*D-1 in decreasing order.*/
	for(int l=2*D-1; l>0; l--) {
		/*Go through all intervals [i,i+l).*/
		for(int i=0; i<D; i++){
			VALUETYPE A = Af[i][l];
			/*Go through all points n ∈ [i,i+l).
			If l>D it suffices to go until i+D-1.*/
			int intervalend = fmin(i+l,i+D);
			int n=i;
			while(n < intervalend){
				/*If the average of the current best interval for n is surely larger than A[i][l] then  the same will be true for all subsequent points in that interval.
				Thus, to save time just skip ahead to the end of that interval.*/
				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{
					/*Since the new average might win set the new best average to be the maximum of the last maximum and the new.*/
					A = maximum(A,Af[Pbest[n%D]][Lbest[n%D]]);
					/*To avoid duplicating computations, look ahead how many subsequent points have as their latest best interval the same interval, and accordingly set their new current best average to be the same.*/
					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;
				}
			}
			/*Since we may be working with errors, the above action of taking maxima may have increased the error.
			That means we may have to update the average of the current interval to carry that increased error.*/
			Af[i][l] = A;
		}
	}
	
	/*Finally store best averages in maximal function.*/
	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, indexed by i, but skip trivial ones, and pick only one representative from each class of translaties.
The return value is the length of the domain, with special values 0 which indicates we skip this function, and -1 which indicates that the index i is too large to represent a function.*/
int generate_each_charf(VALUETYPE* f, index_t i){
	/*In order to find which function corresponds to index i, we efficiently go through all indeces up to i starting at 0.*/
	index_t s=0;
	/*The smallest domain we consider is the circle of length 2.*/
	int d=2;
	/*We skip the functions that are all 0 or all 1s.
	That means every function has a point n such that f[n]=1 and f[n+1]=0.
	Since we also don't want to consider different translates, we may translate each function such that this value n equals d-1.
	That means for each d we have to consider only arrays of length d that begin with 1 and end with 0, of which there exist 2^{d-2} many.*/
	index_t powd = 1UL << d-2;
	/*At this point s is the first index that corresponds to functions on the circle of length d.
	Check if our index i goes past the last index that corresponds to a function on the circle of length d.*/
	while(i-s >= powd){
		/*Increment s,d,powd to correspond to functions of length d+1.*/
		s += powd;
		d++;
		powd = powd << 1;
	}

	/*If our index i corresponds to a function on a circle of length larger than N then return -1 which will causo the program to stop.*/
	if(d>N) return -1;

	/*t encodes a characteristic function on the circle of length d.
	Since we only consider functions with f[0]=0 and f[d-1]=1, set its places accordingly.*/
	index_t t = ( (i-s) << 1 ) + 1UL;

	/*This variable indicates if we want to consider this t.
	Next we will perform test to determine if it should be set to false.*/
	bool is_representative = true;

	/*Check if t encodes a function that is made up of multiple copies of a shorter function.
	We skip those because they and their maximal function are periodic and thus should be considered on a smaller circle instead.*/
	/*Go through all possible lenghts r of that shorter function.*/
	for(int r = 1; r < d; r++){
		/*r has to be a divisor of d.*/
		if(d % r == 0){
			bool is_r_copy = true;
			/*bit mask of length r*/
			index_t ones = 0;
			for(int o = 0; o<r; o++) ones += 1UL <<o;
			/*extract from t an array of length r*/
			index_t tail = t & ones;
			/*check if t is d/r many copies of that array*/
			for(int n = 1; n < d/r; n++) if( ( (t>>n*r) & ones ) != tail ) is_r_copy = false;
			if(is_r_copy) is_representative = false;
		}
	}

	/*We aim to filter out equivalent translates.
	The representative we pick is the one which represents the smallest integer.
	Every function which is a nontrivial translate of itself must in fact consist of multiplie copies of a shorter array (right?).
	Those we already filtered out in the previous step, so every nontrivial translate of that representative indeed represents a strictly smaller integer.
	*/
	if(is_representative){
		index_t ones = 0;
		for(int o = 0; o<d; o++) ones += 1UL << o;
		for(int n=1; n<d; n++) if( t >= ( (t<<n) & ones ) + ((t>>(d-n)) & ones) ) is_representative = false;
	}

	/*Set f to the values encoded in bit array t.*/
	if(is_representative){
		for(int n=0; n<d; n++) f[n] = int_to_valuetype((t >> n) & 1UL);
		return d;
	}
	else return 0;
}

int generate_triangle(VALUETYPE* f, index_t 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, index_t 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, index_t 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
/*Write into pointer s the function corresponding to an index, and that the ratio of the L^p norms of the kth derivative of the maximal function and the function equals r.*/
void format_result(char* s, index_t 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,"& ");
	}
}

/*Given an index compute the ratio of the L^p norms of derivatives up to order k of the maximal function and the function for a given range of exponents p.
If for any order of derivative or exponent p the latest record for that ratio is broken, update the value of the record and the index of the witnessing function.*/
int compute(index_t index, EXPTYPE exponents[P], VALUETYPE (*records_ratio)[K+1][P], index_t (*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(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;
}

/*Save formatted string of current records into pointer text. Each entry is formatted using format_result.*/
void format_results(char* text, EXPTYPE exponents[P], VALUETYPE ratios[K+1][P], index_t 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], index_t 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], index_t 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);
}

/*information given to each thread*/
typedef struct { EXPTYPE exponents[P]; VALUETYPE (*records_ratio)[K+1][P]; index_t (*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 mod 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];
	index_t 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;

	/*Start threads that do the computation.*/
	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);
	}

	/*Capture user input.*/
	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); //Stop once all threads have finished.

	/*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;
}