charf.c:

 - Define function to compute maximal function.
 - Define type of function values as c macro.
 - Use same type also for initiating f.

1 files changed, 96 insertions, 85 deletions

diff --git a/charf.c b/charf.c
index a99508e..5ce7321 100644
--- a/charf.c
+++ b/charf.c
@@ -2,10 +2,18 @@
 #include <stdlib.h>
 #include <math.h>
 
+#define EXACT false
+#define VALUETYPE double
+#define EXPTYPE double
 
-/*given function f on domain [0,D-1], compute derivatives f^{(0)} until f^{(K)} and store them in df*/
-void differentiate(double* f, double** df, int D, int K){
-	df[0] = f;
+#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*/
@@ -17,82 +25,85 @@ void differentiate(double* f, double** df, int D, int K){
 }
 
 /*given function f on domain [0,D-1] compute pth root of integral of |f|^p*/
-double integrate(double* f, double p, int D){
-	double sum = 0.0;
+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);
 }
 
-/*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.*/
-double compute_derivatives(int* f, int** Sf, double** Af, double* Mf, double** df, double** dMf, double p, double* intdf, double* intdMf, int D, int K){
-
-/*Convert integer valued f to float valued df[0]*/
-for(int i=0;i<D;i++) df[0][i] = (double) f[i];
-
-/*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] = 1.0*Sf[i][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] = 1.0*Sf[i][i+n]/(n+1);
-		Af[i+n][i] = Af[i][i+n];
+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];
 	}
-}
-
-/*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];
+	
+	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];
+		}
 	}
-}
-
-/*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");
-//}
-
-/*Compute 0th until Kth derivative of f and Mf*/
-differentiate(df[0],df,D,K);
-differentiate(dMf[0],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]);
+	
+	/*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 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]);
-//}
+	//for(int i=0;i<D;i++){
+		//for(int j=0;j<D;j++) printf("%0.1f ",Af[i][j]);
+		//printf("\n");
+	//}
+}
 
-/*Return ratio of L^p norms*/
-return intdMf[k]/intdf[k];
+/*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() {
@@ -107,26 +118,26 @@ int D=5*N;
 int K=3;
 
 /*exponent p of the L^p norm to consider*/
-double p = 1.0;
+EXPTYPE p = 1.0;
 
 /*allocate memory for f*/
-int* f = malloc(D*sizeof(int));
+VALUETYPE* f = malloc(D*sizeof(VALUETYPE));
 
 /*allocate memory for temporary variables, such as averages*/
-int** Sf = malloc(D*sizeof(int*));
-double** Af = malloc(D*sizeof(double*));
-double* Mf = malloc(D*sizeof(double));
+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(int));
-	Af[i] = malloc(D*sizeof(double));
+	Sf[i] = malloc(D*sizeof(VALUETYPE));
+	Af[i] = malloc(D*sizeof(VALUETYPE));
 }
 
 /*Allocate memory for derivatives*/
-double** df = malloc(D*sizeof(double*));
-double** dMf = malloc(D*sizeof(double*));
+VALUETYPE** df = malloc(D*sizeof(VALUETYPE*));
+VALUETYPE** dMf = malloc(D*sizeof(VALUETYPE*));
 for(int k=0;k<=K;k++){
-	df[k] = malloc(D*sizeof(double));
-	dMf[k] = malloc(D*sizeof(double));
+	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;
@@ -135,21 +146,21 @@ for(int k=0;k<=K;k++){
 dMf[0] = Mf;
 
 /*Allocate memory for integrals*/
-double* intdf = malloc(K*sizeof(double));
-double* intdMf = malloc(K*sizeof(double));
+VALUETYPE* intdf = malloc(K*sizeof(VALUETYPE));
+VALUETYPE* intdMf = malloc(K*sizeof(VALUETYPE));
 
 /*Allocate memory for ||Mf^{(k)}||_p/||f^{(k)}||_p.*/
-double r = 0.0;
+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;
+	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] = (t >> i) & 1;
+		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.*/
@@ -158,7 +169,7 @@ for(int t=1; t<=(1 << N)-1; t++){
 	/*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("%.1d ",f[i]);
+		for(int i=0;i<D;i++) printf("%1.0f ",f[i]);
 		printf("\n");
 		printf("%.4f\n",r);
 	}