diff options
| -rw-r--r-- | double-error.c | 35 | ||||
| -rw-r--r-- | double-error.h | 4 |
2 files changed, 22 insertions, 17 deletions
diff --git a/double-error.c b/double-error.c index 8df0974..63c820f 100644 --- a/double-error.c +++ b/double-error.c @@ -4,13 +4,16 @@ #include <float.h> #include <stdbool.h> +typedef double vtype; +typedef double etype; #define EPS 2*DBL_EPSILON +#define MIN 2*DBL_MIN -typedef struct {double v; double e;} double_error; +typedef struct {vtype v; etype e;} double_error; double_error int_to_valuetype(int i){ double_error de; - de.v = (double) i; + de.v = (vtype) i; if(abs(i) < 9999999) de.e = 0.0; else de.e = EPS*abs(i); return de; @@ -30,7 +33,7 @@ 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 + ( fabs(de1.v) + fabs(de2.v) + de1.e + de2.e )*EPS + DBL_MIN; + de.e = de1.e + de2.e + ( fabs(de1.v) + fabs(de2.v) + de1.e + de2.e )*EPS + MIN; return de; } @@ -42,11 +45,11 @@ double_error difference(double_error de1, double_error de2){ } bool is_greater_certainly(double_error de1, double_error de2){ - return de1.v-de1.e > de2.v+de2.e + EPS*(fabs(de1.v)+fabs(de1.e)+fabs(de2.v)+fabs(de2.e)) + DBL_MIN; + return de1.v-de1.e > de2.v+de2.e + EPS*(fabs(de1.v)+fabs(de1.e)+fabs(de2.v)+fabs(de2.e)) + MIN; } bool is_greater_possibly(double_error de1, double_error de2){ - return de1.v+de1.e + EPS*(fabs(de1.v)+fabs(de1.e)+fabs(de2.v)+fabs(de2.e)) + DBL_MIN > de2.v-de2.e; + return de1.v+de1.e + EPS*(fabs(de1.v)+fabs(de1.e)+fabs(de2.v)+fabs(de2.e)) + MIN > de2.v-de2.e; } double_error maximum(double_error de1, double_error de2){ @@ -65,11 +68,11 @@ 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; - 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; + vtype p12 = fabs(de1.v)*de2.e*(1+EPS) + MIN; + vtype p21 = fabs(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 + ( fabs(de.v) + s )*EPS + MIN; return de; } @@ -77,7 +80,7 @@ double_error ratio(double_error de1, double_error de2){ double_error de2i; 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 + else de2i.e = 2*de2.e/(de2.v*de2.v) + fabs(de2i.v)*EPS + MIN; //2* is extra margin in case de2.e is near de2.v/2 return product(de1,de2i); } @@ -92,11 +95,11 @@ double_error absolute(double_error de){ double_error power(double_error de, double_error p) { double_error dep; 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; + 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); diff --git a/double-error.h b/double-error.h index dbfe283..3f87ab2 100644 --- a/double-error.h +++ b/double-error.h @@ -3,7 +3,9 @@ #include <stdbool.h> -typedef struct {double v; double e;} double_error; +typedef double vtype; +typedef double etype; +typedef struct {vtype v; etype e;} double_error; double_error int_to_valuetype(int); |
