?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
bernstein: add piecewise_lst::is_equal
[barvinok.git] / testlib.cc
blob12f4a7c5b19ea27b2bd0ab1403abc63308c74abe
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 "matrix_read.h"
16
17 using std::cout;
18 using std::cerr;
19 using std::endl;
20
21 template <typename T>
22 void set_from_string(T& v, const char *s)
23 {
24 std::istringstream str(s);
25 str >> v;
26 }
27
28 int test_evalue(struct barvinok_options *options)
29 {
30 unsigned nvar, nparam;
31 char **all_vars;
32 evalue *poly1, poly2;
33
34 poly1 = evalue_read_from_str("(1/4 * n^4 + 1/2 * n^3 + 1/4 * n^2)",
35 NULL, &all_vars, &nvar, &nparam,
36 options->MaxRays);
37 Free_ParamNames(all_vars, nvar+nparam);
38
39 value_init(poly2.d);
40 evalue_copy(&poly2, poly1);
41 evalue_negate(poly1);
42 eadd(&poly2, poly1);
43 reduce_evalue(poly1);
44 assert(EVALUE_IS_ZERO(*poly1));
45 evalue_free(poly1);
46 free_evalue_refs(&poly2);
47 }
48
49 int test_split_periods(struct barvinok_options *options)
50 {
51 unsigned nvar, nparam;
52 char **all_vars;
53 evalue *e;
54
55 e = evalue_read_from_str("U + 2V + 3 >= 0\n- U -2V >= 0\n- U 10 >= 0\n"
56 "U >= 0\n\n({( 1/3 * U + ( 2/3 * V + 0 ))})",
57 NULL, &all_vars, &nvar, &nparam,
58 options->MaxRays);
59 Free_ParamNames(all_vars, nvar+nparam);
60
61 evalue_split_periods(e, 2, options->MaxRays);
62 assert(value_zero_p(e->d));
63 assert(e->x.p->type == partition);
64 assert(e->x.p->size == 4);
65 assert(value_zero_p(e->x.p->arr[1].d));
66 assert(e->x.p->arr[1].x.p->type == polynomial);
67 assert(value_zero_p(e->x.p->arr[3].d));
68 assert(e->x.p->arr[3].x.p->type == polynomial);
69 evalue_free(e);
70 }
71
72 int test_specialization(struct barvinok_options *options)
73 {
74 Value v;
75 value_init(v);
76 mpq_t count;
77 mpq_init(count);
78
79 value_set_si(v, 5);
80 dpoly n(2, v);
81 assert(value_cmp_si(n.coeff->p[0], 1) == 0);
82 assert(value_cmp_si(n.coeff->p[1], 5) == 0);
83 assert(value_cmp_si(n.coeff->p[2], 10) == 0);
84
85 value_set_si(v, 1);
86 dpoly d(2, v, 1);
87 value_set_si(v, 2);
88 dpoly d2(2, v, 1);
89 d *= d2;
90 assert(value_cmp_si(d.coeff->p[0], 2) == 0);
91 assert(value_cmp_si(d.coeff->p[1], 1) == 0);
92 assert(value_cmp_si(d.coeff->p[2], 0) == 0);
93
94 n.div(d, count, 1);
95 mpq_canonicalize(count);
96 assert(value_cmp_si(mpq_numref(count), 31) == 0);
97 assert(value_cmp_si(mpq_denref(count), 8) == 0);
98
99 value_set_si(v, -2);
100 dpoly n2(2, v);
101 assert(value_cmp_si(n2.coeff->p[0], 1) == 0);
102 assert(value_cmp_si(n2.coeff->p[1], -2) == 0);
103 assert(value_cmp_si(n2.coeff->p[2], 3) == 0);
104
105 n2.div(d, count, 1);
106 mpq_canonicalize(count);
107 assert(value_cmp_si(mpq_numref(count), 6) == 0);
108 assert(value_cmp_si(mpq_denref(count), 1) == 0);
109
110 value_clear(v);
111 mpq_clear(count);
112 }
113
114 int test_lattice_points(struct barvinok_options *options)
115 {
116 Param_Vertices V;
117 mat_ZZ tmp;
118 set_from_string(tmp, "[[0 0 0 0 4][0 0 0 0 4][-1 0 1 0 4]]");
119 V.Vertex = zz2matrix(tmp);
120 vec_ZZ lambda;
121 set_from_string(lambda, "[3 5 7]");
122 mat_ZZ rays;
123 set_from_string(rays, "[[0 1 0][4 0 1][0 0 -1]]");
124 term_info num;
125 evalue *point[4];
126
127 unsigned nvar, nparam;
128 char **all_vars;
129 point[0] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
130 "( 7 * {( 1/4 * a + ( 3/4 * c + 3/4 ) ) } + -21/4 ) ) )",
131 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
132 Free_ParamNames(all_vars, nvar+nparam);
133 point[1] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
134 "( 7 * {( 1/4 * a + ( 3/4 * c + 1/2 ) ) } + -1/2 ) ) )",
135 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
136 Free_ParamNames(all_vars, nvar+nparam);
137 point[2] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
138 "( 7 * {( 1/4 * a + ( 3/4 * c + 1/4 ) ) } + 17/4 ) ) )",
139 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
140 Free_ParamNames(all_vars, nvar+nparam);
141 point[3] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
142 "( 7 * {( 1/4 * a + ( 3/4 * c + 0 ) ) } + 9 ) ) )",
143 "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
144 Free_ParamNames(all_vars, nvar+nparam);
145
146 lattice_point(&V, rays, lambda, &num, 4, options);
147 Matrix_Free(V.Vertex);
148
149 for (int i = 0; i < 4; ++i) {
150 assert(eequal(num.E[i], point[i]));
151 evalue_free(point[i]);
152 free_evalue_refs(num.E[i]);
153 delete num.E[i];
154 }
155 delete [] num.E;
156 }
157
158 static int test_icounter(struct barvinok_options *options)
159 {
160 icounter cnt(2);
161 vec_QQ n_coeff;
162 mat_ZZ n_power;
163 mat_ZZ d_power;
164 set_from_string(n_coeff, "[-2/1 1/1]");
165 set_from_string(n_power, "[[2 6][3 6]]");
166 d_power.SetDims(0, 2);
167 cnt.reduce(n_coeff, n_power, d_power);
168 assert(value_cmp_si(mpq_numref(cnt.count), -1) == 0);
169 assert(value_cmp_si(mpq_denref(cnt.count), 1) == 0);
170 }
171
172 static Matrix *matrix_read_from_str(const char *s)
173 {
174 std::istringstream str(s);
175 return Matrix_Read(str);
176 }
177
178 static int test_infinite_counter(struct barvinok_options *options)
179 {
180 Matrix *M = matrix_read_from_str("1 3\n 1 1 0\n");
181 Polyhedron *ctx = Constraints2Polyhedron(M, options->MaxRays);
182 Matrix_Free(M);
183
184 /* (1 -1/2 x^5 - 1/2 x^7)/(1-x) */
185 infinite_counter *cnt = new infinite_counter(1, 1);
186 cnt->init(ctx);
187 vec_QQ n_coeff;
188 mat_ZZ n_power;
189 mat_ZZ d_power;
190 set_from_string(n_coeff, "[1/1 -1/2 -1/2]");
191 set_from_string(n_power, "[[0][5][7]]");
192 set_from_string(d_power, "[[1]]");
193 cnt->reduce(n_coeff, n_power, d_power);
194 assert(value_cmp_si(mpq_numref(cnt->count[0]), 6) == 0);
195 assert(value_cmp_si(mpq_denref(cnt->count[0]), 1) == 0);
196 assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) == 0);
197 assert(value_cmp_si(mpq_denref(cnt->count[1]), 1) == 0);
198 delete cnt;
199 Polyhedron_Free(ctx);
200
201 M = matrix_read_from_str("2 4\n 1 1 0 0\n 1 0 1 0\n");
202 ctx = Constraints2Polyhedron(M, options->MaxRays);
203 Matrix_Free(M);
204
205 /* (1 - xy)/((1-x)(1-xy)) */
206 cnt = new infinite_counter(2, 3);
207 cnt->init(ctx);
208 set_from_string(n_coeff, "[1/1 -1/1]");
209 set_from_string(n_power, "[[0 0][1 1]]");
210 set_from_string(d_power, "[[1 0][1 1]]");
211 cnt->reduce(n_coeff, n_power, d_power);
212 assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) != 0);
213 assert(value_cmp_si(mpq_numref(cnt->count[2]), 0) == 0);
214 assert(value_cmp_si(mpq_denref(cnt->count[2]), 1) == 0);
215 assert(value_cmp_si(mpq_numref(cnt->count[3]), 0) == 0);
216 assert(value_cmp_si(mpq_denref(cnt->count[3]), 1) == 0);
217 delete cnt;
218
219 cnt = new infinite_counter(2, 2);
220 cnt->init(ctx);
221 set_from_string(n_coeff, "[-1/2 1/1 -1/3]");
222 set_from_string(n_power, "[[2 6][3 6]]");
223 d_power.SetDims(0, 2);
224 cnt->reduce(n_coeff, n_power, d_power);
225 assert(value_cmp_si(mpq_numref(cnt->count[0]), 1) == 0);
226 assert(value_cmp_si(mpq_denref(cnt->count[0]), 6) == 0);
227 assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) == 0);
228 assert(value_cmp_si(mpq_denref(cnt->count[1]), 1) == 0);
229 assert(value_cmp_si(mpq_numref(cnt->count[2]), 0) == 0);
230 assert(value_cmp_si(mpq_denref(cnt->count[2]), 1) == 0);
231 delete cnt;
232
233 cnt = new infinite_counter(2, 2);
234 cnt->init(ctx);
235 set_from_string(n_coeff, "[1/1]");
236 set_from_string(n_power, "[[0 11]]");
237 set_from_string(d_power, "[[0 1]]");
238 cnt->reduce(n_coeff, n_power, d_power);
239 assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) != 0);
240 assert(value_cmp_si(mpq_numref(cnt->count[2]), 0) == 0);
241 assert(value_cmp_si(mpq_denref(cnt->count[2]), 1) == 0);
242 delete cnt;
243
244 Polyhedron_Free(ctx);
245
246 return 0;
247 }
248
249 static int test_series(struct barvinok_options *options)
250 {
251 Matrix *M = matrix_read_from_str(
252 "12 11\n"
253 " 0 1 0 0 0 0 0 1 0 0 3 \n"
254 " 0 0 1 0 0 0 0 -1 1 0 -5 \n"
255 " 0 0 0 1 0 0 0 0 -2 -1 6 \n"
256 " 0 0 0 0 1 0 0 1 1 0 5 \n"
257 " 0 0 0 0 0 1 0 0 -1 0 0 \n"
258 " 0 0 0 0 0 0 1 -2 0 -1 -3 \n"
259 " 1 0 0 0 0 0 0 2 0 1 3 \n"
260 " 1 0 0 0 0 0 0 1 -1 0 5 \n"
261 " 1 0 0 0 0 0 0 -1 -1 0 -5 \n"
262 " 1 0 0 0 0 0 0 -1 0 0 -3 \n"
263 " 1 0 0 0 0 0 0 0 2 1 -6 \n"
264 " 1 0 0 0 0 0 0 0 1 0 0 \n");
265 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
266 Matrix_Free(M);
267 Polyhedron *C = Universe_Polyhedron(3);
268 gen_fun *gf = barvinok_series_with_options(P, C, options);
269 Polyhedron_Free(P);
270 Polyhedron_Free(C);
271 delete gf;
272
273 M = matrix_read_from_str(
274 "7 8\n"
275 " 0 1 1 0 0 1 0 2 \n"
276 " 0 0 0 1 0 -2 0 6 \n"
277 " 0 0 0 0 1 -1 0 -1 \n"
278 " 0 0 0 0 0 0 1 0 \n"
279 " 1 0 1 0 0 0 0 0 \n"
280 " 1 0 -1 0 0 -1 0 -2 \n"
281 " 1 0 0 0 0 1 0 -3 \n");
282 P = Constraints2Polyhedron(M, options->MaxRays);
283 Matrix_Free(M);
284 C = Universe_Polyhedron(2);
285 gf = barvinok_series_with_options(P, C, options);
286 Polyhedron_Free(P);
287 Polyhedron_Free(C);
288 delete gf;
289
290 M = matrix_read_from_str(
291 "2 3\n"
292 "1 1 0\n"
293 "1 -1 10\n");
294 P = Constraints2Polyhedron(M, options->MaxRays);
295 Matrix_Free(M);
296 C = Universe_Polyhedron(1);
297 gf = barvinok_series_with_options(P, C, options);
298 Polyhedron_Free(P);
299 Polyhedron_Free(C);
300 gen_fun *sum = gf->summate(1, options);
301 delete gf;
302 delete sum;
303
304 return 0;
305 }
306
307 int test_todd(struct barvinok_options *options)
308 {
309 tcounter t(2, options->max_index);
310 assert(value_cmp_si(t.todd.coeff->p[0], 1) == 0);
311 assert(value_cmp_si(t.todd.coeff->p[1], -3) == 0);
312 assert(value_cmp_si(t.todd.coeff->p[2], 3) == 0);
313 assert(value_cmp_si(t.todd_denom->p[0], 1) == 0);
314 assert(value_cmp_si(t.todd_denom->p[1], 6) == 0);
315 assert(value_cmp_si(t.todd_denom->p[2], 36) == 0);
316
317 vec_ZZ lambda;
318 set_from_string(lambda, "[1 -1]");
319 zz2values(lambda, t.lambda->p);
320
321 mat_ZZ rays;
322 set_from_string(rays, "[[-1 0][-1 1]]");
323
324 QQ one(1, 1);
325
326 vec_ZZ v;
327 set_from_string(v, "[2 0 1]");
328 Vector *vertex = Vector_Alloc(3);
329 zz2values(v, vertex->p);
330
331 t.handle(rays, vertex->p, one, 1, options);
332 assert(value_cmp_si(mpq_numref(t.count), 71) == 0);
333 assert(value_cmp_si(mpq_denref(t.count), 24) == 0);
334
335 set_from_string(rays, "[[0 -1][1 -1]]");
336 set_from_string(v, "[0 2 1]");
337 zz2values(v, vertex->p);
338
339 t.handle(rays, vertex->p, one, 1, options);
340 assert(value_cmp_si(mpq_numref(t.count), 71) == 0);
341 assert(value_cmp_si(mpq_denref(t.count), 12) == 0);
342
343 set_from_string(rays, "[[1 0][0 1]]");
344 set_from_string(v, "[0 0 1]");
345 zz2values(v, vertex->p);
346
347 t.handle(rays, vertex->p, one, 1, options);
348 assert(value_cmp_si(mpq_numref(t.count), 6) == 0);
349 assert(value_cmp_si(mpq_denref(t.count), 1) == 0);
350
351 Vector_Free(vertex);
352 }
353
354 int test_bernoulli(struct barvinok_options *options)
355 {
356 struct bernoulli_coef *bernoulli_coef;
357 struct poly_list *faulhaber, *bernoulli;
358 bernoulli_coef = bernoulli_coef_compute(2);
359 faulhaber = faulhaber_compute(4);
360 bernoulli_coef = bernoulli_coef_compute(8);
361 assert(value_cmp_si(bernoulli_coef->num->p[6], 1) == 0);
362 assert(value_cmp_si(bernoulli_coef->den->p[6], 42) == 0);
363 assert(value_cmp_si(faulhaber->poly[3]->p[0], 0) == 0);
364 assert(value_cmp_si(faulhaber->poly[3]->p[1], 0) == 0);
365 assert(value_cmp_si(faulhaber->poly[3]->p[2], 1) == 0);
366 assert(value_cmp_si(faulhaber->poly[3]->p[3], -2) == 0);
367 assert(value_cmp_si(faulhaber->poly[3]->p[4], 1) == 0);
368 assert(value_cmp_si(faulhaber->poly[3]->p[5], 4) == 0);
369
370 bernoulli = bernoulli_compute(6);
371 assert(value_cmp_si(bernoulli->poly[6]->p[0], 1) == 0);
372 assert(value_cmp_si(bernoulli->poly[6]->p[1], 0) == 0);
373 assert(value_cmp_si(bernoulli->poly[6]->p[2], -21) == 0);
374 assert(value_cmp_si(bernoulli->poly[6]->p[3], 0) == 0);
375 assert(value_cmp_si(bernoulli->poly[6]->p[4], 105) == 0);
376 assert(value_cmp_si(bernoulli->poly[6]->p[5], -126) == 0);
377 assert(value_cmp_si(bernoulli->poly[6]->p[6], 42) == 0);
378 assert(value_cmp_si(bernoulli->poly[6]->p[7], 42) == 0);
379
380 unsigned nvar, nparam;
381 char **all_vars;
382 evalue *base, *sum1, *sum2;
383 base = evalue_read_from_str("(1 * n + 1)", NULL, &all_vars, &nvar, &nparam,
384 options->MaxRays);
385
386 sum1 = evalue_polynomial(faulhaber->poly[3], base);
387 Free_ParamNames(all_vars, nvar+nparam);
388
389 sum2 = evalue_read_from_str("(1/4 * n^4 + 1/2 * n^3 + 1/4 * n^2)",
390 NULL, &all_vars, &nvar, &nparam,
391 options->MaxRays);
392 Free_ParamNames(all_vars, nvar+nparam);
393 assert(eequal(sum1, sum2));
394 evalue_free(base);
395 evalue_free(sum1);
396 evalue_free(sum2);
397 }
398
399 int test_bernoulli_sum(struct barvinok_options *options)
400 {
401 unsigned nvar, nparam;
402 char **all_vars;
403 evalue *e, *sum1, *sum2;
404 e = evalue_read_from_str("i + -1 >= 0\n -i + n >= 0\n\n 1 + (-1 *i) + i^2",
405 "i", &all_vars, &nvar, &nparam,
406 options->MaxRays);
407 Free_ParamNames(all_vars, nvar+nparam);
408
409 sum1 = Bernoulli_sum_evalue(e, 1, options);
410 sum2 = evalue_read_from_str("n -1 >= 0\n\n (1/3 * n^3 + 2/3 * n)",
411 NULL, &all_vars, &nvar, &nparam,
412 options->MaxRays);
413 Free_ParamNames(all_vars, nvar+nparam);
414 evalue_negate(sum1);
415 eadd(sum2, sum1);
416 reduce_evalue(sum1);
417 assert(EVALUE_IS_ZERO(*sum1));
418 evalue_free(e);
419 evalue_free(sum1);
420
421 e = evalue_read_from_str("-i + -1 >= 0\n i + n >= 0\n\n 1 + i + i^2",
422 "i", &all_vars, &nvar, &nparam,
423 options->MaxRays);
424 Free_ParamNames(all_vars, nvar+nparam);
425 sum1 = Bernoulli_sum_evalue(e, 1, options);
426 evalue_negate(sum1);
427 eadd(sum2, sum1);
428 reduce_evalue(sum1);
429 assert(EVALUE_IS_ZERO(*sum1));
430 evalue_free(e);
431 evalue_free(sum1);
432
433 evalue_free(sum2);
434
435 e = evalue_read_from_str("i + 4 >= 0\n -i + n >= 0\n\n i",
436 "i", &all_vars, &nvar, &nparam,
437 options->MaxRays);
438 Free_ParamNames(all_vars, nvar+nparam);
439 sum1 = Bernoulli_sum_evalue(e, 1, options);
440 sum2 = evalue_read_from_str("n + 0 >= 0\n\n (1/2 * n^2 + 1/2 * n + (-10))\n"
441 "n + 4 >= 0\n -n -1 >= 0\n\n (1/2 * n^2 + 1/2 * n + (-10))",
442 NULL, &all_vars, &nvar, &nparam,
443 options->MaxRays);
444 Free_ParamNames(all_vars, nvar+nparam);
445 evalue_negate(sum1);
446 eadd(sum2, sum1);
447 reduce_evalue(sum1);
448 assert(EVALUE_IS_ZERO(*sum1));
449 evalue_free(e);
450 evalue_free(sum1);
451 evalue_free(sum2);
452
453 e = evalue_read_from_str("i -5 >= 0\n -i + n >= 0\n j -1 >= 0\n i -j >= 0\n"
454 "k -1 >= 0\n j -k >= 0\n\n1",
455 "i,j,k", &all_vars, &nvar, &nparam,
456 options->MaxRays);
457 Free_ParamNames(all_vars, nvar+nparam);
458 sum1 = Bernoulli_sum_evalue(e, 3, options);
459 sum2 = evalue_read_from_str("n -5 >= 0\n\n"
460 "1/6 * n^3 + 1/2 * n^2 + 1/3 * n + -20",
461 NULL, &all_vars, &nvar, &nparam,
462 options->MaxRays);
463 Free_ParamNames(all_vars, nvar+nparam);
464 evalue_negate(sum1);
465 eadd(sum2, sum1);
466 reduce_evalue(sum1);
467 assert(EVALUE_IS_ZERO(*sum1));
468 evalue_free(e);
469 evalue_free(sum1);
470 evalue_free(sum2);
471 }
472
473 int test_hilbert(struct barvinok_options *options)
474 {
475 Matrix *M = matrix_read_from_str(
476 "2 4\n"
477 " 1 4 -3 0 \n"
478 " 1 3 2 0 \n");
479 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
480 Matrix_Free(M);
481
482 M = Cone_Hilbert_Basis(P, options->MaxRays);
483 assert(M->NbRows = 5);
484 assert(M->NbColumns = 3);
485 Matrix_Free(M);
486
487 M = Cone_Integer_Hull(P, options);
488 assert(M->NbRows = 4);
489 assert(M->NbColumns = 3);
490 Matrix_Free(M);
491
492 Polyhedron_Free(P);
493 }
494
495 int main(int argc, char **argv)
496 {
497 struct barvinok_options *options = barvinok_options_new_with_defaults();
498 test_evalue(options);
499 test_split_periods(options);
500 test_specialization(options);
501 test_lattice_points(options);
502 test_icounter(options);
503 test_infinite_counter(options);
504 test_series(options);
505 test_todd(options);
506 test_bernoulli(options);
507 test_bernoulli_sum(options);
508 test_hilbert(options);
509 barvinok_options_free(options);
510 }
lattice point enumeration library