add Polyhedron_Integer_Minimum for computing the integer minimum of a polyhedron
[barvinok.git] / testlib.cc
blob405e02af670b14c68129d751b98a209756bca517
1 #include <assert.h>
2 #include <sstream>
3 #include <barvinok/barvinok.h>
4 #include <barvinok/set.h>
5 #include <barvinok/options.h>
6 #include <barvinok/evalue.h>
7 #include <barvinok/util.h>
8 #include "conversion.h"
9 #include "evalue_read.h"
10 #include "dpoly.h"
11 #include "lattice_point.h"
12 #include "counter.h"
13 #include "bernoulli.h"
14 #include "hilbert.h"
15 #include "ilp.h"
16 #include "matrix_read.h"
18 using std::cout;
19 using std::cerr;
20 using std::endl;
22 template <typename T>
23 void set_from_string(T& v, const char *s)
25 std::istringstream str(s);
26 str >> v;
29 int test_evalue(struct barvinok_options *options)
31 unsigned nvar, nparam;
32 char **all_vars;
33 evalue *poly1, poly2;
35 poly1 = evalue_read_from_str("(1/4 * n^4 + 1/2 * n^3 + 1/4 * n^2)",
36 NULL, &all_vars, &nvar, &nparam,
37 options->MaxRays);
38 Free_ParamNames(all_vars, nvar+nparam);
40 value_init(poly2.d);
41 evalue_copy(&poly2, poly1);
42 evalue_negate(poly1);
43 eadd(&poly2, poly1);
44 reduce_evalue(poly1);
45 assert(EVALUE_IS_ZERO(*poly1));
46 evalue_free(poly1);
47 free_evalue_refs(&poly2);
50 int test_split_periods(struct barvinok_options *options)
52 unsigned nvar, nparam;
53 char **all_vars;
54 evalue *e;
56 e = evalue_read_from_str("U + 2V + 3 >= 0\n- U -2V >= 0\n- U 10 >= 0\n"
57 "U >= 0\n\n({( 1/3 * U + ( 2/3 * V + 0 ))})",
58 NULL, &all_vars, &nvar, &nparam,
59 options->MaxRays);
60 Free_ParamNames(all_vars, nvar+nparam);
62 evalue_split_periods(e, 2, options->MaxRays);
63 assert(value_zero_p(e->d));
64 assert(e->x.p->type == partition);
65 assert(e->x.p->size == 4);
66 assert(value_zero_p(e->x.p->arr[1].d));
67 assert(e->x.p->arr[1].x.p->type == polynomial);
68 assert(value_zero_p(e->x.p->arr[3].d));
69 assert(e->x.p->arr[3].x.p->type == polynomial);
70 evalue_free(e);
73 int test_specialization(struct barvinok_options *options)
75 Value v;
76 value_init(v);
77 mpq_t count;
78 mpq_init(count);
80 value_set_si(v, 5);
81 dpoly n(2, v);
82 assert(value_cmp_si(n.coeff->p[0], 1) == 0);
83 assert(value_cmp_si(n.coeff->p[1], 5) == 0);
84 assert(value_cmp_si(n.coeff->p[2], 10) == 0);
86 value_set_si(v, 1);
87 dpoly d(2, v, 1);
88 value_set_si(v, 2);
89 dpoly d2(2, v, 1);
90 d *= d2;
91 assert(value_cmp_si(d.coeff->p[0], 2) == 0);
92 assert(value_cmp_si(d.coeff->p[1], 1) == 0);
93 assert(value_cmp_si(d.coeff->p[2], 0) == 0);
95 n.div(d, count, 1);
96 mpq_canonicalize(count);
97 assert(value_cmp_si(mpq_numref(count), 31) == 0);
98 assert(value_cmp_si(mpq_denref(count), 8) == 0);
100 value_set_si(v, -2);
101 dpoly n2(2, v);
102 assert(value_cmp_si(n2.coeff->p[0], 1) == 0);
103 assert(value_cmp_si(n2.coeff->p[1], -2) == 0);
104 assert(value_cmp_si(n2.coeff->p[2], 3) == 0);
106 n2.div(d, count, 1);
107 mpq_canonicalize(count);
108 assert(value_cmp_si(mpq_numref(count), 6) == 0);
109 assert(value_cmp_si(mpq_denref(count), 1) == 0);
111 value_clear(v);
112 mpq_clear(count);
115 int test_lattice_points(struct barvinok_options *options)
117 Param_Vertices V;
118 mat_ZZ tmp;
119 set_from_string(tmp, "[[0 0 0 0 4][0 0 0 0 4][-1 0 1 0 4]]");
120 V.Vertex = zz2matrix(tmp);
121 vec_ZZ lambda;
122 set_from_string(lambda, "[3 5 7]");
123 mat_ZZ rays;
124 set_from_string(rays, "[[0 1 0][4 0 1][0 0 -1]]");
125 term_info num;
126 evalue *point[4];
128 unsigned nvar, nparam;
129 char **all_vars;
130 point[0] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
131 "( 7 * {( 1/4 * a + ( 3/4 * c + 3/4 ) ) } + -21/4 ) ) )",
132 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
133 Free_ParamNames(all_vars, nvar+nparam);
134 point[1] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
135 "( 7 * {( 1/4 * a + ( 3/4 * c + 1/2 ) ) } + -1/2 ) ) )",
136 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
137 Free_ParamNames(all_vars, nvar+nparam);
138 point[2] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
139 "( 7 * {( 1/4 * a + ( 3/4 * c + 1/4 ) ) } + 17/4 ) ) )",
140 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
141 Free_ParamNames(all_vars, nvar+nparam);
142 point[3] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
143 "( 7 * {( 1/4 * a + ( 3/4 * c + 0 ) ) } + 9 ) ) )",
144 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
145 Free_ParamNames(all_vars, nvar+nparam);
147 lattice_point(&V, rays, lambda, &num, 4, options);
148 Matrix_Free(V.Vertex);
150 for (int i = 0; i < 4; ++i) {
151 assert(eequal(num.E[i], point[i]));
152 evalue_free(point[i]);
153 free_evalue_refs(num.E[i]);
154 delete num.E[i];
156 delete [] num.E;
159 static int test_icounter(struct barvinok_options *options)
161 icounter cnt(2);
162 vec_QQ n_coeff;
163 mat_ZZ n_power;
164 mat_ZZ d_power;
165 set_from_string(n_coeff, "[-2/1 1/1]");
166 set_from_string(n_power, "[[2 6][3 6]]");
167 d_power.SetDims(0, 2);
168 cnt.reduce(n_coeff, n_power, d_power);
169 assert(value_cmp_si(mpq_numref(cnt.count), -1) == 0);
170 assert(value_cmp_si(mpq_denref(cnt.count), 1) == 0);
173 static Matrix *matrix_read_from_str(const char *s)
175 std::istringstream str(s);
176 return Matrix_Read(str);
179 static int test_infinite_counter(struct barvinok_options *options)
181 Matrix *M = matrix_read_from_str("1 3\n 1 1 0\n");
182 Polyhedron *ctx = Constraints2Polyhedron(M, options->MaxRays);
183 Matrix_Free(M);
185 /* (1 -1/2 x^5 - 1/2 x^7)/(1-x) */
186 infinite_counter *cnt = new infinite_counter(1, 1);
187 cnt->init(ctx);
188 vec_QQ n_coeff;
189 mat_ZZ n_power;
190 mat_ZZ d_power;
191 set_from_string(n_coeff, "[1/1 -1/2 -1/2]");
192 set_from_string(n_power, "[[0][5][7]]");
193 set_from_string(d_power, "[[1]]");
194 cnt->reduce(n_coeff, n_power, d_power);
195 assert(value_cmp_si(mpq_numref(cnt->count[0]), 6) == 0);
196 assert(value_cmp_si(mpq_denref(cnt->count[0]), 1) == 0);
197 assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) == 0);
198 assert(value_cmp_si(mpq_denref(cnt->count[1]), 1) == 0);
199 delete cnt;
200 Polyhedron_Free(ctx);
202 M = matrix_read_from_str("2 4\n 1 1 0 0\n 1 0 1 0\n");
203 ctx = Constraints2Polyhedron(M, options->MaxRays);
204 Matrix_Free(M);
206 /* (1 - xy)/((1-x)(1-xy)) */
207 cnt = new infinite_counter(2, 3);
208 cnt->init(ctx);
209 set_from_string(n_coeff, "[1/1 -1/1]");
210 set_from_string(n_power, "[[0 0][1 1]]");
211 set_from_string(d_power, "[[1 0][1 1]]");
212 cnt->reduce(n_coeff, n_power, d_power);
213 assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) != 0);
214 assert(value_cmp_si(mpq_numref(cnt->count[2]), 0) == 0);
215 assert(value_cmp_si(mpq_denref(cnt->count[2]), 1) == 0);
216 assert(value_cmp_si(mpq_numref(cnt->count[3]), 0) == 0);
217 assert(value_cmp_si(mpq_denref(cnt->count[3]), 1) == 0);
218 delete cnt;
220 cnt = new infinite_counter(2, 2);
221 cnt->init(ctx);
222 set_from_string(n_coeff, "[-1/2 1/1 -1/3]");
223 set_from_string(n_power, "[[2 6][3 6]]");
224 d_power.SetDims(0, 2);
225 cnt->reduce(n_coeff, n_power, d_power);
226 assert(value_cmp_si(mpq_numref(cnt->count[0]), 1) == 0);
227 assert(value_cmp_si(mpq_denref(cnt->count[0]), 6) == 0);
228 assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) == 0);
229 assert(value_cmp_si(mpq_denref(cnt->count[1]), 1) == 0);
230 assert(value_cmp_si(mpq_numref(cnt->count[2]), 0) == 0);
231 assert(value_cmp_si(mpq_denref(cnt->count[2]), 1) == 0);
232 delete cnt;
234 cnt = new infinite_counter(2, 2);
235 cnt->init(ctx);
236 set_from_string(n_coeff, "[1/1]");
237 set_from_string(n_power, "[[0 11]]");
238 set_from_string(d_power, "[[0 1]]");
239 cnt->reduce(n_coeff, n_power, d_power);
240 assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) != 0);
241 assert(value_cmp_si(mpq_numref(cnt->count[2]), 0) == 0);
242 assert(value_cmp_si(mpq_denref(cnt->count[2]), 1) == 0);
243 delete cnt;
245 Polyhedron_Free(ctx);
247 return 0;
250 static int test_series(struct barvinok_options *options)
252 Matrix *M = matrix_read_from_str(
253 "12 11\n"
254 " 0 1 0 0 0 0 0 1 0 0 3 \n"
255 " 0 0 1 0 0 0 0 -1 1 0 -5 \n"
256 " 0 0 0 1 0 0 0 0 -2 -1 6 \n"
257 " 0 0 0 0 1 0 0 1 1 0 5 \n"
258 " 0 0 0 0 0 1 0 0 -1 0 0 \n"
259 " 0 0 0 0 0 0 1 -2 0 -1 -3 \n"
260 " 1 0 0 0 0 0 0 2 0 1 3 \n"
261 " 1 0 0 0 0 0 0 1 -1 0 5 \n"
262 " 1 0 0 0 0 0 0 -1 -1 0 -5 \n"
263 " 1 0 0 0 0 0 0 -1 0 0 -3 \n"
264 " 1 0 0 0 0 0 0 0 2 1 -6 \n"
265 " 1 0 0 0 0 0 0 0 1 0 0 \n");
266 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
267 Matrix_Free(M);
268 Polyhedron *C = Universe_Polyhedron(3);
269 gen_fun *gf = barvinok_series_with_options(P, C, options);
270 Polyhedron_Free(P);
271 Polyhedron_Free(C);
272 delete gf;
274 M = matrix_read_from_str(
275 "7 8\n"
276 " 0 1 1 0 0 1 0 2 \n"
277 " 0 0 0 1 0 -2 0 6 \n"
278 " 0 0 0 0 1 -1 0 -1 \n"
279 " 0 0 0 0 0 0 1 0 \n"
280 " 1 0 1 0 0 0 0 0 \n"
281 " 1 0 -1 0 0 -1 0 -2 \n"
282 " 1 0 0 0 0 1 0 -3 \n");
283 P = Constraints2Polyhedron(M, options->MaxRays);
284 Matrix_Free(M);
285 C = Universe_Polyhedron(2);
286 gf = barvinok_series_with_options(P, C, options);
287 Polyhedron_Free(P);
288 Polyhedron_Free(C);
289 delete gf;
291 M = matrix_read_from_str(
292 "2 3\n"
293 "1 1 0\n"
294 "1 -1 10\n");
295 P = Constraints2Polyhedron(M, options->MaxRays);
296 Matrix_Free(M);
297 C = Universe_Polyhedron(1);
298 gf = barvinok_series_with_options(P, C, options);
299 Polyhedron_Free(P);
300 Polyhedron_Free(C);
301 gen_fun *sum = gf->summate(1, options);
302 delete gf;
303 delete sum;
305 return 0;
308 int test_todd(struct barvinok_options *options)
310 tcounter t(2, options->max_index);
311 assert(value_cmp_si(t.todd.coeff->p[0], 1) == 0);
312 assert(value_cmp_si(t.todd.coeff->p[1], -3) == 0);
313 assert(value_cmp_si(t.todd.coeff->p[2], 3) == 0);
314 assert(value_cmp_si(t.todd_denom->p[0], 1) == 0);
315 assert(value_cmp_si(t.todd_denom->p[1], 6) == 0);
316 assert(value_cmp_si(t.todd_denom->p[2], 36) == 0);
318 vec_ZZ lambda;
319 set_from_string(lambda, "[1 -1]");
320 zz2values(lambda, t.lambda->p);
322 mat_ZZ rays;
323 set_from_string(rays, "[[-1 0][-1 1]]");
325 QQ one(1, 1);
327 vec_ZZ v;
328 set_from_string(v, "[2 0 1]");
329 Vector *vertex = Vector_Alloc(3);
330 zz2values(v, vertex->p);
332 t.handle(rays, vertex->p, one, 1, options);
333 assert(value_cmp_si(mpq_numref(t.count), 71) == 0);
334 assert(value_cmp_si(mpq_denref(t.count), 24) == 0);
336 set_from_string(rays, "[[0 -1][1 -1]]");
337 set_from_string(v, "[0 2 1]");
338 zz2values(v, vertex->p);
340 t.handle(rays, vertex->p, one, 1, options);
341 assert(value_cmp_si(mpq_numref(t.count), 71) == 0);
342 assert(value_cmp_si(mpq_denref(t.count), 12) == 0);
344 set_from_string(rays, "[[1 0][0 1]]");
345 set_from_string(v, "[0 0 1]");
346 zz2values(v, vertex->p);
348 t.handle(rays, vertex->p, one, 1, options);
349 assert(value_cmp_si(mpq_numref(t.count), 6) == 0);
350 assert(value_cmp_si(mpq_denref(t.count), 1) == 0);
352 Vector_Free(vertex);
355 int test_bernoulli(struct barvinok_options *options)
357 struct bernoulli_coef *bernoulli_coef;
358 struct poly_list *faulhaber, *bernoulli;
359 bernoulli_coef = bernoulli_coef_compute(2);
360 faulhaber = faulhaber_compute(4);
361 bernoulli_coef = bernoulli_coef_compute(8);
362 assert(value_cmp_si(bernoulli_coef->num->p[6], 1) == 0);
363 assert(value_cmp_si(bernoulli_coef->den->p[6], 42) == 0);
364 assert(value_cmp_si(faulhaber->poly[3]->p[0], 0) == 0);
365 assert(value_cmp_si(faulhaber->poly[3]->p[1], 0) == 0);
366 assert(value_cmp_si(faulhaber->poly[3]->p[2], 1) == 0);
367 assert(value_cmp_si(faulhaber->poly[3]->p[3], -2) == 0);
368 assert(value_cmp_si(faulhaber->poly[3]->p[4], 1) == 0);
369 assert(value_cmp_si(faulhaber->poly[3]->p[5], 4) == 0);
371 bernoulli = bernoulli_compute(6);
372 assert(value_cmp_si(bernoulli->poly[6]->p[0], 1) == 0);
373 assert(value_cmp_si(bernoulli->poly[6]->p[1], 0) == 0);
374 assert(value_cmp_si(bernoulli->poly[6]->p[2], -21) == 0);
375 assert(value_cmp_si(bernoulli->poly[6]->p[3], 0) == 0);
376 assert(value_cmp_si(bernoulli->poly[6]->p[4], 105) == 0);
377 assert(value_cmp_si(bernoulli->poly[6]->p[5], -126) == 0);
378 assert(value_cmp_si(bernoulli->poly[6]->p[6], 42) == 0);
379 assert(value_cmp_si(bernoulli->poly[6]->p[7], 42) == 0);
381 unsigned nvar, nparam;
382 char **all_vars;
383 evalue *base, *sum1, *sum2;
384 base = evalue_read_from_str("(1 * n + 1)", NULL, &all_vars, &nvar, &nparam,
385 options->MaxRays);
387 sum1 = evalue_polynomial(faulhaber->poly[3], base);
388 Free_ParamNames(all_vars, nvar+nparam);
390 sum2 = evalue_read_from_str("(1/4 * n^4 + 1/2 * n^3 + 1/4 * n^2)",
391 NULL, &all_vars, &nvar, &nparam,
392 options->MaxRays);
393 Free_ParamNames(all_vars, nvar+nparam);
394 assert(eequal(sum1, sum2));
395 evalue_free(base);
396 evalue_free(sum1);
397 evalue_free(sum2);
400 int test_bernoulli_sum(struct barvinok_options *options)
402 unsigned nvar, nparam;
403 char **all_vars;
404 evalue *e, *sum1, *sum2;
405 e = evalue_read_from_str("i + -1 >= 0\n -i + n >= 0\n\n 1 + (-1 *i) + i^2",
406 "i", &all_vars, &nvar, &nparam,
407 options->MaxRays);
408 Free_ParamNames(all_vars, nvar+nparam);
410 sum1 = Bernoulli_sum_evalue(e, 1, options);
411 sum2 = evalue_read_from_str("n -1 >= 0\n\n (1/3 * n^3 + 2/3 * n)",
412 NULL, &all_vars, &nvar, &nparam,
413 options->MaxRays);
414 Free_ParamNames(all_vars, nvar+nparam);
415 evalue_negate(sum1);
416 eadd(sum2, sum1);
417 reduce_evalue(sum1);
418 assert(EVALUE_IS_ZERO(*sum1));
419 evalue_free(e);
420 evalue_free(sum1);
422 e = evalue_read_from_str("-i + -1 >= 0\n i + n >= 0\n\n 1 + i + i^2",
423 "i", &all_vars, &nvar, &nparam,
424 options->MaxRays);
425 Free_ParamNames(all_vars, nvar+nparam);
426 sum1 = Bernoulli_sum_evalue(e, 1, options);
427 evalue_negate(sum1);
428 eadd(sum2, sum1);
429 reduce_evalue(sum1);
430 assert(EVALUE_IS_ZERO(*sum1));
431 evalue_free(e);
432 evalue_free(sum1);
434 evalue_free(sum2);
436 e = evalue_read_from_str("i + 4 >= 0\n -i + n >= 0\n\n i",
437 "i", &all_vars, &nvar, &nparam,
438 options->MaxRays);
439 Free_ParamNames(all_vars, nvar+nparam);
440 sum1 = Bernoulli_sum_evalue(e, 1, options);
441 sum2 = evalue_read_from_str("n + 0 >= 0\n\n (1/2 * n^2 + 1/2 * n + (-10))\n"
442 "n + 4 >= 0\n -n -1 >= 0\n\n (1/2 * n^2 + 1/2 * n + (-10))",
443 NULL, &all_vars, &nvar, &nparam,
444 options->MaxRays);
445 Free_ParamNames(all_vars, nvar+nparam);
446 evalue_negate(sum1);
447 eadd(sum2, sum1);
448 reduce_evalue(sum1);
449 assert(EVALUE_IS_ZERO(*sum1));
450 evalue_free(e);
451 evalue_free(sum1);
452 evalue_free(sum2);
454 e = evalue_read_from_str("i -5 >= 0\n -i + n >= 0\n j -1 >= 0\n i -j >= 0\n"
455 "k -1 >= 0\n j -k >= 0\n\n1",
456 "i,j,k", &all_vars, &nvar, &nparam,
457 options->MaxRays);
458 Free_ParamNames(all_vars, nvar+nparam);
459 sum1 = Bernoulli_sum_evalue(e, 3, options);
460 sum2 = evalue_read_from_str("n -5 >= 0\n\n"
461 "1/6 * n^3 + 1/2 * n^2 + 1/3 * n + -20",
462 NULL, &all_vars, &nvar, &nparam,
463 options->MaxRays);
464 Free_ParamNames(all_vars, nvar+nparam);
465 evalue_negate(sum1);
466 eadd(sum2, sum1);
467 reduce_evalue(sum1);
468 assert(EVALUE_IS_ZERO(*sum1));
469 evalue_free(e);
470 evalue_free(sum1);
471 evalue_free(sum2);
474 int test_hilbert(struct barvinok_options *options)
476 Matrix *M = matrix_read_from_str(
477 "2 4\n"
478 " 1 4 -3 0 \n"
479 " 1 3 2 0 \n");
480 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
481 Matrix_Free(M);
483 M = Cone_Hilbert_Basis(P, options->MaxRays);
484 assert(M->NbRows = 5);
485 assert(M->NbColumns = 3);
486 Matrix_Free(M);
488 M = Cone_Integer_Hull(P, options);
489 assert(M->NbRows = 4);
490 assert(M->NbColumns = 3);
491 Matrix_Free(M);
493 Polyhedron_Free(P);
496 int test_ilp(struct barvinok_options *options)
498 Matrix *M = matrix_read_from_str(
499 "2 4\n"
500 " 1 4 -3 0 \n"
501 " 1 3 2 0 \n");
502 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
503 Matrix_Free(M);
504 Vector *obj = Vector_Alloc(2);
505 value_set_si(obj->p[0], 7);
506 value_set_si(obj->p[1], -1);
507 Value min, max;
508 value_init(min);
509 value_init(max);
511 value_set_si(min, 1);
512 value_set_si(max, 17);
513 Vector *opt = Polyhedron_Integer_Minimum(P, obj->p, min, max,
514 NULL, 0, options);
515 assert(opt);
516 assert(value_cmp_si(opt->p[0], 1) == 0);
517 assert(value_cmp_si(opt->p[1], 1) == 0);
518 assert(value_cmp_si(opt->p[2], 1) == 0);
519 Vector_Free(opt);
521 value_clear(min);
522 value_clear(max);
523 Vector_Free(obj);
524 Polyhedron_Free(P);
527 int main(int argc, char **argv)
529 struct barvinok_options *options = barvinok_options_new_with_defaults();
530 test_evalue(options);
531 test_split_periods(options);
532 test_specialization(options);
533 test_lattice_points(options);
534 test_icounter(options);
535 test_infinite_counter(options);
536 test_series(options);
537 test_todd(options);
538 test_bernoulli(options);
539 test_bernoulli_sum(options);
540 test_hilbert(options);
541 test_ilp(options);
542 barvinok_options_free(options);