Define variables const to allow to use statically allocated arrays, but still have to allocate matrix for Mf computation dynamically because it doesn't fit in stack or something.

1 files changed, 47 insertions, 45 deletions

diff --git a/charf.c b/charf.c
index f52d41e..1508c61 100644
--- a/charf.c
+++ b/charf.c
@@ -8,7 +8,7 @@
 #define EXACT false
 #endif
 
-#define NUM_THREADS 4
+#define NUM_THREADS 6
 
 #if EXACT
 #include "ratio.h"
@@ -20,9 +20,18 @@
 #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, int D, int K){
+void differentiate(VALUETYPE* f, VALUETYPE* df){
 	/*Set zeroth derivative to be f.*/
 	for(int i=0; i<D; i++){
 		df[i] = f[i];
@@ -38,23 +47,33 @@ void differentiate(VALUETYPE* f, VALUETYPE* df, int D, int K){
 			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, int D){
+VALUETYPE integratep(VALUETYPE* f, EXPTYPE p){
 	VALUETYPE integralp = convert_int(0);
 	for(int i=0;i<D;i++){
-		integralp = sum(integralp,power(absolute(f[i]),p));
+		VALUETYPE padd = power(absolute(f[i]),p);
+		integralp = sum(integralp,padd);
 	}
 	return integralp;
 }
 
-void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){
+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 = malloc(D*sizeof(VALUETYPE*));
+	//VALUETYPE Sf[D][D];
 	/*Af[i][j] will be the average of f on [min(i,j),max(i,j)]*/
-	VALUETYPE** Af = malloc(D*sizeof(VALUETYPE*));
-	for(int i=0;i<D;i++){
+	//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));
 	}
@@ -82,7 +101,7 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){
 			if(is_greater(Af[i][j], Mf[i])) Mf[i] = Af[i][j];
 		}
 	}
-	
+
 	/*Print computed functions and averages*/
 	/*
 	printf("Mf    ");
@@ -93,37 +112,37 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){
 		printf("\n");
 	}
 	*/
-	for(int i=0;i<D;i++){
+
+	for(int i=0; i<D; i++){
 		free(Sf[i]);
 		free(Af[i]);
 	}
-	free(Sf);
-	free(Af);
 }
 
-void compute(int N, int K, int D, EXPTYPE p, int t){
+void compute(EXPTYPE p, int t){
 	/*allocate memory for f*/
-	VALUETYPE* f = malloc(D*sizeof(VALUETYPE));
+	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+K))/2+K+i] = convert_int((t >> i) & 1);
+		f[(D-(N+2*K))/2+K+i] = convert_int((t >> i) & 1);
 		//if(i%3==0) f[2*N+i+K/2] = 1;
 	}
 
-	VALUETYPE* Mf = malloc(D*sizeof(VALUETYPE));
-	compute_maximalfunction(f, Mf, D);
+	VALUETYPE Mf[D];
+	compute_maximalfunction(f, Mf);
 
 	/*Allocate memory for derivatives.*/
-	VALUETYPE* df = malloc(D*sizeof(VALUETYPE));
-	VALUETYPE* dMf = malloc(D*sizeof(VALUETYPE));
+	VALUETYPE df[D];
+	VALUETYPE dMf[D];
 
 	/*Compute Kth derivative of f and Mf*/
-	differentiate(f,df,D,K);
-	differentiate(Mf,dMf,D,K);
+	differentiate(f,df);
+	differentiate(Mf,dMf);
+
 	
 	/*Print derivatives*/
 	/*
@@ -138,69 +157,52 @@ void compute(int N, int K, int D, EXPTYPE p, int t){
 	*/
 
 	/*Compute L^p norm of derivatives*/
-	VALUETYPE intdfp = integratep(df,p,D);
-	VALUETYPE intdMfp = integratep(dMf,p,D);
+	VALUETYPE intdfp = integratep(df,p);
+	VALUETYPE intdMfp = integratep(dMf,p);
 	//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);
 
-	free(df);
-	free(dMf);
-
 	//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)>.6)
+	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));
 	}
-	free(f);
-	free(Mf);
 }
 
-struct Args {int* N; int* K; int* D; EXPTYPE* p; int* t; };
+struct Args {EXPTYPE* p; int* t; };
 
 void* perform_work(void* arguments){
 	struct Args *args = arguments;
-	int N = *(args->N);
-	int K = *(args->K);
-	int D = *(args->D);
 	EXPTYPE p = *(args->p);
 	int* s = args->t;
 
 	while(*s <= (1 << N)-1){
 		int t = *s;
 		*s = t+1;
-		compute(N,K,D,p,t);
+		compute(p,t);
 	}
 	return NULL;
 }
 
 int main() {
-	/*length of the support of f*/
-	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.*/
-	int K=10;
-
-	/*length of the domain, minimum N+K to make sure the support of f^{(k)} belongs to the domain*/
-	int D=N+K+N/2;
-
 	/*exponent p of the L^p norm to consider*/
 	EXPTYPE p = 1;
 
-	int t = 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 = {&N, &K, &D, &p, &t};
+	struct Args args = {&p, &t};
 
 	pthread_t threads[NUM_THREADS];
 	int result_code;