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