1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
|
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <float.h>
#define EPS 2*DBL_EPSILON
typedef struct {double v; double e;} double_error;
double_error convert_int(int i){
double_error de;
de.v = (double) i;
de.e = EPS*abs(i);
return de;
}
double_error to_exptype(double d){
double_error de;
de.v = d;
de.e = EPS*fabs(d) + DBL_MIN;
return de;
}
bool to_string(char* s, double_error de){
sprintf(s,"%f +/- %6.1e",de.v,de.e);
return true;
}
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;
return de;
}
double_error difference(double_error de1, double_error de2){
double_error de2m;
de2m.v = -de2.v;
de2m.e = de2.e;
return sum(de1,de2m);
}
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.v + de1.e;
double max2 = de2.v + de2.e;
if( max1 >= max2 ){
de.v = de1.v;
de.e = de1.e + max1*EPS + DBL_MIN;
} else {
de.v = de2.v;
de.e = de2.e + max2*EPS + DBL_MIN;
}
return de;
}
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;
return de;
}
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
return product(de1,de2i);
}
double_error absolute(double_error de){
double_error dea;
dea.v = fabs(de.v);
dea.e = de.e;
return dea;
}
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;
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.v; }
bool root_to_string(char* s, double_error de, double_error p){
double_error de1p = power(de,ratio(convert_int(1),p));
to_string(s,de1p);
return true;
}
|
