?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Polyhedron_Sample: rename "dim" variable to "space"
[barvinok.git] / testlib.cc
blob0342d8fa1996e82300caae9d203c8bb4da253238
1 #include <assert.h>
2 #include <stdlib.h>
3 #include <sstream>
4 #include <NTL/vec_ZZ.h>
5 #include <NTL/mat_ZZ.h>
6 #include <barvinok/NTL_QQ.h>
7 #include <barvinok/polylib.h>
8 #include <barvinok/barvinok.h>
9 #include <barvinok/set.h>
10 #include <barvinok/options.h>
11 #include <barvinok/basis_reduction.h>
12 #include <barvinok/evalue.h>
13 #include <barvinok/util.h>
14 #include "conversion.h"
15 #include "evalue_read.h"
16 #include "dpoly.h"
17 #include "lattice_point.h"
18 #include "counter.h"
19 #include "bernoulli.h"
20 #include "hilbert.h"
21 #include "hull.h"
22 #include "ilp.h"
23 #include "laurent.h"
24 #include "matrix_read.h"
25 #include "remove_equalities.h"
26 #include "config.h"
27
28 using std::cout;
29 using std::cerr;
30 using std::endl;
31
32 template <typename T>
33 void set_from_string(T& v, const char *s)
34 {
35 std::istringstream str(s);
36 str >> v;
37 }
38
39 static Matrix *matrix_read_from_str(const char *s)
40 {
41 std::istringstream str(s);
42 return Matrix_Read(str);
43 }
44
45 static int test_equalities(struct barvinok_options *options)
46 {
47 Matrix *M = matrix_read_from_str(
48 "11 11\n"
49 " 0 23 0 0 -10 0 0 0 7 -44 -8 \n"
50 " 0 0 23 0 5 0 0 0 -15 114 27 \n"
51 " 0 0 0 1 0 0 0 0 0 -1 2 \n"
52 " 0 0 0 0 0 1 0 0 -1 8 0 \n"
53 " 0 0 0 0 0 0 1 0 0 -1 2 \n"
54 " 0 0 0 0 0 0 0 1 0 -1 -1 \n"
55 " 1 0 0 0 10 0 0 0 -7 44 8 \n"
56 " 1 0 0 0 -5 0 0 0 15 -114 -27 \n"
57 " 1 0 0 0 1 0 0 0 0 0 0 \n"
58 " 1 0 0 0 0 0 0 0 1 -8 0 \n"
59 " 1 0 0 0 0 0 0 0 0 1 -2 \n");
60 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
61 Matrix_Free(M);
62 Polyhedron *Q = P;
63 remove_all_equalities(&P, NULL, NULL, NULL, 2, options->MaxRays);
64 assert(P->NbEq == 0);
65 if (P != Q)
66 Polyhedron_Free(Q);
67 Polyhedron_Free(P);
68 return 0;
69 }
70
71 int test_evalue_read(struct barvinok_options *options)
72 {
73 unsigned nvar, nparam;
74 const char **all_vars;
75 evalue *e1, *e2;
76
77 e1 = evalue_read_from_str("(1 * aa + 2 * a)",
78 NULL, &all_vars, &nvar, &nparam, options->MaxRays);
79 Free_ParamNames(all_vars, nvar+nparam);
80 e2 = evalue_read_from_str("(3 * aa)",
81 NULL, &all_vars, &nvar, &nparam, options->MaxRays);
82 Free_ParamNames(all_vars, nvar+nparam);
83 assert(!eequal(e1, e2));
84 evalue_free(e1);
85 evalue_free(e2);
86 return 0;
87 }
88
89 static void evalue_check_disjoint(evalue *e)
90 {
91 int i, j;
92
93 if (!e)
94 return;
95 if (EVALUE_IS_ZERO(*e))
96 return;
97 for (i = 0; i < e->x.p->size/2; ++i) {
98 Polyhedron *A = EVALUE_DOMAIN(e->x.p->arr[2*i]);
99 for (j = i+1; j < e->x.p->size/2; ++j) {
100 Polyhedron *B = EVALUE_DOMAIN(e->x.p->arr[2*j]);
101 Polyhedron *I = DomainIntersection(A, B, 0);
102 assert(emptyQ(I));
103 Polyhedron_Free(I);
104 }
105 }
106 }
107
108 static int test_eadd(struct barvinok_options *options)
109 {
110 unsigned nvar, nparam;
111 const char **all_vars;
112 evalue *e1, *e2;
113
114 e1 = evalue_read_from_str(" d -1 = 0\n"
115 " h -3 >= 0\n"
116 " - h + 100 >= 0\n"
117 "\n"
118 "(1)\n",
119 "d,h", &all_vars, &nvar, &nparam,
120 options->MaxRays);
121 Free_ParamNames(all_vars, nvar+nparam);
122 e2 = evalue_read_from_str(
123 " h -3 = 0\n"
124 " d -1 >= 0\n"
125 " - d + 3 >= 0\n"
126 "\n"
127 "(0)\n"
128 " d -1 >= 0\n"
129 " - d + 3 >= 0\n"
130 " h -4 >= 0\n"
131 " - h + 100 >= 0\n"
132 "\n"
133 "(1)\n",
134 "d,h", &all_vars, &nvar, &nparam,
135 options->MaxRays);
136 Free_ParamNames(all_vars, nvar+nparam);
137 eadd(e2, e1);
138 evalue_check_disjoint(e1);
139 evalue_free(e1);
140 evalue_free(e2);
141 return 0;
142 }
143
144 int test_evalue(struct barvinok_options *options)
145 {
146 unsigned nvar, nparam;
147 const char **all_vars;
148 evalue *poly1, poly2;
149
150 poly1 = evalue_read_from_str("(1/4 * n^4 + 1/2 * n^3 + 1/4 * n^2)",
151 NULL, &all_vars, &nvar, &nparam,
152 options->MaxRays);
153 Free_ParamNames(all_vars, nvar+nparam);
154
155 value_init(poly2.d);
156 evalue_copy(&poly2, poly1);
157 evalue_negate(poly1);
158 eadd(&poly2, poly1);
159 reduce_evalue(poly1);
160 assert(EVALUE_IS_ZERO(*poly1));
161 evalue_free(poly1);
162 free_evalue_refs(&poly2);
163 return 0;
164 }
165
166 int test_substitute(struct barvinok_options *options)
167 {
168 unsigned nvar, nparam;
169 const char **all_vars;
170 const char *vars = "a,b";
171 evalue *e1, *e2;
172 evalue *subs[2];
173
174 e1 = evalue_read_from_str("[ { 1/3 * a } = 0 ] * \n"
175 " ([ { 1/5 * b + 2/5 } = 0 ] * 5) + \n"
176 "[ { 1/3 * a } != 0 ] * 42",
177 vars, &all_vars, &nvar, &nparam,
178 options->MaxRays);
179 Free_ParamNames(all_vars, nvar+nparam);
180
181 subs[0] = evalue_read_from_str("(2 * b + 5)",
182 vars, &all_vars, &nvar, &nparam,
183 options->MaxRays);
184 Free_ParamNames(all_vars, nvar+nparam);
185 subs[1] = evalue_read_from_str("(a + 1)",
186 vars, &all_vars, &nvar, &nparam,
187 options->MaxRays);
188 Free_ParamNames(all_vars, nvar+nparam);
189
190 evalue_substitute(e1, subs);
191 evalue_free(subs[0]);
192 evalue_free(subs[1]);
193 reduce_evalue(e1);
194
195 e2 = evalue_read_from_str("[ { 2/3 * b + 2/3 } = 0 ] * \n"
196 " ([ { 1/5 * a + 3/5 } = 0 ] * 5) + \n"
197 "[ { 2/3 * b + 2/3 } != 0 ] * 42",
198 vars, &all_vars, &nvar, &nparam,
199 options->MaxRays);
200 Free_ParamNames(all_vars, nvar+nparam);
201 reduce_evalue(e2);
202
203 assert(eequal(e1, e2));
204
205 evalue_free(e1);
206 evalue_free(e2);
207 return 0;
208 }
209
210 int test_specialization(struct barvinok_options *options)
211 {
212 Value v;
213 value_init(v);
214 mpq_t count;
215 mpq_init(count);
216
217 value_set_si(v, 5);
218 dpoly n(2, v);
219 assert(value_cmp_si(n.coeff->p[0], 1) == 0);
220 assert(value_cmp_si(n.coeff->p[1], 5) == 0);
221 assert(value_cmp_si(n.coeff->p[2], 10) == 0);
222
223 value_set_si(v, 1);
224 dpoly d(2, v, 1);
225 value_set_si(v, 2);
226 dpoly d2(2, v, 1);
227 d *= d2;
228 assert(value_cmp_si(d.coeff->p[0], 2) == 0);
229 assert(value_cmp_si(d.coeff->p[1], 1) == 0);
230 assert(value_cmp_si(d.coeff->p[2], 0) == 0);
231
232 n.div(d, count, 1);
233 mpq_canonicalize(count);
234 assert(value_cmp_si(mpq_numref(count), 31) == 0);
235 assert(value_cmp_si(mpq_denref(count), 8) == 0);
236
237 value_set_si(v, -2);
238 dpoly n2(2, v);
239 assert(value_cmp_si(n2.coeff->p[0], 1) == 0);
240 assert(value_cmp_si(n2.coeff->p[1], -2) == 0);
241 assert(value_cmp_si(n2.coeff->p[2], 3) == 0);
242
243 n2.div(d, count, 1);
244 mpq_canonicalize(count);
245 assert(value_cmp_si(mpq_numref(count), 6) == 0);
246 assert(value_cmp_si(mpq_denref(count), 1) == 0);
247
248 value_clear(v);
249 mpq_clear(count);
250 return 0;
251 }
252
253 int test_lattice_points(struct barvinok_options *options)
254 {
255 Param_Vertices V;
256 mat_ZZ tmp;
257 set_from_string(tmp, "[[0 0 0 0 4][0 0 0 0 4][-1 0 1 0 4]]");
258 V.Vertex = zz2matrix(tmp);
259 vec_ZZ lambda;
260 set_from_string(lambda, "[3 5 7]");
261 mat_ZZ rays;
262 set_from_string(rays, "[[0 1 0][4 0 1][0 0 -1]]");
263 term_info num;
264 evalue *point[4];
265
266 unsigned nvar, nparam;
267 const char **all_vars;
268 point[0] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
269 "( 7 * {( 1/4 * a + ( 3/4 * c + 3/4 ) ) } + -21/4 ) ) )",
270 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
271 Free_ParamNames(all_vars, nvar+nparam);
272 point[1] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
273 "( 7 * {( 1/4 * a + ( 3/4 * c + 1/2 ) ) } + -1/2 ) ) )",
274 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
275 Free_ParamNames(all_vars, nvar+nparam);
276 point[2] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
277 "( 7 * {( 1/4 * a + ( 3/4 * c + 1/4 ) ) } + 17/4 ) ) )",
278 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
279 Free_ParamNames(all_vars, nvar+nparam);
280 point[3] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
281 "( 7 * {( 1/4 * a + ( 3/4 * c + 0 ) ) } + 9 ) ) )",
282 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
283 Free_ParamNames(all_vars, nvar+nparam);
284
285 lattice_point(&V, rays, lambda, &num, 4, options);
286 Matrix_Free(V.Vertex);
287
288 for (int i = 0; i < 4; ++i) {
289 assert(eequal(num.E[i], point[i]));
290 evalue_free(point[i]);
291 evalue_free(num.E[i]);
292 }
293 delete [] num.E;
294 return 0;
295 }
296
297 static int test_icounter(struct barvinok_options *options)
298 {
299 icounter cnt(2);
300 vec_QQ n_coeff;
301 mat_ZZ n_power;
302 mat_ZZ d_power;
303 set_from_string(n_coeff, "[-2/1 1/1]");
304 set_from_string(n_power, "[[2 6][3 6]]");
305 d_power.SetDims(0, 2);
306 cnt.reduce(n_coeff, n_power, d_power);
307 assert(value_cmp_si(mpq_numref(cnt.count), -1) == 0);
308 assert(value_cmp_si(mpq_denref(cnt.count), 1) == 0);
309 return 0;
310 }
311
312 static int test_infinite_counter(struct barvinok_options *options)
313 {
314 Matrix *M = matrix_read_from_str("1 3\n 1 1 0\n");
315 Polyhedron *ctx = Constraints2Polyhedron(M, options->MaxRays);
316 Matrix_Free(M);
317
318 /* (1 -1/2 x^5 - 1/2 x^7)/(1-x) */
319 infinite_counter *cnt = new infinite_counter(1, 1);
320 cnt->init(ctx, 0);
321 vec_QQ n_coeff;
322 mat_ZZ n_power;
323 mat_ZZ d_power;
324 set_from_string(n_coeff, "[1/1 -1/2 -1/2]");
325 set_from_string(n_power, "[[0][5][7]]");
326 set_from_string(d_power, "[[1]]");
327 cnt->reduce(n_coeff, n_power, d_power);
328 assert(value_cmp_si(mpq_numref(cnt->count[0]), 6) == 0);
329 assert(value_cmp_si(mpq_denref(cnt->count[0]), 1) == 0);
330 assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) == 0);
331 assert(value_cmp_si(mpq_denref(cnt->count[1]), 1) == 0);
332 delete cnt;
333 Polyhedron_Free(ctx);
334
335 M = matrix_read_from_str("2 4\n 1 1 0 0\n 1 0 1 0\n");
336 ctx = Constraints2Polyhedron(M, options->MaxRays);
337 Matrix_Free(M);
338
339 /* (1 - xy)/((1-x)(1-xy)) */
340 cnt = new infinite_counter(2, 3);
341 cnt->init(ctx, 0);
342 set_from_string(n_coeff, "[1/1 -1/1]");
343 set_from_string(n_power, "[[0 0][1 1]]");
344 set_from_string(d_power, "[[1 0][1 1]]");
345 cnt->reduce(n_coeff, n_power, d_power);
346 assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) != 0);
347 assert(value_cmp_si(mpq_numref(cnt->count[2]), 0) == 0);
348 assert(value_cmp_si(mpq_denref(cnt->count[2]), 1) == 0);
349 assert(value_cmp_si(mpq_numref(cnt->count[3]), 0) == 0);
350 assert(value_cmp_si(mpq_denref(cnt->count[3]), 1) == 0);
351 delete cnt;
352
353 cnt = new infinite_counter(2, 2);
354 cnt->init(ctx, 0);
355 set_from_string(n_coeff, "[-1/2 1/1 -1/3]");
356 set_from_string(n_power, "[[2 6][3 6]]");
357 d_power.SetDims(0, 2);
358 cnt->reduce(n_coeff, n_power, d_power);
359 assert(value_cmp_si(mpq_numref(cnt->count[0]), 1) == 0);
360 assert(value_cmp_si(mpq_denref(cnt->count[0]), 6) == 0);
361 assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) == 0);
362 assert(value_cmp_si(mpq_denref(cnt->count[1]), 1) == 0);
363 assert(value_cmp_si(mpq_numref(cnt->count[2]), 0) == 0);
364 assert(value_cmp_si(mpq_denref(cnt->count[2]), 1) == 0);
365 delete cnt;
366
367 cnt = new infinite_counter(2, 2);
368 cnt->init(ctx, 0);
369 set_from_string(n_coeff, "[1/1]");
370 set_from_string(n_power, "[[0 11]]");
371 set_from_string(d_power, "[[0 1]]");
372 cnt->reduce(n_coeff, n_power, d_power);
373 assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) != 0);
374 assert(value_cmp_si(mpq_numref(cnt->count[2]), 0) == 0);
375 assert(value_cmp_si(mpq_denref(cnt->count[2]), 1) == 0);
376 delete cnt;
377
378 Polyhedron_Free(ctx);
379
380 return 0;
381 }
382
383 static int test_series(struct barvinok_options *options)
384 {
385 Matrix *M;
386 Polyhedron *P, *C;
387 gen_fun *gf;
388
389 M = matrix_read_from_str(
390 "2 3\n"
391 "1 1 0\n"
392 "1 -1 10\n");
393 P = Constraints2Polyhedron(M, options->MaxRays);
394 Matrix_Free(M);
395 C = Universe_Polyhedron(1);
396 gf = barvinok_series_with_options(P, C, options);
397 Polyhedron_Free(P);
398 Polyhedron_Free(C);
399 gen_fun *sum = gf->summate(1, options);
400 delete gf;
401 delete sum;
402
403 return 0;
404 }
405
406 int test_todd(struct barvinok_options *options)
407 {
408 tcounter t(2, options->max_index);
409 assert(value_cmp_si(t.todd.coeff->p[0], 1) == 0);
410 assert(value_cmp_si(t.todd.coeff->p[1], -3) == 0);
411 assert(value_cmp_si(t.todd.coeff->p[2], 3) == 0);
412 assert(value_cmp_si(t.todd_denom->p[0], 1) == 0);
413 assert(value_cmp_si(t.todd_denom->p[1], 6) == 0);
414 assert(value_cmp_si(t.todd_denom->p[2], 36) == 0);
415
416 vec_ZZ lambda;
417 set_from_string(lambda, "[1 -1]");
418 zz2values(lambda, t.lambda->p);
419
420 mat_ZZ rays;
421 set_from_string(rays, "[[-1 0][-1 1]]");
422
423 QQ one(1, 1);
424
425 vec_ZZ v;
426 set_from_string(v, "[2 0 1]");
427 Vector *vertex = Vector_Alloc(3);
428 zz2values(v, vertex->p);
429
430 t.handle(rays, vertex->p, one, 1, options);
431 assert(value_cmp_si(mpq_numref(t.count), 71) == 0);
432 assert(value_cmp_si(mpq_denref(t.count), 24) == 0);
433
434 set_from_string(rays, "[[0 -1][1 -1]]");
435 set_from_string(v, "[0 2 1]");
436 zz2values(v, vertex->p);
437
438 t.handle(rays, vertex->p, one, 1, options);
439 assert(value_cmp_si(mpq_numref(t.count), 71) == 0);
440 assert(value_cmp_si(mpq_denref(t.count), 12) == 0);
441
442 set_from_string(rays, "[[1 0][0 1]]");
443 set_from_string(v, "[0 0 1]");
444 zz2values(v, vertex->p);
445
446 t.handle(rays, vertex->p, one, 1, options);
447 assert(value_cmp_si(mpq_numref(t.count), 6) == 0);
448 assert(value_cmp_si(mpq_denref(t.count), 1) == 0);
449
450 Vector_Free(vertex);
451 return 0;
452 }
453
454 int test_bernoulli(struct barvinok_options *options)
455 {
456 struct bernoulli_coef *bernoulli_coef;
457 struct poly_list *faulhaber, *bernoulli;
458 bernoulli_coef = bernoulli_coef_compute(2);
459 faulhaber = faulhaber_compute(4);
460 bernoulli_coef = bernoulli_coef_compute(8);
461 assert(value_cmp_si(bernoulli_coef->num->p[6], 1) == 0);
462 assert(value_cmp_si(bernoulli_coef->den->p[6], 42) == 0);
463 assert(value_cmp_si(faulhaber->poly[3]->p[0], 0) == 0);
464 assert(value_cmp_si(faulhaber->poly[3]->p[1], 0) == 0);
465 assert(value_cmp_si(faulhaber->poly[3]->p[2], 1) == 0);
466 assert(value_cmp_si(faulhaber->poly[3]->p[3], -2) == 0);
467 assert(value_cmp_si(faulhaber->poly[3]->p[4], 1) == 0);
468 assert(value_cmp_si(faulhaber->poly[3]->p[5], 4) == 0);
469
470 bernoulli = bernoulli_compute(6);
471 assert(value_cmp_si(bernoulli->poly[6]->p[0], 1) == 0);
472 assert(value_cmp_si(bernoulli->poly[6]->p[1], 0) == 0);
473 assert(value_cmp_si(bernoulli->poly[6]->p[2], -21) == 0);
474 assert(value_cmp_si(bernoulli->poly[6]->p[3], 0) == 0);
475 assert(value_cmp_si(bernoulli->poly[6]->p[4], 105) == 0);
476 assert(value_cmp_si(bernoulli->poly[6]->p[5], -126) == 0);
477 assert(value_cmp_si(bernoulli->poly[6]->p[6], 42) == 0);
478 assert(value_cmp_si(bernoulli->poly[6]->p[7], 42) == 0);
479
480 unsigned nvar, nparam;
481 const char **all_vars;
482 evalue *base, *sum1, *sum2;
483 base = evalue_read_from_str("(1 * n + 1)", NULL, &all_vars, &nvar, &nparam,
484 options->MaxRays);
485
486 sum1 = evalue_polynomial(faulhaber->poly[3], base);
487 Free_ParamNames(all_vars, nvar+nparam);
488
489 sum2 = evalue_read_from_str("(1/4 * n^4 + 1/2 * n^3 + 1/4 * n^2)",
490 NULL, &all_vars, &nvar, &nparam,
491 options->MaxRays);
492 Free_ParamNames(all_vars, nvar+nparam);
493 assert(eequal(sum1, sum2));
494 evalue_free(base);
495 evalue_free(sum1);
496 evalue_free(sum2);
497 return 0;
498 }
499
500 int test_bernoulli_sum(struct barvinok_options *options)
501 {
502 int summation = options->summation;
503 options->summation = BV_SUM_BERNOULLI;
504
505 unsigned nvar, nparam;
506 const char **all_vars;
507 evalue *e, *sum1, *sum2;
508 e = evalue_read_from_str("i + -1 >= 0\n -i + n >= 0\n\n 1 + (-1 *i) + i^2",
509 "i", &all_vars, &nvar, &nparam,
510 options->MaxRays);
511 Free_ParamNames(all_vars, nvar+nparam);
512
513 sum1 = barvinok_summate(e, 1, options);
514 sum2 = evalue_read_from_str("n -1 >= 0\n\n (1/3 * n^3 + 2/3 * n)",
515 NULL, &all_vars, &nvar, &nparam,
516 options->MaxRays);
517 Free_ParamNames(all_vars, nvar+nparam);
518 evalue_negate(sum1);
519 eadd(sum2, sum1);
520 reduce_evalue(sum1);
521 assert(EVALUE_IS_ZERO(*sum1));
522 evalue_free(e);
523 evalue_free(sum1);
524
525 e = evalue_read_from_str("-i + -1 >= 0\n i + n >= 0\n\n 1 + i + i^2",
526 "i", &all_vars, &nvar, &nparam,
527 options->MaxRays);
528 Free_ParamNames(all_vars, nvar+nparam);
529 sum1 = barvinok_summate(e, 1, options);
530 evalue_negate(sum1);
531 eadd(sum2, sum1);
532 reduce_evalue(sum1);
533 assert(EVALUE_IS_ZERO(*sum1));
534 evalue_free(e);
535 evalue_free(sum1);
536
537 evalue_free(sum2);
538
539 e = evalue_read_from_str("i + 4 >= 0\n -i + n >= 0\n\n i",
540 "i", &all_vars, &nvar, &nparam,
541 options->MaxRays);
542 Free_ParamNames(all_vars, nvar+nparam);
543 sum1 = barvinok_summate(e, 1, options);
544 sum2 = evalue_read_from_str("n + 0 >= 0\n\n (1/2 * n^2 + 1/2 * n + (-10))\n"
545 "n + 4 >= 0\n -n -1 >= 0\n\n (1/2 * n^2 + 1/2 * n + (-10))",
546 NULL, &all_vars, &nvar, &nparam,
547 options->MaxRays);
548 Free_ParamNames(all_vars, nvar+nparam);
549 evalue_negate(sum1);
550 eadd(sum2, sum1);
551 reduce_evalue(sum1);
552 assert(EVALUE_IS_ZERO(*sum1));
553 evalue_free(e);
554 evalue_free(sum1);
555 evalue_free(sum2);
556
557 e = evalue_read_from_str("i -5 >= 0\n -i + n >= 0\n j -1 >= 0\n i -j >= 0\n"
558 "k -1 >= 0\n j -k >= 0\n\n1",
559 "i,j,k", &all_vars, &nvar, &nparam,
560 options->MaxRays);
561 Free_ParamNames(all_vars, nvar+nparam);
562 sum1 = barvinok_summate(e, 3, options);
563 sum2 = evalue_read_from_str("n -5 >= 0\n\n"
564 "1/6 * n^3 + 1/2 * n^2 + 1/3 * n + -20",
565 NULL, &all_vars, &nvar, &nparam,
566 options->MaxRays);
567 Free_ParamNames(all_vars, nvar+nparam);
568 evalue_negate(sum1);
569 eadd(sum2, sum1);
570 reduce_evalue(sum1);
571 assert(EVALUE_IS_ZERO(*sum1));
572 evalue_free(e);
573 evalue_free(sum1);
574 evalue_free(sum2);
575
576 options->summation = summation;
577 return 0;
578 }
579
580 int test_hilbert(struct barvinok_options *options)
581 {
582 #ifdef USE_ZSOLVE
583 Matrix *M = matrix_read_from_str(
584 "2 4\n"
585 " 1 4 -3 0 \n"
586 " 1 3 2 0 \n");
587 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
588 Matrix_Free(M);
589
590 M = Cone_Hilbert_Basis(P, options->MaxRays);
591 assert(M->NbRows == 5);
592 assert(M->NbColumns == 3);
593 Matrix_Free(M);
594
595 M = Cone_Integer_Hull(P, NULL, 0, options);
596 assert(M->NbRows == 4);
597 assert(M->NbColumns == 3);
598 Matrix_Free(M);
599
600 Polyhedron_Free(P);
601 #endif
602 return 0;
603 }
604
605 int test_ilp(struct barvinok_options *options)
606 {
607 Matrix *M = matrix_read_from_str(
608 "2 4\n"
609 " 1 4 -3 0 \n"
610 " 1 3 2 0 \n");
611 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
612 Matrix_Free(M);
613 Vector *obj = Vector_Alloc(2);
614 value_set_si(obj->p[0], 7);
615 value_set_si(obj->p[1], -1);
616 Value min, max;
617 value_init(min);
618 value_init(max);
619
620 value_set_si(min, 1);
621 value_set_si(max, 17);
622 Vector *opt = Polyhedron_Integer_Minimum(P, obj->p, min, max,
623 NULL, 0, options);
624 assert(opt);
625 assert(value_cmp_si(opt->p[0], 1) == 0);
626 assert(value_cmp_si(opt->p[1], 1) == 0);
627 assert(value_cmp_si(opt->p[2], 1) == 0);
628 Vector_Free(opt);
629
630 value_clear(min);
631 value_clear(max);
632 Vector_Free(obj);
633 Polyhedron_Free(P);
634 return 0;
635 }
636
637 int test_hull(struct barvinok_options *options)
638 {
639 Matrix *M = matrix_read_from_str(
640 "4 4\n"
641 "1 32 -20 7\n"
642 "1 8 -44 187\n"
643 "1 -48 -4 285\n"
644 "1 8 68 -199\n");
645 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
646 Matrix_Free(M);
647
648 Matrix *hull = Polyhedron_Integer_Hull(P, options);
649 Polyhedron_Free(P);
650 assert(hull->NbRows == 4);
651 M = Matrix_Alloc(hull->NbRows, 1+hull->NbColumns);
652 for (int i = 0; i < hull->NbRows; ++i) {
653 value_set_si(M->p[i][0], 1);
654 Vector_Copy(hull->p[i], M->p[i]+1, hull->NbColumns);
655 }
656 Matrix_Free(hull);
657 Polyhedron *H = Constraints2Polyhedron(M, options->MaxRays);
658 Matrix_Free(M);
659
660 M = matrix_read_from_str(
661 "4 4\n"
662 "1 2 3 1 \n"
663 "1 3 4 1 \n"
664 "1 5 3 1 \n"
665 "1 5 5 1 \n");
666 P = Constraints2Polyhedron(M, options->MaxRays);
667 Matrix_Free(M);
668 assert(PolyhedronIncludes(P, H) && PolyhedronIncludes(H, P));
669 Polyhedron_Free(P);
670 Polyhedron_Free(H);
671
672 M = matrix_read_from_str(
673 "3 4\n"
674 "1 2 6 -3 \n"
675 "1 2 -6 3 \n"
676 "1 -2 0 3 \n");
677 P = Constraints2Polyhedron(M, options->MaxRays);
678 Matrix_Free(M);
679 assert(!emptyQ(P));
680 hull = Polyhedron_Integer_Hull(P, options);
681 Polyhedron_Free(P);
682 assert(hull->NbRows == 0);
683 Matrix_Free(hull);
684 return 0;
685 }
686
687 static int test_laurent(struct barvinok_options *options)
688 {
689 unsigned nvar, nparam;
690 const char **all_vars;
691 evalue *e, *sum, *res;
692
693 e = evalue_read_from_str(" x1 >= 0\n"
694 " x2 >= 0\n"
695 " -x1 -x2 + 2 >= 0\n"
696 "\n"
697 "(N + M * x2)\n",
698 "x1,x2", &all_vars, &nvar, &nparam,
699 options->MaxRays);
700 Free_ParamNames(all_vars, nvar+nparam);
701
702 int summation = options->summation;
703 options->summation = BV_SUM_LAURENT;
704 sum = barvinok_summate(e, nvar, options);
705 options->summation = summation;
706
707 res = evalue_read_from_str("(6 * N + 4 * M)",
708 "", &all_vars, &nvar, &nparam,
709 options->MaxRays);
710 Free_ParamNames(all_vars, nvar+nparam);
711
712 assert(value_zero_p(sum->d));
713 assert(sum->x.p->type == partition);
714 assert(sum->x.p->size == 2);
715
716 assert(eequal(res, &sum->x.p->arr[1]));
717
718 evalue_free(e);
719 evalue_free(sum);
720 evalue_free(res);
721 return 0;
722 }
723
724 /* Check that Polyhedron_Reduced_Basis produces a result
725 * of the expected dimensions (without crashing).
726 */
727 static int test_basis_reduction(struct barvinok_options *options)
728 {
729 Matrix *M;
730 Polyhedron *P;
731
732 M = matrix_read_from_str(
733 "4 4\n"
734 "1 1 0 0 \n"
735 "1 0 1 0 \n"
736 "1 -1 0 1 \n"
737 "1 0 -1 1 \n");
738 P = Constraints2Polyhedron(M, options->MaxRays);
739 Matrix_Free(M);
740
741 M = Polyhedron_Reduced_Basis(P, options);
742
743 assert(M);
744 assert(M->NbRows == 2);
745 assert(M->NbColumns == 2);
746
747 Polyhedron_Free(P);
748 Matrix_Free(M);
749
750 return 0;
751 }
752
753 int main(int argc, char **argv)
754 {
755 struct barvinok_options *options = barvinok_options_new_with_defaults();
756 test_equalities(options);
757 test_evalue_read(options);
758 test_eadd(options);
759 test_evalue(options);
760 test_substitute(options);
761 test_specialization(options);
762 test_lattice_points(options);
763 test_icounter(options);
764 test_infinite_counter(options);
765 test_series(options);
766 test_todd(options);
767 test_bernoulli(options);
768 test_bernoulli_sum(options);
769 test_hilbert(options);
770 test_ilp(options);
771 test_hull(options);
772 test_laurent(options);
773 test_basis_reduction(options);
774 barvinok_options_free(options);
775
776 return EXIT_SUCCESS;
777 }
lattice point enumeration library