forked from LLNL/fpzip
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfpe.inl
61 lines (58 loc) · 1.39 KB
/
fpe.inl
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
template <typename T, typename U, uint mbits, uint ebits>
FPE<T, U, mbits, ebits>
FPE<T, U, mbits, ebits>::operator +(FPE<T, U, mbits, ebits> x) const
{
// addition
if (equalsign(x))
return u > x.u ? add(x.u) : x.add(u);
else {
x = -x;
return u > x.u ? sub(x.u) : u < x.u ? -x.sub(u) : FPE(U(0));
}
}
template <typename T, typename U, uint mbits, uint ebits>
FPE<T, U, mbits, ebits>
FPE<T, U, mbits, ebits>::operator -(FPE<T, U, mbits, ebits> x) const
{
// subtraction
if (equalsign(x))
return u > x.u ? sub(x.u) : u < x.u ? -x.sub(u) : FPE(U(0));
else {
x = -x;
return u > x.u ? add(x.u) : x.add(u);
}
}
template <typename T, typename U, uint mbits, uint ebits>
FPE<T, U, mbits, ebits>
FPE<T, U, mbits, ebits>::add(U g) const
{
// effective addition f + g, |f| >= |g|
U f = u;
uint k = (f >> m) - (g >> m);
if (k <= m) {
U c = (e - (f & (e - 1)));
U d = (e + (g & (e - 1))) >> k;
if (c < d)
d = (c + d) >> 1;
f += d;
}
return FPE(f);
}
template <typename T, typename U, uint mbits, uint ebits>
FPE<T, U, mbits, ebits>
FPE<T, U, mbits, ebits>::sub(U g) const
{
// effective subtraction f - g, |f| >= |g|
U f = u;
uint k = (f >> m) - (g >> m);
if (k <= m) {
U c = ((0 + (f & (e - 1))) << 1);
U d = ((e + (g & (e - 1))) << 1) >> k;
while (c < d) {
d += d - c;
c += e << 1;
}
f -= d >> 1;
}
return FPE(f);
}