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