Make error bounds precise and correct (hopefully).

2 files changed, 46 insertions, 35 deletions

diff --git a/double-error.c b/double-error.c
index a918cff..591a820 100644
--- a/double-error.c
+++ b/double-error.c
@@ -5,47 +5,52 @@
 
 #define EPS 2*DBL_EPSILON 
 
-typedef struct {double d; double e;} double_error;
+typedef struct {double v; double e;} double_error;
 
 double_error convert_int(int i){
 	double_error de;
-	de.d = (double) i;
+	de.v = (double) i;
 	de.e = EPS*abs(i);
 	return de;
 }
 
 double_error to_exptype(double d){
 	double_error de;
-	de.d = d;
-	de.e = max(EPS*fabs(d),DBL_MIN);
+	de.v = d;
+	de.e = EPS*fabs(d) + DBL_MIN;
 	return de;
 }
 
+bool to_string(char* s, double_error de){
+	sprintf(s,"%f +/- %.16f",de.v,de.e);
+	return true;
+}
+
 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 + DBL_MIN;
+	de.v = de1.v + de2.v;
+	de.e = de1.e + de2.e + ( fabs(de1.v) + fabs(de2.v) + de1.e + de2.e )*EPS + DBL_MIN;
 	return de;
 }
 
 double_error difference(double_error de1, double_error de2){
 	double_error de2m;
-	de2m.d = -de2.d;
+	de2m.v = -de2.v;
 	de2m.e = de2.e;
 	return sum(de1,de2m);
 }
 
-bool is_greater(double_error de1, double_error de2){ return (de1.d > de2.d); }
+bool is_greater(double_error de1, double_error de2){ return (de1.v > de2.v); }
 
 double_error maximum(double_error de1, double_error de2){
 	double_error de;
-	double max1 = de1.d+de1.e;
-	double max2 = de2.d+de2.e;
+	double max1 = de1.v + de1.e;
+	double max2 = de2.v + de2.e;
 	if( max1 >= max2 ){
-		de.d = de1.d;
+		de.v = de1.v;
 		de.e = de1.e + max1*EPS + DBL_MIN;
 	} else {
-		de.d = de2.d;
+		de.v = de2.v;
 		de.e = de2.e + max2*EPS + DBL_MIN;
 	}
 
@@ -54,46 +59,52 @@ double_error maximum(double_error de1, double_error de2){
 
 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 + DBL_MIN;
+	de.v = de1.v * de2.v;
+	double p12 = fabs(de1.v)*de2.e*(1+EPS) + DBL_MIN;
+	double p21 = fabs(de2.v)*de1.e*(1+EPS) + DBL_MIN;
+	double p22 = de1.e*de2.e*(1+EPS) + DBL_MIN;
+	double s = (p12 + p21 + p22)*(1+2*EPS);
+	de.e = s + ( fabs(de.v) + s )*EPS + DBL_MIN;
 	return de;
 }
 
 double_error ratio(double_error de1, double_error de2){
 	double_error de2i;
-	de2i.d = 1/de2.d;
-	if(de2.e > de2.d/2) de2i.e = INFINITY;
-	else de2i.e = 2*de2.e/(de2.d*de2.d) + fabs(de2i.d)*EPS + DBL_MIN; //2* is extra margin in case de2.e is near de2.d/2
+	de2i.v = 1/de2.v;
+	if(de2.e > fabs(de2.v)/2) de2i.e = INFINITY;
+	else de2i.e = 2*de2.e/(de2.v*de2.v) + fabs(de2i.v)*EPS + DBL_MIN; //2* is extra margin in case de2.e is near de2.v/2
 	return product(de1,de2i);
 
 }
 
 double_error absolute(double_error de){
 	double_error dea;
-	dea.d = fabs(de.d);
+	dea.v = fabs(de.v);
 	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;
+	dep.v = pow(de.v,p.v);
+	double dmin = (de.v - de.e)*(1-EPS);
+	double dmax = (de.v + de.e)*(1+EPS);
+	double pmin = (p.v - p.e)*(1-EPS);
+	double pmax = (p.v + p.e)*(1+EPS);
+	double rmin,rmax;
+	if(dmin >= 1){
+		rmin = pow(dmin,pmin);
+		rmax = pow(dmax,pmax);
+	} else if(dmax >= 1) {
+		rmin = pow(dmin,pmax);
+		rmax = pow(dmax,pmax);
+	} else {
+		rmin = pow(dmin,pmax);
+		rmax = pow(dmax,pmin);
 	}
+	if(dmin >= 0) dep.e = rmax-rmin + 2*EPS*rmax;
+	else dep.e = (1+2*EPS)*rmax;
 	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;
-}
+double to_double(double_error de){ return de.v; }
diff --git a/double-error.h b/double-error.h
index c4dadb9..bce8cd5 100644
--- a/double-error.h
+++ b/double-error.h
@@ -1,4 +1,4 @@
-typedef struct {double d; double e;} double_error;
+typedef struct {double v; double e;} double_error;
 
 double_error convert_int(int);