Define safe sums and products and finish exact version.

1 files changed, 53 insertions, 47 deletions

diff --git a/charf.c b/charf.c
index ccbc265..97f9e34 100644
--- a/charf.c
+++ b/charf.c
@@ -2,11 +2,13 @@
 #include <stdlib.h>
 #include <math.h>
 
+#ifndef EXACT
 #define EXACT false
+#endif
 
 #if EXACT
 #include "ratio.h"
-#define VALUETYPE ratio
+#define VALUETYPE rational
 #define EXPTYPE unsigned int
 #else
 #include "double.h"
@@ -22,18 +24,18 @@ void differentiate(VALUETYPE** df, int D, int 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];
+			df[k][i] = difference(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 integralp = 0.0;
+VALUETYPE integratep(VALUETYPE* f, EXPTYPE p, int D){
+	VALUETYPE integralp = convert_int(0);
 	for(int i=0;i<D;i++){
-		integralp += pow(fabs(f[i]),p);
+		integralp = sum(integralp,power(absolute(f[i]),p));
 	}
-	return pow(integralp,1/p);
+	return integralp;
 }
 
 void compute_maximalfunction(VALUETYPE* f, VALUETYPE** Sf, VALUETYPE** Af, VALUETYPE* Mf, int D){
@@ -47,10 +49,10 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE** Sf, VALUETYPE** Af, VALUE
 	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][i+n] = sum(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][i+n] = ratio(Sf[i][i+n],convert_int(n+1));
 			Af[i+n][i] = Af[i][i+n];
 		}
 	}
@@ -59,23 +61,25 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE** Sf, VALUETYPE** Af, VALUE
 	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];
+			if(is_greater(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");
-	//}
+	/*
+	printf("Mf    ");
+	for(int i=0;i<D;i++) printf("%+0.1f ",to_double(Mf[i]));
+	printf("\n");
+	for(int i=0;i<D;i++){
+		for(int j=0;j<D;j++) printf("%0.1f ",to_double(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){
+VALUETYPE compute_derivatives(VALUETYPE* f, VALUETYPE** Sf, VALUETYPE** Af, VALUETYPE* Mf, VALUETYPE** df, VALUETYPE** dMf, EXPTYPE p, VALUETYPE* intdfp, VALUETYPE* intdMfp, int D, int K){
 
 	/*Convert integer valued f to float valued df[0]*/
 	df[0] = f;
@@ -87,41 +91,43 @@ VALUETYPE compute_derivatives(VALUETYPE* f, VALUETYPE** Sf, VALUETYPE** Af, VALU
 	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++){
+		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]);
+	intdfp[k] = integratep(df[k],p,D);
+	intdMfp[k] = integratep(dMf[k],p,D);
+	//printf("%d: %f / %f = %f\n",k,intdMfp[k],intdfp[k],intdMfp[k]/intdfp[k]);
 	//}
 	
 	/*Return ratio of L^p norms*/
-	return intdMf[k]/intdf[k];
+	return ratio(intdMfp[k], intdfp[k]);
 }
 
 int main() {
 
 /*length of the support of f*/
-int N=16;
-
-/*length of the domain*/
-int D=5*N;
+int N=14;
 
 /*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;
+int K=14;
+
+/*length of the domain*/
+int D=N+K+1;
 
 /*exponent p of the L^p norm to consider*/
-EXPTYPE p = 1.0;
+EXPTYPE p = 1;
 
 /*allocate memory for f*/
 VALUETYPE* f = malloc(D*sizeof(VALUETYPE));
@@ -142,41 +148,41 @@ 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;
+		df[k][i] = convert_int(0);
+		dMf[k][i] = convert_int(0);
 	}
 }
 dMf[0] = Mf;
 
 /*Allocate memory for integrals*/
-VALUETYPE* intdf = malloc(K*sizeof(VALUETYPE));
-VALUETYPE* intdMf = malloc(K*sizeof(VALUETYPE));
+VALUETYPE* intdfp = malloc(K*sizeof(VALUETYPE));
+VALUETYPE* intdMfp = malloc(K*sizeof(VALUETYPE));
 
 /*Allocate memory for ||Mf^{(k)}||_p/||f^{(k)}||_p.*/
-VALUETYPE r = 0.0;
+VALUETYPE r = convert_int(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<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/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);
+		f[(D-N)/2+i+K/2] = convert_int((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);
+	r = compute_derivatives(f, Sf, Af, Mf, df, dMf, p, intdfp, intdMfp, 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){
+	if(to_double(r)>.4997){
 		printf("f: ");
-		for(int i=0;i<D;i++) printf("%1.0f ",f[i]);
+		for(int i=0;i<D;i++) printf("%1.0f ",to_double(f[i]));
 		printf("\n");
-		printf("%.4f\n",r);
+		printf("%.4f\n",to_double(r));
 	}
 }
-
+	
 return 0;
 }