testlib.cc

   1 #include <assert.h>
   2 #include <sstream>
   3 #include <barvinok/set.h>
   4 #include <barvinok/options.h>
   5 #include <barvinok/evalue.h>
   6 #include <barvinok/util.h>
   7 #include "conversion.h"
   8 #include "evalue_read.h"
   9 #include "dpoly.h"
  10 #include "lattice_point.h"
  11 #include "tcounter.h"
  12 #include "bernoulli.h"
  13
  14 using std::cout;
  15 using std::cerr;
  16 using std::endl;
  17
  18 template <typename T>
  19 void set_from_string(T& v, const char *s)
  20 {
  21     std::istringstream str(s);
  22     str >> v;
  23 }
  24
  25 int test_evalue(struct barvinok_options *options)
  26 {
  27     unsigned nvar, nparam;
  28     char **all_vars;
  29     evalue *poly1, poly2;
  30
  31     poly1 = evalue_read_from_str("(1/4 * n^4 + 1/2 * n^3 + 1/4 * n^2)",
  32                                  NULL, &all_vars, &nvar, &nparam,
  33                                  options->MaxRays);
  34     Free_ParamNames(all_vars, nvar+nparam);
  35
  36     value_init(poly2.d);
  37     evalue_copy(&poly2, poly1);
  38     evalue_negate(poly1);
  39     eadd(&poly2, poly1);
  40     reduce_evalue(poly1);
  41     assert(EVALUE_IS_ZERO(*poly1));
  42     evalue_free(poly1);
  43     free_evalue_refs(&poly2);
  44 }
  45
  46 int test_split_periods(struct barvinok_options *options)
  47 {
  48     unsigned nvar, nparam;
  49     char **all_vars;
  50     evalue *e;
  51
  52     e = evalue_read_from_str("U + 2V + 3 >= 0\n- U -2V  >= 0\n- U  10 >= 0\n"
  53                              "U  >= 0\n\n({( 1/3 * U + ( 2/3 * V + 0 ))})",
  54                              NULL, &all_vars, &nvar, &nparam,
  55                              options->MaxRays);
  56     Free_ParamNames(all_vars, nvar+nparam);
  57
  58     evalue_split_periods(e, 2, options->MaxRays);
  59     assert(value_zero_p(e->d));
  60     assert(e->x.p->type == partition);
  61     assert(e->x.p->size == 4);
  62     assert(value_zero_p(e->x.p->arr[1].d));
  63     assert(e->x.p->arr[1].x.p->type == polynomial);
  64     assert(value_zero_p(e->x.p->arr[3].d));
  65     assert(e->x.p->arr[3].x.p->type == polynomial);
  66     evalue_free(e);
  67 }
  68
  69 int test_specialization(struct barvinok_options *options)
  70 {
  71     Value v;
  72     value_init(v);
  73     mpq_t count;
  74     mpq_init(count);
  75
  76     value_set_si(v, 5);
  77     dpoly n(2, v);
  78     assert(value_cmp_si(n.coeff->p[0], 1) == 0);
  79     assert(value_cmp_si(n.coeff->p[1], 5) == 0);
  80     assert(value_cmp_si(n.coeff->p[2], 10) == 0);
  81
  82     value_set_si(v, 1);
  83     dpoly d(2, v, 1);
  84     value_set_si(v, 2);
  85     dpoly d2(2, v, 1);
  86     d *= d2;
  87     assert(value_cmp_si(d.coeff->p[0], 2) == 0);
  88     assert(value_cmp_si(d.coeff->p[1], 1) == 0);
  89     assert(value_cmp_si(d.coeff->p[2], 0) == 0);
  90
  91     n.div(d, count, 1);
  92     mpq_canonicalize(count);
  93     assert(value_cmp_si(mpq_numref(count), 31) == 0);
  94     assert(value_cmp_si(mpq_denref(count), 8) == 0);
  95
  96     value_set_si(v, -2);
  97     dpoly n2(2, v);
  98     assert(value_cmp_si(n2.coeff->p[0], 1) == 0);
  99     assert(value_cmp_si(n2.coeff->p[1], -2) == 0);
 100     assert(value_cmp_si(n2.coeff->p[2], 3) == 0);
 101
 102     n2.div(d, count, 1);
 103     mpq_canonicalize(count);
 104     assert(value_cmp_si(mpq_numref(count), 6) == 0);
 105     assert(value_cmp_si(mpq_denref(count), 1) == 0);
 106
 107     value_clear(v);
 108     mpq_clear(count);
 109 }
 110
 111 int test_lattice_points(struct barvinok_options *options)
 112 {
 113     Param_Vertices V;
 114     mat_ZZ tmp;
 115     set_from_string(tmp, "[[0 0 0 0 4][0 0 0 0 4][-1 0 1 0 4]]");
 116     V.Vertex = zz2matrix(tmp);
 117     vec_ZZ lambda;
 118     set_from_string(lambda, "[3 5 7]");
 119     mat_ZZ rays;
 120     set_from_string(rays, "[[0 1 0][4 0 1][0 0 -1]]");
 121     term_info num;
 122     evalue *point[4];
 123
 124     unsigned nvar, nparam;
 125     char **all_vars;
 126     point[0] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
 127                     "( 7 * {( 1/4 * a + ( 3/4 * c + 3/4 ) ) } + -21/4 ) ) )",
 128                     "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
 129     Free_ParamNames(all_vars, nvar+nparam);
 130     point[1] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
 131                     "( 7 * {( 1/4 * a + ( 3/4 * c + 1/2 ) ) } + -1/2 ) ) )",
 132                     "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
 133     Free_ParamNames(all_vars, nvar+nparam);
 134     point[2] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
 135                     "( 7 * {( 1/4 * a + ( 3/4 * c + 1/4 ) ) } + 17/4 ) ) )",
 136                     "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
 137     Free_ParamNames(all_vars, nvar+nparam);
 138     point[3] = evalue_read_from_str("( -7/4 * a + ( 7/4 * c + "
 139                     "( 7 * {( 1/4 * a + ( 3/4 * c + 0 ) ) } + 9 ) ) )",
 140                     "a,b,c", &all_vars, &nvar, &nparam, options->MaxRays);
 141     Free_ParamNames(all_vars, nvar+nparam);
 142
 143     lattice_point(&V, rays, lambda, &num, 4, options);
 144     Matrix_Free(V.Vertex);
 145
 146     for (int i = 0; i < 4; ++i) {
 147         assert(eequal(num.E[i], point[i]));
 148         evalue_free(point[i]);
 149         free_evalue_refs(num.E[i]);
 150         delete num.E[i];
 151     }
 152     delete [] num.E;
 153 }
 154
 155 int test_todd(struct barvinok_options *options)
 156 {
 157     tcounter t(2, options->max_index);
 158     assert(value_cmp_si(t.todd.coeff->p[0], 1) == 0);
 159     assert(value_cmp_si(t.todd.coeff->p[1], -3) == 0);
 160     assert(value_cmp_si(t.todd.coeff->p[2], 3) == 0);
 161     assert(value_cmp_si(t.todd_denom->p[0], 1) == 0);
 162     assert(value_cmp_si(t.todd_denom->p[1], 6) == 0);
 163     assert(value_cmp_si(t.todd_denom->p[2], 36) == 0);
 164
 165     vec_ZZ lambda;
 166     set_from_string(lambda, "[1 -1]");
 167     zz2values(lambda, t.lambda->p);
 168
 169     mat_ZZ rays;
 170     set_from_string(rays, "[[-1 0][-1 1]]");
 171
 172     QQ one(1, 1);
 173
 174     vec_ZZ v;
 175     set_from_string(v, "[2 0 1]");
 176     Vector *vertex = Vector_Alloc(3);
 177     zz2values(v, vertex->p);
 178
 179     t.handle(rays, vertex->p, one, 1, options);
 180     assert(value_cmp_si(mpq_numref(t.count), 71) == 0);
 181     assert(value_cmp_si(mpq_denref(t.count), 24) == 0);
 182
 183     set_from_string(rays, "[[0 -1][1 -1]]");
 184     set_from_string(v, "[0 2 1]");
 185     zz2values(v, vertex->p);
 186
 187     t.handle(rays, vertex->p, one, 1, options);
 188     assert(value_cmp_si(mpq_numref(t.count), 71) == 0);
 189     assert(value_cmp_si(mpq_denref(t.count), 12) == 0);
 190
 191     set_from_string(rays, "[[1 0][0 1]]");
 192     set_from_string(v, "[0 0 1]");
 193     zz2values(v, vertex->p);
 194
 195     t.handle(rays, vertex->p, one, 1, options);
 196     assert(value_cmp_si(mpq_numref(t.count), 6) == 0);
 197     assert(value_cmp_si(mpq_denref(t.count), 1) == 0);
 198
 199     Vector_Free(vertex);
 200 }
 201
 202 int test_bernoulli(struct barvinok_options *options)
 203 {
 204     struct bernoulli_coef *bernoulli_coef;
 205     struct poly_list *faulhaber, *bernoulli;
 206     bernoulli_coef = bernoulli_coef_compute(2);
 207     faulhaber = faulhaber_compute(4);
 208     bernoulli_coef = bernoulli_coef_compute(8);
 209     assert(value_cmp_si(bernoulli_coef->num->p[6], 1) == 0);
 210     assert(value_cmp_si(bernoulli_coef->den->p[6], 42) == 0);
 211     assert(value_cmp_si(faulhaber->poly[3]->p[0], 0) == 0);
 212     assert(value_cmp_si(faulhaber->poly[3]->p[1], 0) == 0);
 213     assert(value_cmp_si(faulhaber->poly[3]->p[2], 1) == 0);
 214     assert(value_cmp_si(faulhaber->poly[3]->p[3], -2) == 0);
 215     assert(value_cmp_si(faulhaber->poly[3]->p[4], 1) == 0);
 216     assert(value_cmp_si(faulhaber->poly[3]->p[5], 4) == 0);
 217
 218     bernoulli = bernoulli_compute(6);
 219     assert(value_cmp_si(bernoulli->poly[6]->p[0], 1) == 0);
 220     assert(value_cmp_si(bernoulli->poly[6]->p[1], 0) == 0);
 221     assert(value_cmp_si(bernoulli->poly[6]->p[2], -21) == 0);
 222     assert(value_cmp_si(bernoulli->poly[6]->p[3], 0) == 0);
 223     assert(value_cmp_si(bernoulli->poly[6]->p[4], 105) == 0);
 224     assert(value_cmp_si(bernoulli->poly[6]->p[5], -126) == 0);
 225     assert(value_cmp_si(bernoulli->poly[6]->p[6], 42) == 0);
 226     assert(value_cmp_si(bernoulli->poly[6]->p[7], 42) == 0);
 227
 228     unsigned nvar, nparam;
 229     char **all_vars;
 230     evalue *base, *sum1, *sum2;
 231     base = evalue_read_from_str("(1 * n + 1)", NULL, &all_vars, &nvar, &nparam,
 232                                 options->MaxRays);
 233
 234     sum1 = evalue_polynomial(faulhaber->poly[3], base);
 235     Free_ParamNames(all_vars, nvar+nparam);
 236
 237     sum2 = evalue_read_from_str("(1/4 * n^4 + 1/2 * n^3 + 1/4 * n^2)",
 238                                 NULL, &all_vars, &nvar, &nparam,
 239                                 options->MaxRays);
 240     Free_ParamNames(all_vars, nvar+nparam);
 241     assert(eequal(sum1, sum2));
 242     evalue_free(base);
 243     evalue_free(sum1);
 244     evalue_free(sum2);
 245 }
 246
 247 int test_bernoulli_sum(struct barvinok_options *options)
 248 {
 249     unsigned nvar, nparam;
 250     char **all_vars;
 251     evalue *e, *sum1, *sum2;
 252     e = evalue_read_from_str("i + -1 >= 0\n -i + n >= 0\n\n 1 + (-1 *i) + i^2",
 253                              "i", &all_vars, &nvar, &nparam,
 254                                 options->MaxRays);
 255     Free_ParamNames(all_vars, nvar+nparam);
 256
 257     sum1 = Bernoulli_sum_evalue(e, 1, options);
 258     sum2 = evalue_read_from_str("n -1 >= 0\n\n (1/3 * n^3 + 2/3 * n)",
 259                                 NULL, &all_vars, &nvar, &nparam,
 260                                 options->MaxRays);
 261     Free_ParamNames(all_vars, nvar+nparam);
 262     evalue_negate(sum1);
 263     eadd(sum2, sum1);
 264     reduce_evalue(sum1);
 265     assert(EVALUE_IS_ZERO(*sum1));
 266     evalue_free(e);
 267     evalue_free(sum1);
 268
 269     e = evalue_read_from_str("-i + -1 >= 0\n i + n >= 0\n\n 1 + i + i^2",
 270                              "i", &all_vars, &nvar, &nparam,
 271                                 options->MaxRays);
 272     Free_ParamNames(all_vars, nvar+nparam);
 273     sum1 = Bernoulli_sum_evalue(e, 1, options);
 274     evalue_negate(sum1);
 275     eadd(sum2, sum1);
 276     reduce_evalue(sum1);
 277     assert(EVALUE_IS_ZERO(*sum1));
 278     evalue_free(e);
 279     evalue_free(sum1);
 280
 281     evalue_free(sum2);
 282
 283     e = evalue_read_from_str("i + 4 >= 0\n -i + n >= 0\n\n i",
 284                              "i", &all_vars, &nvar, &nparam,
 285                                 options->MaxRays);
 286     Free_ParamNames(all_vars, nvar+nparam);
 287     sum1 = Bernoulli_sum_evalue(e, 1, options);
 288     sum2 = evalue_read_from_str("n + 0 >= 0\n\n (1/2 * n^2 + 1/2 * n + (-10))\n"
 289                     "n + 4 >= 0\n -n -1 >= 0\n\n (1/2 * n^2 + 1/2 * n + (-10))",
 290                                 NULL, &all_vars, &nvar, &nparam,
 291                                 options->MaxRays);
 292     Free_ParamNames(all_vars, nvar+nparam);
 293     evalue_negate(sum1);
 294     eadd(sum2, sum1);
 295     reduce_evalue(sum1);
 296     assert(EVALUE_IS_ZERO(*sum1));
 297     evalue_free(e);
 298     evalue_free(sum1);
 299     evalue_free(sum2);
 300
 301     e = evalue_read_from_str("i -5 >= 0\n -i + n >= 0\n j -1 >= 0\n i -j >= 0\n"
 302                              "k -1 >= 0\n j -k >= 0\n\n1",
 303                              "i,j,k", &all_vars, &nvar, &nparam,
 304                                 options->MaxRays);
 305     Free_ParamNames(all_vars, nvar+nparam);
 306     sum1 = Bernoulli_sum_evalue(e, 3, options);
 307     sum2 = evalue_read_from_str("n -5 >= 0\n\n"
 308                                 "1/6 * n^3 + 1/2 * n^2 + 1/3 * n + -20",
 309                                 NULL, &all_vars, &nvar, &nparam,
 310                                 options->MaxRays);
 311     Free_ParamNames(all_vars, nvar+nparam);
 312     evalue_negate(sum1);
 313     eadd(sum2, sum1);
 314     reduce_evalue(sum1);
 315     assert(EVALUE_IS_ZERO(*sum1));
 316     evalue_free(e);
 317     evalue_free(sum1);
 318     evalue_free(sum2);
 319 }
 320
 321 int main(int argc, char **argv)
 322 {
 323     struct barvinok_options *options = barvinok_options_new_with_defaults();
 324     test_evalue(options);
 325     test_split_periods(options);
 326     test_specialization(options);
 327     test_lattice_points(options);
 328     test_todd(options);
 329     test_bernoulli(options);
 330     test_bernoulli_sum(options);
 331     barvinok_options_free(options);
 332 }