Catch float overflow of float operations most of the time.

1 files changed, 30 insertions, 22 deletions

diff --git a/double-error.c b/double-error.c
index 4daac9e..46f73aa 100644
--- a/double-error.c
+++ b/double-error.c
@@ -8,6 +8,7 @@ typedef long double vtype;
 typedef long double etype;
 #define EPS 2*LDBL_EPSILON
 #define MIN 2*LDBL_MIN
+#define MAX LDBL_MAX
 
 typedef struct {vtype v; etype e;} double_error;
 
@@ -38,7 +39,8 @@ bool exptype_is_infinite(double_error d){ return d.v == INFINITY; }
 double_error sum(double_error de1, double_error de2){
 	double_error de;
 	de.v = de1.v + de2.v;
-	de.e = de1.e + de2.e + ( vabsolute(de1.v) + vabsolute(de2.v) + de1.e + de2.e )*EPS + MIN;
+	if(MAX-vabsolute(de1.v) <= vabsolute(de2.v)) de.e = INFINITY; //This does not account for the error. Should be improved.
+	else de.e = de1.e + de2.e + ( vabsolute(de1.v) + vabsolute(de2.v) + de1.e + de2.e )*EPS + MIN;
 	return de;
 }
 
@@ -73,11 +75,14 @@ double_error maximum(double_error de1, double_error de2){
 double_error product(double_error de1, double_error de2){
 	double_error de;
 	de.v = de1.v * de2.v;
-	vtype p12 = vabsolute(de1.v)*de2.e*(1+EPS) + MIN;
-	vtype p21 = vabsolute(de2.v)*de1.e*(1+EPS) + MIN;
-	vtype p22 = de1.e*de2.e*(1+EPS) + MIN;
-	etype s = (p12 + p21 + p22)*(1+2*EPS);
-	de.e = s + ( vabsolute(de.v) + s )*EPS + MIN;
+	if(MAX/vabsolute(de1.v) <= vabsolute(de2.v)) de.e = INFINITY; //This does not account for the error. Should be improved.
+	else{
+		vtype p12 = vabsolute(de1.v)*de2.e*(1+EPS) + MIN;
+		vtype p21 = vabsolute(de2.v)*de1.e*(1+EPS) + MIN;
+		vtype p22 = de1.e*de2.e*(1+EPS) + MIN;
+		etype s = (p12 + p21 + p22)*(1+2*EPS);
+		de.e = s + ( vabsolute(de.v) + s )*EPS + MIN;
+	}
 	return de;
 }
 
@@ -111,23 +116,26 @@ double_error power(double_error de, double_error p) {
 	if(p.e == 0.0 && ((vtype) ((int) p.v) == p.v) && (int) p.v == 1) return de;
 	double_error dep;
 	dep.v = pow(de.v,p.v);
-	vtype dmin = (de.v - de.e)*(1-EPS);
-	vtype dmax = (de.v + de.e)*(1+EPS);
-	vtype pmin = (p.v - p.e)*(1-EPS);
-	vtype pmax = (p.v + p.e)*(1+EPS);
-	vtype 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(logl(MAX) < logl(de.v)*p.v) dep.e = INFINITY; //This should be checked in two step and does not account for the error. Should be improved.
+	else{
+		vtype dmin = (de.v - de.e)*(1-EPS);
+		vtype dmax = (de.v + de.e)*(1+EPS);
+		vtype pmin = (p.v - p.e)*(1-EPS);
+		vtype pmax = (p.v + p.e)*(1+EPS);
+		vtype 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;
 	}
-	if(dmin >= 0) dep.e = rmax-rmin + 2*EPS*rmax;
-	else dep.e = (1+2*EPS)*rmax;
 	return dep;
 }