Define safe sums and products and finish exact version.

4 files changed, 144 insertions, 72 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;
 }
diff --git a/double.c b/double.c
index 170e56a..9edaffd 100644
--- a/double.c
+++ b/double.c
@@ -1,6 +1,8 @@
 #include <stdlib.h>
 #include <math.h>
 
+double convert_int(int i){ return (double) i; }
+
 double sum(double d1, double d2){ return d1+d2; }
 
 double difference(double d1, double d2){ return d1-d2; }
@@ -11,7 +13,7 @@ double product(double d1, double d2){ return d1*d2; }
 
 double ratio(double d1, double d2){ return d1/d2; }
 
-double absolute(double d){ return abs(d); }
+double absolute(double d){ return fabs(d); }
 
 double power(double d, double p) { return pow(d,p); }
 
diff --git a/ratio.c b/ratio.c
index 8c67538..70a198a 100644
--- a/ratio.c
+++ b/ratio.c
@@ -1,64 +1,125 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include <limits.h>
 
-typedef struct {long int n; long unsigned int d} rational;
+typedef unsigned long long num;
+
+typedef struct {bool s; num n; num d} rational;
+
+unsigned int safe_sum(num n1, num n2){
+	if(n2 > ULLONG_MAX-n1){
+		printf("Sum overflow: Adding %llu and %llu\n",n1,n2);
+		exit(1);
+	}
+	else {
+		return n1+n2;
+	}
+}
+
+num safe_product(num n1, num n2){
+	if(n1 != 0 && n2 > ULLONG_MAX/n1){
+		printf("Product overflow: Multiplying %llu by %llu\n",n1,n2);
+		exit(1);
+	}
+	else {
+		return n1*n2;
+	}
+}
 
 rational convert_int(int i){ 
-	rational r = {i,1};
+	rational r;
+	r.d = 1;
+	r.n = abs(i);
+	if(i >= 0) { r.s=false; } else { r.s=true; }
 	return r;
 }
 
-rational cancel(rational r){
-	long unsigned int a = abs(r.n), b = r.d, c;
+num gcd(num a, num b){
+	num c;
 	while (b) {
 		c = a % b;
 		a = b;
 		b = c;
 	}
-	rational res = {r.n/a,r.d/a};
+	return a;
+}
+
+rational cancel(rational r){
+	num a = gcd(r.n,r.d);
+	rational res = {r.s,r.n/a,r.d/a};
 	return res;
 }
 
 rational sum(rational r1, rational r2){
-	rational r = {r1.n*r2.d+r2.n*r1.d,r1.d*r2.d};
+	num a = gcd(r1.d,r2.d);
+	num r1da = r1.d/a;
+	num r2da = r2.d/a;
+	rational r;
+	r.d = safe_product(r1da,r2.d);
+	num n1 = safe_product(r1.n,r2da);
+	num n2 = safe_product(r2.n,r1da);
+	if(r1.s == r2.s){
+		r.n = safe_sum(n1,n2);
+		r.s = r1.s;
+	}
+	else {
+		if(n1 >= n2) {
+			r.n = n1 - n2;
+			r.s = r1.s;
+		}
+		else {
+			r.n = n2 - n1;
+			r.s = r2.s;
+		}
+	}
 	return cancel(r);
 }
 
 rational difference(rational r1, rational r2){
-	rational r = {r1.n*r2.d-r2.n*r1.d,r1.d*r2.d};
-	return cancel(r);
+	r2.s = !r2.s;
+	return sum(r1,r2);
 }
 
-bool is_greater(rational r1, rational r2){ return (difference(r1,r2).n > 0); }
+bool is_greater(rational r1, rational r2){ return !(difference(r1,r2).s); }
 
 rational product(rational r1, rational r2){
-	rational r = {r1.n*r2.n,r1.d*r2.d};
+	rational r;
+	rational s1 = {r1.s, r1.n, r2.d};
+	rational s2 = {r2.s, r2.n, r1.d};
+	rational t1 = cancel(s1);
+	rational t2 = cancel(s2);
+	r.s = t1.s^t2.s;
+	r.n = safe_product(t1.n,t2.n);
+	r.d = safe_product(t1.d,t2.d);
 	return cancel(r);
 }
 
 rational ratio(rational r1, rational r2){
-	rational r = {0,1};
-	if (r2.n>0){
-		rational r = {0,1};
-	}
-	else{
-		rational r = {-r1.n*r2.d,r1.d*(-r2.n)};
-	}
-	return cancel(r);
+	rational r2i = {r2.s,r2.d,r2.n};
+	return product(r1,r2i);
 }
 
-double to_double(rational r){ return ((double)r.n)/((double)r.d); }
-
 rational absolute(rational r){
-	rational s = {abs(r.n),r.d};
+	rational s = {false,r.n,r.d};
 	return s;
 }
 
 rational power(rational r, unsigned int p){
-	rational s = {1,1};
+	rational s = {0,1,1};
 	for(int i = 1; i<= p; i++){
-		rational s = product(r,s);
+		s = product(r,s);
 	}
 	return s;
 }
+
+double to_double(rational r){
+	double i;
+	if(r.s) { i=-1.0; } else { i=1.0; }
+	return i*((double)r.n)/((double)r.d);
+}
+
+bool to_string(char* s,rational r){
+	sprintf(s,"%llu / %llu",r.n,r.d);
+	return true;
+}
diff --git a/ratio.h b/ratio.h
index 26ed1c2..bd6a1c7 100644
--- a/ratio.h
+++ b/ratio.h
@@ -2,7 +2,7 @@
 #include <stdlib.h>
 #include <math.h>
 
-typedef struct {long int n; long unsigned int d}rational;
+typedef struct {bool s; unsigned long long n; unsigned long long d}rational;
 
 rational convert_int(int);
 
@@ -20,4 +20,7 @@ rational absolute(rational);
 
 rational power(rational,unsigned int);
 
+
 double to_double(rational);
+
+bool to_string(char*,rational r);