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
|
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <float.h>
#define EPS 2*DBL_EPSILON
typedef struct {double d; double e;} double_error;
double_error convert_int(int i){
double_error de;
de.d = (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);
return de;
}
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;
return de;
}
double_error difference(double_error de1, double_error de2){
double_error de2m;
de2m.d = -de2.d;
de2m.e = de2.e;
return sum(de1,de2m);
}
bool is_greater(double_error de1, double_error de2){ return (de1.d > de2.d); }
double_error maximum(double_error de1, double_error de2){
double_error de;
double max1 = de1.d+de1.e;
double max2 = de2.d+de2.e;
if( max1 >= max2 ){
de.d = de1.d;
de.e = de1.e + max1*EPS + DBL_MIN;
} else {
de.d = de2.d;
de.e = de2.e + max2*EPS + DBL_MIN;
}
return de;
}
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;
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
return product(de1,de2i);
}
double_error absolute(double_error de){
double_error dea;
dea.d = fabs(de.d);
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;
}
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;
}
|
