-
Notifications
You must be signed in to change notification settings - Fork 1
/
unum2.hpp
200 lines (166 loc) · 6.13 KB
/
unum2.hpp
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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
#pragma once
#include <cstdint>
#include <iostream>
#include <ratio>
// T = support type (enough for 1<<N)
// exact = exact values ordered increasingly
template <class T, int... exacts>
class Punum
{
public:
enum { N = sizeof...(exacts), N2 = N<<3, MASK = N2-1 };
static constexpr int values[] = { exacts...};
using index_t=T;
T v; // index in the N2 space
static constexpr Punum exactindex(int index) { return Punum((N2>>2)+(index<<1)); }; // index*2 + one being the first and zero-based
constexpr Punum next() const { return Punum((v+1) & MASK ); }
constexpr Punum prev() const { return Punum((v-1) & MASK ); }
constexpr Punum abs() const { return isstrictlynegative() ? -(*this): *this; } // could be >= infinity because infinity is sign symmetric
constexpr Punum neg() const { return Punum(-v); };
constexpr Punum inv() const { return Punum(-(v+(N2 >> 1))); }
constexpr bool isinf() const { return v == (N2>>1); }
constexpr bool iszero() const { return v == 0; }
constexpr bool isone() const { return v == (N2>>2); }
constexpr bool isexact() const { return (v&1) == 0; }
//C++14 constexpr bool isfractional() const { Punum w = abs(); return (w.v > 0) && (w.v < (N2>>2)); }
constexpr bool isfractional() const { return v > 0 && (abs().v < (N2>>2)); } // (0 < x < 1) or (-1 < x < 0) == (-1,1) removing 0
constexpr bool isstrictlynegative() const { return v > (N2>>1); } // -inf < x < 0
static constexpr Punum zero() { return Punum(0); }
static constexpr Punum one() { return Punum(N2>>2); }
static constexpr Punum inf() { return Punum(N2>>1); }
// custom operators
constexpr Punum operator-() const { return neg(); }
constexpr Punum operator~() const { return inv(); }
constexpr Punum operator-(const Punum & other) const { return (*this) + (-other); } // uses +
constexpr Punum operator/(const Punum & other) const { return (*this) * other.inv(); } // uses *
constexpr Punum operator+(const Punum & other) const { return *this; } // TBD
constexpr Punum operator*(const Punum & other) const { return *this; } // TBD
// slowproduct
// slowsum
// exactvalue
// iostream
// sqrt
// exp
// conversion
// max
static Punum convert(float x);
static Punum convert(int x);
template<std::intmax_t N,std::intmax_t D>
static constexpr Punum convert(std::ratio<N,D> r);
private:
explicit Punum(int x) : v(x) {} // MAYBE private
};
#if 0
template <class T, int... exacts>
template<std::intmax_t N,std::intmax_t D>
constexpr Punum<T,exacts...> convert(std::ratio<N,D> x)
{
if(x.num == 0)
return zero();
else if(x.num < 0)
return -convert(td::ratio<-N,D>());
int lo = 0;
int hi = N;
if (x.num < x.den)
{
while(true)
{
auto mid = lo + ((hi - lo) >> 1);
if (mid == lo || mid == hi)
break;
// inv(etable[mid]) > x
if(values[mid]*x.den > x.num)
{
lo = mid;
}
else
{
hi = mid;
}
}
if (lo >= 0 && x.num == x.den*values[lo])
return inv(fromexactsindex(lo)); // FIX
else if (hi < N && x.num == x.den*values[hi]) // FIX as in the loop
return inv(fromexactsindex(hi)); // FIX
else if (lo == 0):
return one().prev();
else if hi >= N:
return zero().next();
else
return inv((fromexactsindex(lo).next())); // FIX
}
else
{
while(true)
{
auto mid = lo + ((hi - lo) >> 1);
if (mid == lo || mid == hi)
break;
if(values[mid]*x.den < x.num)
{
lo = mid;
}
else
{
hi = mid;
}
}
if (lo >= 0 && (etable[lo]*x.den == x.num))
return fromexactsindex(lo);
else if (hi < N && (etable[hi]x.den == x.num))
return (fromexactsindex(hi));
else if (lo == 0):
return one().next();
else if hi >= N:
return inf().prev();
else
return fromexactsindex(lo).next();
}
}
#endif
// functionals if needed
template <class T, int... exacts>
constexpr Punum<T,exacts...> neg(Punum<T,exacts...> x) { return -x; }
template <class T, int... exacts>
constexpr Punum<T,exacts...> inv(Punum<T,exacts...> x) { return ~x; }
// Bound object [first..last] or everything or empty
template <class APunum>
class Pbound
{
APunum first,last;
bool empty_ = true;
Pbound() {}
Pbound(APunum afirst, APunum alast): first(afirst),last(alast),empty_(false) {}
Pbound(APunum ax) : first(ax) ,last(ax),empty_(false) {}
constexpr bool isempty() const { return empty_; }
constexpr bool iseverything() const { return !empty_ && first == last.next(); }
constexpr bool issingle() const { return !empty_ && first == last; }
constexpr bool isexact() const { return issingle() && first.isexact(); }
constexpr bool iszero() const { return !empty_ && first.iszero(); }
constexpr bool isone() const { return !empty_ && first.isone(); }
constexpr bool isinf() const { return !empty_ && first.isinf(); }
static constexpr Pbound zero() { return Pbound(APunum::zero()); }
static constexpr Pbound one() { return Pbound(APunum::one()); }
static constexpr Pbound inf() { return Pbound(APunum::inf()); }
static constexpr Pbound everything() { return Pbound(APunum::zero(),APunum::zero().prev()); }
static constexpr Pbound empty() { return Pbound(); }
constexpr Pbound inv() const { return isempty() ? *this : Pbound(~last,~first); }
constexpr Pbound neg() const { return isempty() ? *this : Pbound(-last,-first); }
constexpr Pbound complement() const { return empty_ ? everything() : iseverything() ? empty() : Pbound(last.next(),first.prev()); }
constexpr Pbound operator-(const Pbound & other) const { return (*this) + (-other); } // uses +
constexpr Pbound operator/(const Pbound & other) const { return (*this) * other.inv(); } // uses *
constexpr Pbound operator+(const Pbound & other) const { return *this; } // TBD
constexpr Pbound operator*(const Pbound & other) const { return *this; } // TBD
// in
// intersect
// shortestcover
// outer when operation is CONVEX
// finiteplus
// +
// *
// ==
// ^
// exp
// sqrt
// maximum
};