1
2 /*++
3 Copyright (c) 2015 Microsoft Corporation
4
5 --*/
6 #include "util/rlimit.h"
7 #include "math/hilbert/hilbert_basis.h"
8
9 /*
10 Test generation of linear congruences a la Karr.
11
12 */
13 namespace karr {
14
15 struct matrix {
16 vector<vector<rational> > A;
17 vector<rational> b;
18
sizekarr::matrix19 unsigned size() const { return A.size(); }
20
resetkarr::matrix21 void reset() {
22 A.reset();
23 b.reset();
24 }
25
operator =karr::matrix26 matrix& operator=(matrix const& other) {
27 reset();
28 append(other);
29 return *this;
30 }
31
appendkarr::matrix32 void append(matrix const& other) {
33 A.append(other.A);
34 b.append(other.b);
35 }
36
displaykarr::matrix37 void display(std::ostream& out) {
38 for (unsigned i = 0; i < A.size(); ++i) {
39 for (unsigned j = 0; j < A[i].size(); ++j) {
40 out << A[i][j] << " ";
41 }
42 out << " = " << -b[i] << "\n";
43 }
44 }
45 };
46
47 // treat src as a homogeneous matrix.
dualizeH(matrix & dst,matrix const & src)48 void dualizeH(matrix& dst, matrix const& src) {
49 reslimit rl;
50 hilbert_basis hb(rl);
51 for (unsigned i = 0; i < src.size(); ++i) {
52 vector<rational> v(src.A[i]);
53 v.push_back(src.b[i]);
54 hb.add_eq(v, rational(0));
55 }
56 for (unsigned i = 0; i < 1 + src.A[0].size(); ++i) {
57 hb.set_is_int(i);
58 }
59 lbool is_sat = hb.saturate();
60 hb.display(std::cout);
61 VERIFY(is_sat == l_true);
62 dst.reset();
63 unsigned basis_size = hb.get_basis_size();
64 for (unsigned i = 0; i < basis_size; ++i) {
65 bool is_initial;
66 vector<rational> soln;
67 hb.get_basis_solution(i, soln, is_initial);
68 if (!is_initial) {
69 dst.b.push_back(soln.back());
70 soln.pop_back();
71 dst.A.push_back(soln);
72 }
73 }
74 }
75
76 // treat src as an inhomegeneous matrix.
dualizeI(matrix & dst,matrix const & src)77 void dualizeI(matrix& dst, matrix const& src) {
78 reslimit rl;
79 hilbert_basis hb(rl);
80 for (unsigned i = 0; i < src.size(); ++i) {
81 hb.add_eq(src.A[i], -src.b[i]);
82 }
83 for (unsigned i = 0; i < src.A[0].size(); ++i) {
84 hb.set_is_int(i);
85 }
86 lbool is_sat = hb.saturate();
87 hb.display(std::cout);
88 VERIFY(is_sat == l_true);
89 dst.reset();
90 unsigned basis_size = hb.get_basis_size();
91 bool first_initial = true;
92 for (unsigned i = 0; i < basis_size; ++i) {
93 bool is_initial;
94 vector<rational> soln;
95 hb.get_basis_solution(i, soln, is_initial);
96 if (is_initial && first_initial) {
97 dst.A.push_back(soln);
98 dst.b.push_back(rational(1));
99 first_initial = false;
100 }
101 else if (!is_initial) {
102 dst.A.push_back(soln);
103 dst.b.push_back(rational(0));
104 }
105 }
106 }
107
juxtapose(matrix & dst,matrix const & M,matrix const & N)108 void juxtapose(matrix& dst, matrix const& M, matrix const& N) {
109 dst = M;
110 dst.append(N);
111 }
112
join(matrix & dst,matrix const & M,matrix const & N)113 void join(matrix& dst, matrix const& M, matrix const& N) {
114 matrix MD, ND, dstD;
115 dualizeI(MD, M);
116 dualizeI(ND, N);
117 juxtapose(dstD, MD, ND);
118 dualizeH(dst, dstD);
119 }
120
joinD(matrix & dst,matrix const & MD,matrix const & ND)121 void joinD(matrix& dst, matrix const& MD, matrix const& ND) {
122 matrix dstD;
123 juxtapose(dstD, MD, ND);
124 dualizeH(dst, dstD);
125 }
126
transition(matrix & dst,matrix const & src,matrix const & Ab)127 void transition(
128 matrix& dst,
129 matrix const& src,
130 matrix const& Ab) {
131 matrix T;
132 // length of rows in Ab are twice as long as
133 // length of rows in src.
134 ENSURE(2*src.A[0].size() == Ab.A[0].size());
135 vector<rational> zeros;
136 for (unsigned i = 0; i < src.A[0].size(); ++i) {
137 zeros.push_back(rational(0));
138 }
139 for (unsigned i = 0; i < src.size(); ++i) {
140 T.A.push_back(src.A[i]);
141 T.A.back().append(zeros);
142 }
143 T.b.append(src.b);
144 T.append(Ab);
145
146 T.display(std::cout << "T:\n");
147 matrix TD;
148 dualizeI(TD, T);
149 TD.display(std::cout << "TD:\n");
150 for (unsigned i = 0; i < TD.size(); ++i) {
151 vector<rational> v;
152 v.append(src.size(), TD.A[i].data() + src.size());
153 dst.A.push_back(v);
154 dst.b.push_back(TD.b[i]);
155 }
156 dst.display(std::cout << "dst\n");
157 }
158
V(int i,int j)159 static vector<rational> V(int i, int j) {
160 vector<rational> v;
161 v.push_back(rational(i));
162 v.push_back(rational(j));
163 return v;
164 }
165
V(int i,int j,int k,int l)166 static vector<rational> V(int i, int j, int k, int l) {
167 vector<rational> v;
168 v.push_back(rational(i));
169 v.push_back(rational(j));
170 v.push_back(rational(k));
171 v.push_back(rational(l));
172 return v;
173 }
174
175 #if 0
176 static vector<rational> V(int i, int j, int k, int l, int m) {
177 vector<rational> v;
178 v.push_back(rational(i));
179 v.push_back(rational(j));
180 v.push_back(rational(k));
181 v.push_back(rational(l));
182 v.push_back(rational(m));
183 return v;
184 }
185 #endif
186
V(int i,int j,int k,int l,int x,int y,int z)187 static vector<rational> V(int i, int j, int k, int l, int x, int y, int z) {
188 vector<rational> v;
189 v.push_back(rational(i));
190 v.push_back(rational(j));
191 v.push_back(rational(k));
192 v.push_back(rational(l));
193 v.push_back(rational(x));
194 v.push_back(rational(y));
195 v.push_back(rational(z));
196 return v;
197 }
198
199 #define R(_x_) rational(_x_)
200
201
tst1()202 static void tst1() {
203 matrix Theta;
204 matrix Ab;
205
206 //
207 Theta.A.push_back(V(1, 0));
208 Theta.b.push_back(R(0));
209 Theta.A.push_back(V(0, 1));
210 Theta.b.push_back(R(-2));
211
212 Theta.display(std::cout << "Theta\n");
213
214 Ab.A.push_back(V(-1, 0, 1, 0));
215 Ab.b.push_back(R(1));
216 Ab.A.push_back(V(-1, -2, 0, 1));
217 Ab.b.push_back(R(1));
218
219 Ab.display(std::cout << "Ab\n");
220
221 matrix ThetaD;
222 dualizeI(ThetaD, Theta);
223 ThetaD.display(std::cout);
224
225 matrix t1D, e1;
226 transition(t1D, Theta, Ab);
227 joinD(e1, t1D, ThetaD);
228
229 t1D.display(std::cout << "t1D\n");
230 e1.display(std::cout << "e1\n");
231
232 matrix t2D, e2;
233 transition(t2D, e1, Ab);
234 joinD(e2, t2D, ThetaD);
235
236 t2D.display(std::cout << "t2D\n");
237 e2.display(std::cout << "e2\n");
238 }
239
tst2()240 void tst2() {
241 /**
242 0 0 0 0 0 0 0 = 0
243 0 0 0 0 0 0 0 = 0
244 0 0 0 0 0 0 0 = 0
245 0 0 0 0 0 0 0 = 0
246 0 0 0 0 1 0 0 = 0
247 0 0 0 0 -1 0 0 = 0
248 0 1 0 0 0 0 0 = 0
249 0 -1 0 0 0 0 0 = 0
250 0 0 0 2 0 0 0 = 0
251 0 0 0 -2 0 0 0 = 0
252 */
253
254 matrix ND;
255 ND.A.push_back(V(0,0,0,0,1,0,0)); ND.b.push_back(R(0));
256 ND.A.push_back(V(0,0,0,0,-1,0,0)); ND.b.push_back(R(0));
257 ND.A.push_back(V(0,1,0,0,0,0,0)); ND.b.push_back(R(0));
258 ND.A.push_back(V(0,-1,0,0,0,0,0)); ND.b.push_back(R(0));
259 ND.A.push_back(V(0,0,0,2,0,0,0)); ND.b.push_back(R(0));
260 ND.A.push_back(V(0,0,0,-2,0,0,0)); ND.b.push_back(R(0));
261
262 ND.display(std::cout << "ND\n");
263
264 matrix N;
265 dualizeH(N, ND);
266
267 N.display(std::cout << "N\n");
268
269
270 }
271
tst3()272 void tst3() {
273 /**
274 0 0 0 0 1 0 0 = 0
275 0 0 0 0 -1 0 0 = 0
276 0 1 0 0 0 0 0 = 0
277 0 -1 0 0 0 0 0 = 0
278 0 0 0 2 0 0 0 = 0
279 0 0 0 -2 0 0 0 = 0
280 */
281
282 matrix ND;
283 ND.A.push_back(V(1,0)); ND.b.push_back(R(0));
284 ND.A.push_back(V(0,2)); ND.b.push_back(R(0));
285
286 ND.display(std::cout << "ND\n");
287
288 matrix N;
289 dualizeH(N, ND);
290
291 N.display(std::cout << "N\n");
292
293
294 }
295
296 };
297
tst_karr()298 void tst_karr() {
299 karr::tst3();
300 return;
301 karr::tst1();
302 }
303