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

#define EXACT false
#define VALUETYPE double
#define EXPTYPE double

#if EXACT
#define VALUETYPE ratio
#define EXPTYPE int
#endif


/*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** df, int D, int K){
	int 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(i=0; i<D-k; i++){
			df[k][i] = df[k-1][i+1] - df[k-1][i];
		}
	}
}

/*given function f on domain [0,D-1] compute pth root of integral of |f|^p*/
VALUETYPE integrate(VALUETYPE* f, EXPTYPE p, int D){
	VALUETYPE sum = 0.0;
	for(int i=0;i<D;i++){
		sum += pow(fabs(f[i]),p);
	}
	return pow(sum,1/p);
}

void compute_maximalfunction(VALUETYPE* f, VALUETYPE** Sf, VALUETYPE** Af, VALUETYPE* Mf, int D){
	/*Sf[i][j] will be the integral of f on [min(i,j),max(i,j)]*/
	/*Af[i][j] will be the average of f on [min(i,j),max(i,j)]*/
	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] = Sf[i][i+n-1]+f[i+n];
			Sf[i+n][i] = Sf[i][i+n];
	
			Af[i][i+n] = Sf[i][i+n]/(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(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 ",Mf[i]);
	//printf("\n");
	//for(int i=0;i<D;i++){
		//for(int j=0;j<D;j++) printf("%0.1f ",Af[i][j]);
		//printf("\n");
	//}
}

/*given integer valued function f on domain D, compute ||Mf^{(K)}||_p/||f^{(K)}||_p*/
/*All the other arguments are pointers to storage for intermediate variables. The purpose is that we do not have to allocate new storage with each invocation. Probably there's a more user friendly way but I don't know much about garbage collection.*/
VALUETYPE compute_derivatives(VALUETYPE* f, VALUETYPE** Sf, VALUETYPE** Af, VALUETYPE* Mf, VALUETYPE** df, VALUETYPE** dMf, EXPTYPE p, VALUETYPE* intdf, VALUETYPE* intdMf, int D, int K){

	/*Convert integer valued f to float valued df[0]*/
	df[0] = f;
	
	compute_maximalfunction(f, Sf, Af, Mf, D);
	
	/*Compute 0th until Kth derivative of f and Mf*/
	differentiate(df,D,K);
	differentiate(dMf,D,K);
	
	/*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");
	//}
	
	//for(int k=0;k<=K;k++){
	
	/*Compute L^p norm of derivatives*/
	int k=K;
	intdf[k] = integrate(df[k],p,D);
	intdMf[k] = integrate(dMf[k],p,D);
	//printf("%d: %f / %f = %f\n",k,intdMf[k],intdf[k],intdMf[k]/intdf[k]);
	//}
	
	/*Return ratio of L^p norms*/
	return intdMf[k]/intdf[k];
}

int main() {

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

/*length of the domain*/
int D=5*N;

/*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.*/
int K=3;

/*exponent p of the L^p norm to consider*/
EXPTYPE p = 1.0;

/*allocate memory for f*/
VALUETYPE* f = malloc(D*sizeof(VALUETYPE));

/*allocate memory for temporary variables, such as averages*/
VALUETYPE** Sf = malloc(D*sizeof(VALUETYPE*));
VALUETYPE** Af = malloc(D*sizeof(VALUETYPE*));
VALUETYPE* Mf = malloc(D*sizeof(VALUETYPE));
for(int i=0;i<D;i++){
	Sf[i] = malloc(D*sizeof(VALUETYPE));
	Af[i] = malloc(D*sizeof(VALUETYPE));
}

/*Allocate memory for derivatives*/
VALUETYPE** df = malloc(D*sizeof(VALUETYPE*));
VALUETYPE** dMf = malloc(D*sizeof(VALUETYPE*));
for(int k=0;k<=K;k++){
	df[k] = malloc(D*sizeof(VALUETYPE));
	dMf[k] = malloc(D*sizeof(VALUETYPE));
	for(int i=0;i<D;i++){
		df[k][i] = 0;
		dMf[k][i] = 0;
	}
}
dMf[0] = Mf;

/*Allocate memory for integrals*/
VALUETYPE* intdf = malloc(K*sizeof(VALUETYPE));
VALUETYPE* intdMf = malloc(K*sizeof(VALUETYPE));

/*Allocate memory for ||Mf^{(k)}||_p/||f^{(k)}||_p.*/
VALUETYPE r = 0.0;

/*Iterate over all strings of 0s and 1s with length N. Those will represent f.*/
for(int t=1; t<=(1 << N)-1; t++){
	/*Initiate f to be zero everywhere.*/
	for(int i=0;i<D;i++) f[i]=0.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/2 so that the most interesting part of f^{(K)} and Mf^{(K)} will be around the center of the domain*/
		f[2*N+i+K/2] = (VALUETYPE) ((t >> i) & 1);
		//if(i%3==0) f[2*N+i+K/2] = 1;
	}
	/*Compute ||Mf^{(k)}||_p/||f^{(k)}||_p.*/
	r = compute_derivatives(f, Sf, Af, Mf, df, dMf, p, intdf, intdMf, D, K);
	//printf("%.3d: %.3f  \n",t,r);
	/*Print f and ||Mf^{(k)}||_p/||f^{(k)}||_p if the latter is close to 1/2.*/
	if(r>.4997){
		printf("f: ");
		for(int i=0;i<D;i++) printf("%1.0f ",f[i]);
		printf("\n");
		printf("%.4f\n",r);
	}
}

return 0;
}