?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
doc: clean up "exponential substitution" section
[barvinok/uuh.git] / testlib.cc
blob88345b0ad35b9f7b604a1ee34205f83bca35737c
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 "remove_equalities.h"
19 #include "config.h"
20
21 using std::cout;
22 using std::cerr;
23 using std::endl;
24
25 template <typename T>
26 void set_from_string(T& v, const char *s)
27 {
28 std::istringstream str(s);
29 str >> v;
30 }
31
32 static Matrix *matrix_read_from_str(const char *s)
33 {
34 std::istringstream str(s);
35 return Matrix_Read(str);
36 }
37
38 static int test_equalities(struct barvinok_options *options)
39 {
40 Matrix *M = matrix_read_from_str(
41 "11 11\n"
42 " 0 23 0 0 -10 0 0 0 7 -44 -8 \n"
43 " 0 0 23 0 5 0 0 0 -15 114 27 \n"
44 " 0 0 0 1 0 0 0 0 0 -1 2 \n"
45 " 0 0 0 0 0 1 0 0 -1 8 0 \n"
46 " 0 0 0 0 0 0 1 0 0 -1 2 \n"
47 " 0 0 0 0 0 0 0 1 0 -1 -1 \n"
48 " 1 0 0 0 10 0 0 0 -7 44 8 \n"
49 " 1 0 0 0 -5 0 0 0 15 -114 -27 \n"
50 " 1 0 0 0 1 0 0 0 0 0 0 \n"
51 " 1 0 0 0 0 0 0 0 1 -8 0 \n"
52 " 1 0 0 0 0 0 0 0 0 1 -2 \n");
53 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
54 Matrix_Free(M);
55 Polyhedron *Q = P;
56 remove_all_equalities(&P, NULL, NULL, NULL, 2, options->MaxRays);
57 assert(P->NbEq == 0);
58 if (P != Q)
59 Polyhedron_Free(Q);
60 Polyhedron_Free(P);
61 }
62
63 int test_evalue_read(struct barvinok_options *options)
64 {
65 unsigned nvar, nparam;
66 char **all_vars;
67 evalue *e1, *e2;
68
69 e1 = evalue_read_from_str("(1 * aa + 2 * a)",
70 NULL, &all_vars, &nvar, &nparam, options->MaxRays);
71 Free_ParamNames(all_vars, nvar+nparam);
72 e2 = evalue_read_from_str("(3 * aa)",
73 NULL, &all_vars, &nvar, &nparam, options->MaxRays);
74 Free_ParamNames(all_vars, nvar+nparam);
75 assert(!eequal(e1, e2));
76 evalue_free(e1);
77 evalue_free(e2);
78 }
79
80 static void evalue_check_disjoint(evalue *e)
81 {
82 int i, j;
83
84 if (!e)
85 return;
86 if (EVALUE_IS_ZERO(*e))
87 return;
88 for (i = 0; i < e->x.p->size/2; ++i) {
89 Polyhedron *A = EVALUE_DOMAIN(e->x.p->arr[2*i]);
90 for (j = i+1; j < e->x.p->size/2; ++j) {
91 Polyhedron *B = EVALUE_DOMAIN(e->x.p->arr[2*j]);
92 Polyhedron *I = DomainIntersection(A, B, 0);
93 assert(emptyQ(I));
94 Polyhedron_Free(I);
95 }
96 }
97 }
98
99 static int test_eadd(struct barvinok_options *options)
100 {
101 unsigned nvar, nparam;
102 char **all_vars;
103 evalue *e1, *e2;
104
105 e1 = evalue_read_from_str(" d -1 = 0\n"
106 " h -3 >= 0\n"
107 " - h + 100 >= 0\n"
108 "\n"
109 "(1)\n",
110 "d,h", &all_vars, &nvar, &nparam,
111 options->MaxRays);
112 Free_ParamNames(all_vars, nvar+nparam);
113 e2 = evalue_read_from_str(
114 " h -3 = 0\n"
115 " d -1 >= 0\n"
116 " - d + 3 >= 0\n"
117 "\n"
118 "(0)\n"
119 " d -1 >= 0\n"
120 " - d + 3 >= 0\n"
121 " h -4 >= 0\n"
122 " - h + 100 >= 0\n"
123 "\n"
124 "(1)\n",
125 "d,h", &all_vars, &nvar, &nparam,
126 options->MaxRays);
127 Free_ParamNames(all_vars, nvar+nparam);
128 eadd(e2, e1);
129 evalue_check_disjoint(e1);
130 evalue_free(e1);
131 evalue_free(e2);
132 }
133
134 int test_evalue(struct barvinok_options *options)
135 {
136 unsigned nvar, nparam;
137 char **all_vars;
138 evalue *poly1, poly2;
139
140 poly1 = evalue_read_from_str("(1/4 * n^4 + 1/2 * n^3 + 1/4 * n^2)",
141 NULL, &all_vars, &nvar, &nparam,
142 options->MaxRays);
143 Free_ParamNames(all_vars, nvar+nparam);
144
145 value_init(poly2.d);
146 evalue_copy(&poly2, poly1);
147 evalue_negate(poly1);
148 eadd(&poly2, poly1);
149 reduce_evalue(poly1);
150 assert(EVALUE_IS_ZERO(*poly1));
151 evalue_free(poly1);
152 free_evalue_refs(&poly2);
153 }
154
155 int test_substitute(struct barvinok_options *options)
156 {
157 unsigned nvar, nparam;
158 char **all_vars;
159 const char *vars = "a,b";
160 evalue *e1, *e2;
161 evalue *subs[2];
162
163 e1 = evalue_read_from_str("[ { 1/3 * a } = 0 ] * \n"
164 " ([ { 1/5 * b + 2/5 } = 0 ] * 5) + \n"
165 "[ { 1/3 * a } != 0 ] * 42",
166 vars, &all_vars, &nvar, &nparam,
167 options->MaxRays);
168 Free_ParamNames(all_vars, nvar+nparam);
169
170 subs[0] = evalue_read_from_str("(2 * b + 5)",
171 vars, &all_vars, &nvar, &nparam,
172 options->MaxRays);
173 Free_ParamNames(all_vars, nvar+nparam);
174 subs[1] = evalue_read_from_str("(a + 1)",
175 vars, &all_vars, &nvar, &nparam,
176 options->MaxRays);
177 Free_ParamNames(all_vars, nvar+nparam);
178
179 evalue_substitute(e1, subs);
180 evalue_free(subs[0]);
181 evalue_free(subs[1]);
182 reduce_evalue(e1);
183
184 e2 = evalue_read_from_str("[ { 2/3 * b + 2/3 } = 0 ] * \n"
185 " ([ { 1/5 * a + 3/5 } = 0 ] * 5) + \n"
186 "[ { 2/3 * b + 2/3 } != 0 ] * 42",
187 vars, &all_vars, &nvar, &nparam,
188 options->MaxRays);
189 Free_ParamNames(all_vars, nvar+nparam);
190 reduce_evalue(e2);
191
192 assert(eequal(e1, e2));
193
194 evalue_free(e1);
195 evalue_free(e2);
196 }
197
198 int test_split_periods(struct barvinok_options *options)
199 {
200 unsigned nvar, nparam;
201 char **all_vars;
202 evalue *e;
203
204 e = evalue_read_from_str("U + 2V + 3 >= 0\n- U -2V >= 0\n- U 10 >= 0\n"
205 "U >= 0\n\n({( 1/3 * U + ( 2/3 * V + 0 ))})",
206 NULL, &all_vars, &nvar, &nparam,
207 options->MaxRays);
208 Free_ParamNames(all_vars, nvar+nparam);
209
210 evalue_split_periods(e, 2, options->MaxRays);
211 assert(value_zero_p(e->d));
212 assert(e->x.p->type == partition);
213 assert(e->x.p->size == 4);
214 assert(value_zero_p(e->x.p->arr[1].d));
215 assert(e->x.p->arr[1].x.p->type == polynomial);
216 assert(value_zero_p(e->x.p->arr[3].d));
217 assert(e->x.p->arr[3].x.p->type == polynomial);
218 evalue_free(e);
219 }
220
221 int test_specialization(struct barvinok_options *options)
222 {
223 Value v;
224 value_init(v);
225 mpq_t count;
226 mpq_init(count);
227
228 value_set_si(v, 5);
229 dpoly n(2, v);
230 assert(value_cmp_si(n.coeff->p[0], 1) == 0);
231 assert(value_cmp_si(n.coeff->p[1], 5) == 0);
232 assert(value_cmp_si(n.coeff->p[2], 10) == 0);
233
234 value_set_si(v, 1);
235 dpoly d(2, v, 1);
236 value_set_si(v, 2);
237 dpoly d2(2, v, 1);
238 d *= d2;
239 assert(value_cmp_si(d.coeff->p[0], 2) == 0);
240 assert(value_cmp_si(d.coeff->p[1], 1) == 0);
241 assert(value_cmp_si(d.coeff->p[2], 0) == 0);
242
243 n.div(d, count, 1);
244 mpq_canonicalize(count);
245 assert(value_cmp_si(mpq_numref(count), 31) == 0);
246 assert(value_cmp_si(mpq_denref(count), 8) == 0);
247
248 value_set_si(v, -2);
249 dpoly n2(2, v);
250 assert(value_cmp_si(n2.coeff->p[0], 1) == 0);
251 assert(value_cmp_si(n2.coeff->p[1], -2) == 0);
252 assert(value_cmp_si(n2.coeff->p[2], 3) == 0);
253
254 n2.div(d, count, 1);
255 mpq_canonicalize(count);
256 assert(value_cmp_si(mpq_numref(count), 6) == 0);
257 assert(value_cmp_si(mpq_denref(count), 1) == 0);
258
259 value_clear(v);
260 mpq_clear(count);
261 }
262
263 int test_lattice_points(struct barvinok_options *options)
264 {
265 Param_Vertices V;
266 mat_ZZ tmp;
267 set_from_string(tmp, "[[0 0 0 0 4][0 0 0 0 4][-1 0 1 0 4]]");
268 V.Vertex = zz2matrix(tmp);
269 vec_ZZ lambda;
270 set_from_string(lambda, "[3 5 7]");
271 mat_ZZ rays;
272 set_from_string(rays, "[[0 1 0][4 0 1][0 0 -1]]");
273 term_info num;
274 evalue *point[4];
275
276 unsigned nvar, nparam;
277 char **all_vars;
278 point[0] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
279 "( 7 * {( 1/4 * a + ( 3/4 * c + 3/4 ) ) } + -21/4 ) ) )",
280 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
281 Free_ParamNames(all_vars, nvar+nparam);
282 point[1] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
283 "( 7 * {( 1/4 * a + ( 3/4 * c + 1/2 ) ) } + -1/2 ) ) )",
284 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
285 Free_ParamNames(all_vars, nvar+nparam);
286 point[2] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
287 "( 7 * {( 1/4 * a + ( 3/4 * c + 1/4 ) ) } + 17/4 ) ) )",
288 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
289 Free_ParamNames(all_vars, nvar+nparam);
290 point[3] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
291 "( 7 * {( 1/4 * a + ( 3/4 * c + 0 ) ) } + 9 ) ) )",
292 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
293 Free_ParamNames(all_vars, nvar+nparam);
294
295 lattice_point(&V, rays, lambda, &num, 4, options);
296 Matrix_Free(V.Vertex);
297
298 for (int i = 0; i < 4; ++i) {
299 assert(eequal(num.E[i], point[i]));
300 evalue_free(point[i]);
301 free_evalue_refs(num.E[i]);
302 delete 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 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 unsigned nvar, nparam;
545 char **all_vars;
546 evalue *e, *sum1, *sum2;
547 e = evalue_read_from_str("i + -1 >= 0\n -i + n >= 0\n\n 1 + (-1 *i) + i^2",
548 "i", &all_vars, &nvar, &nparam,
549 options->MaxRays);
550 Free_ParamNames(all_vars, nvar+nparam);
551
552 sum1 = Bernoulli_sum_evalue(e, 1, options);
553 sum2 = evalue_read_from_str("n -1 >= 0\n\n (1/3 * n^3 + 2/3 * n)",
554 NULL, &all_vars, &nvar, &nparam,
555 options->MaxRays);
556 Free_ParamNames(all_vars, nvar+nparam);
557 evalue_negate(sum1);
558 eadd(sum2, sum1);
559 reduce_evalue(sum1);
560 assert(EVALUE_IS_ZERO(*sum1));
561 evalue_free(e);
562 evalue_free(sum1);
563
564 e = evalue_read_from_str("-i + -1 >= 0\n i + n >= 0\n\n 1 + i + i^2",
565 "i", &all_vars, &nvar, &nparam,
566 options->MaxRays);
567 Free_ParamNames(all_vars, nvar+nparam);
568 sum1 = Bernoulli_sum_evalue(e, 1, options);
569 evalue_negate(sum1);
570 eadd(sum2, sum1);
571 reduce_evalue(sum1);
572 assert(EVALUE_IS_ZERO(*sum1));
573 evalue_free(e);
574 evalue_free(sum1);
575
576 evalue_free(sum2);
577
578 e = evalue_read_from_str("i + 4 >= 0\n -i + n >= 0\n\n i",
579 "i", &all_vars, &nvar, &nparam,
580 options->MaxRays);
581 Free_ParamNames(all_vars, nvar+nparam);
582 sum1 = Bernoulli_sum_evalue(e, 1, options);
583 sum2 = evalue_read_from_str("n + 0 >= 0\n\n (1/2 * n^2 + 1/2 * n + (-10))\n"
584 "n + 4 >= 0\n -n -1 >= 0\n\n (1/2 * n^2 + 1/2 * n + (-10))",
585 NULL, &all_vars, &nvar, &nparam,
586 options->MaxRays);
587 Free_ParamNames(all_vars, nvar+nparam);
588 evalue_negate(sum1);
589 eadd(sum2, sum1);
590 reduce_evalue(sum1);
591 assert(EVALUE_IS_ZERO(*sum1));
592 evalue_free(e);
593 evalue_free(sum1);
594 evalue_free(sum2);
595
596 e = evalue_read_from_str("i -5 >= 0\n -i + n >= 0\n j -1 >= 0\n i -j >= 0\n"
597 "k -1 >= 0\n j -k >= 0\n\n1",
598 "i,j,k", &all_vars, &nvar, &nparam,
599 options->MaxRays);
600 Free_ParamNames(all_vars, nvar+nparam);
601 sum1 = Bernoulli_sum_evalue(e, 3, options);
602 sum2 = evalue_read_from_str("n -5 >= 0\n\n"
603 "1/6 * n^3 + 1/2 * n^2 + 1/3 * n + -20",
604 NULL, &all_vars, &nvar, &nparam,
605 options->MaxRays);
606 Free_ParamNames(all_vars, nvar+nparam);
607 evalue_negate(sum1);
608 eadd(sum2, sum1);
609 reduce_evalue(sum1);
610 assert(EVALUE_IS_ZERO(*sum1));
611 evalue_free(e);
612 evalue_free(sum1);
613 evalue_free(sum2);
614 }
615
616 int test_hilbert(struct barvinok_options *options)
617 {
618 #ifdef USE_ZSOLVE
619 Matrix *M = matrix_read_from_str(
620 "2 4\n"
621 " 1 4 -3 0 \n"
622 " 1 3 2 0 \n");
623 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
624 Matrix_Free(M);
625
626 M = Cone_Hilbert_Basis(P, options->MaxRays);
627 assert(M->NbRows = 5);
628 assert(M->NbColumns = 3);
629 Matrix_Free(M);
630
631 M = Cone_Integer_Hull(P, NULL, 0, options);
632 assert(M->NbRows = 4);
633 assert(M->NbColumns = 3);
634 Matrix_Free(M);
635
636 Polyhedron_Free(P);
637 #endif
638 }
639
640 int test_ilp(struct barvinok_options *options)
641 {
642 Matrix *M = matrix_read_from_str(
643 "2 4\n"
644 " 1 4 -3 0 \n"
645 " 1 3 2 0 \n");
646 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
647 Matrix_Free(M);
648 Vector *obj = Vector_Alloc(2);
649 value_set_si(obj->p[0], 7);
650 value_set_si(obj->p[1], -1);
651 Value min, max;
652 value_init(min);
653 value_init(max);
654
655 value_set_si(min, 1);
656 value_set_si(max, 17);
657 Vector *opt = Polyhedron_Integer_Minimum(P, obj->p, min, max,
658 NULL, 0, options);
659 assert(opt);
660 assert(value_cmp_si(opt->p[0], 1) == 0);
661 assert(value_cmp_si(opt->p[1], 1) == 0);
662 assert(value_cmp_si(opt->p[2], 1) == 0);
663 Vector_Free(opt);
664
665 value_clear(min);
666 value_clear(max);
667 Vector_Free(obj);
668 Polyhedron_Free(P);
669 }
670
671 int test_hull(struct barvinok_options *options)
672 {
673 Matrix *M = matrix_read_from_str(
674 "4 4\n"
675 "1 32 -20 7\n"
676 "1 8 -44 187\n"
677 "1 -48 -4 285\n"
678 "1 8 68 -199\n");
679 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
680 Matrix_Free(M);
681
682 Matrix *hull = Polyhedron_Integer_Hull(P, options);
683 Polyhedron_Free(P);
684 assert(hull->NbRows == 4);
685 M = Matrix_Alloc(hull->NbRows, 1+hull->NbColumns);
686 for (int i = 0; i < hull->NbRows; ++i) {
687 value_set_si(M->p[i][0], 1);
688 Vector_Copy(hull->p[i], M->p[i]+1, hull->NbColumns);
689 }
690 Matrix_Free(hull);
691 Polyhedron *H = Constraints2Polyhedron(M, options->MaxRays);
692 Matrix_Free(M);
693
694 M = matrix_read_from_str(
695 "4 4\n"
696 "1 2 3 1 \n"
697 "1 3 4 1 \n"
698 "1 5 3 1 \n"
699 "1 5 5 1 \n");
700 P = Constraints2Polyhedron(M, options->MaxRays);
701 Matrix_Free(M);
702 assert(PolyhedronIncludes(P, H) && PolyhedronIncludes(H, P));
703 Polyhedron_Free(P);
704 Polyhedron_Free(H);
705
706 M = matrix_read_from_str(
707 "3 4\n"
708 "1 2 6 -3 \n"
709 "1 2 -6 3 \n"
710 "1 -2 0 3 \n");
711 P = Constraints2Polyhedron(M, options->MaxRays);
712 Matrix_Free(M);
713 assert(!emptyQ(P));
714 hull = Polyhedron_Integer_Hull(P, options);
715 Polyhedron_Free(P);
716 assert(hull->NbRows == 0);
717 Matrix_Free(hull);
718 }
719
720 int main(int argc, char **argv)
721 {
722 struct barvinok_options *options = barvinok_options_new_with_defaults();
723 test_equalities(options);
724 test_evalue_read(options);
725 test_eadd(options);
726 test_evalue(options);
727 test_substitute(options);
728 test_split_periods(options);
729 test_specialization(options);
730 test_lattice_points(options);
731 test_icounter(options);
732 test_infinite_counter(options);
733 test_series(options);
734 test_todd(options);
735 test_bernoulli(options);
736 test_bernoulli_sum(options);
737 test_hilbert(options);
738 test_ilp(options);
739 test_hull(options);
740 barvinok_options_free(options);
741 }
patches by uuh to consider upstream