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