testlib.cc

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