Implement double arithmetic with error bounds.

7 files changed, 132 insertions, 6 deletions

diff --git a/charf.c b/charf.c
index 5aabe94..8911ffe 100644
--- a/charf.c
+++ b/charf.c
@@ -4,13 +4,21 @@
 #include <pthread.h>
 #include <assert.h>
 
-#ifndef EXACT
-#define EXACT false
+#define DOUBLEMODE 0
+#define DOUBLEERRORMODE 1
+#define RATIOMODE 2
+
+#ifndef MODE
+#define MODE DOUBLEMODE
 #endif
 
 #define NUM_THREADS 6
 
-#if EXACT
+#if MODE == DOUBLEERRORMODE
+#include "double-error.h"
+#define VALUETYPE double_error
+#define EXPTYPE double_error
+#elif MODE == RATIOMODE
 #include "ratio.h"
 #define VALUETYPE rational
 #define EXPTYPE unsigned int
@@ -173,7 +181,9 @@ void compute(EXPTYPE p, int t){
 		printf("f: ");
 		for(int i=0;i<D;i++) printf("%1.0f ",to_double(f[i]));
 		printf("\n");
-		printf("%.4f\n",to_double(r));
+		char s[128];
+		to_string(s,r);
+		printf("%.4f  (%s)\n",to_double(r),s);
 	}
 }
 
@@ -194,7 +204,7 @@ void* perform_work(void* arguments){
 
 int main() {
 	/*exponent p of the L^p norm to consider*/
-	EXPTYPE p = 1;
+	EXPTYPE p = to_exptype(1);
 
 	/*Iterate over all strings of 0s and 1s with length N. Those will represent f.*/
 	int t = 1;
diff --git a/double-error.c b/double-error.c
new file mode 100644
index 0000000..aaaa405
--- /dev/null
+++ b/double-error.c
@@ -0,0 +1,83 @@
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <float.h>
+
+#define EPS 2*DBL_EPSILON 
+
+typedef struct {double d; double e;} double_error;
+
+double_error convert_int(int i){
+	double_error de;
+	de.d = (double) i;
+	de.e = EPS*abs(i);
+	return de;
+}
+
+double_error to_exptype(double d){
+	double_error de;
+	de.d = d;
+	de.e = EPS*fabs(d);
+	return de;
+}
+
+double_error sum(double_error de1, double_error de2){
+	double_error de;
+	de.d = de1.d + de2.d;
+	de.e = de1.e + de2.e + ( fabs(de1.d) + fabs(de2.d) + de1.e + de2.e )*EPS;
+	return de;
+}
+
+double_error difference(double_error de1, double_error de2){
+	double_error de2m;
+	de2m.d = -de2.d;
+	de2m.e = de2.e;
+	return sum(de1,de2m);
+}
+
+bool is_greater(double_error de1, double_error de2){ return (de1.d > de2.d); }
+
+double_error product(double_error de1, double_error de2){
+	double_error de;
+	de.d = de1.d * de2.d;
+	double p12 = fabs(de1.d)*de2.e;
+	double p21 = fabs(de2.d)*de1.e;
+	double p22 = de1.e+de2.e;
+	double s = p12 + p21 + p22; 
+	de.e = s + ( fabs(de.d) + s )*EPS;
+	return de;
+}
+
+double_error ratio(double_error de1, double_error de2){
+	double_error de2i;
+	de2i.d = 1/de2.d;
+	de2i.e = 2*de2.e/(de2.d*de2.d) + fabs(de2i.d)*EPS; //2* is random extra margin in case error is large
+	return product(de1,de2i);
+
+}
+
+double_error absolute(double_error de){
+	double_error dea;
+	dea.d = fabs(de.d);
+	dea.e = de.e;
+	return dea;
+}
+
+double_error power(double_error de, double_error p) {
+	double_error dep;
+	dep.d = pow(de.d,p.d);
+	if(de.d > 0 && p.d > 0){
+		dep.e = 2*fabs(dep.d)*( de.e/fabs(de.d) + p.e/fabs(p.d) + EPS ); //again random 2*
+	}
+	else {
+		dep.e = de.e + EPS;
+	}
+	return dep;
+}
+
+double to_double(double_error de){ return de.d; }
+
+bool to_string(char* s, double_error de){
+	sprintf(s,"%f +/- %.16f",de.d,de.e);
+	return true;
+}
diff --git a/double-error.h b/double-error.h
new file mode 100644
index 0000000..413d0dc
--- /dev/null
+++ b/double-error.h
@@ -0,0 +1,23 @@
+typedef struct {double d; double e;} double_error;
+
+double_error convert_int(int);
+
+double_error to_exptype(double);
+
+bool is_greater(double_error,double_error);
+
+double_error sum(double_error,double_error);
+
+double_error difference(double_error,double_error);
+
+double_error product(double_error,double_error);
+
+double_error ratio(double_error,double_error);
+
+double_error absolute(double_error);
+
+double_error power(double_error,double_error);
+
+double to_double(double_error);
+
+bool to_string(char*,double_error);
diff --git a/double.c b/double.c
index d31de70..a587a99 100644
--- a/double.c
+++ b/double.c
@@ -4,6 +4,8 @@
 
 double convert_int(int i){ return (double) i; }
 
+double to_exptype(double d){ return d; }
+
 double sum(double d1, double d2){ return d1+d2; }
 
 double difference(double d1, double d2){ return d1-d2; }
diff --git a/double.h b/double.h
index d035f47..56065ac 100644
--- a/double.h
+++ b/double.h
@@ -1,5 +1,7 @@
 double convert_int(int);
 
+double to_exptype(double);
+
 bool is_greater(double,double);
 
 double sum(double,double);
diff --git a/ratio.c b/ratio.c
index 8f86f93..384f944 100644
--- a/ratio.c
+++ b/ratio.c
@@ -2,6 +2,7 @@
 #include <stdlib.h>
 #include <math.h>
 #include <limits.h>
+#include <float.h>
 
 typedef unsigned long long num;
 
@@ -35,6 +36,8 @@ rational convert_int(int i){
 	return r;
 }
 
+unsigned int to_exptype(unsigned int i){ return i; } 
+
 num gcd(num a, num b){
 	num c;
 	while (b) {
@@ -120,6 +123,7 @@ double to_double(rational r){
 }
 
 bool to_string(char* s,rational r){
-	sprintf(s,"%llu / %llu",r.n,r.d);
+	double f = to_double(r);
+	sprintf(s,"%llu / %llu = %f +- %.24f",r.n,r.d,f,f*DBL_EPSILON);
 	return true;
 }
diff --git a/ratio.h b/ratio.h
index cf04c6c..56701f3 100644
--- a/ratio.h
+++ b/ratio.h
@@ -6,6 +6,8 @@ typedef struct {bool s; unsigned long long n; unsigned long long d;} rational;
 
 rational convert_int(int);
 
+unsigned int to_exptype(unsigned int);
+
 bool is_greater(rational,rational);
 
 rational sum(rational,rational);