diff options
| -rw-r--r-- | double-error.c | 79 | ||||
| -rw-r--r-- | double-error.h | 2 |
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); |