summate.c

   1 #include <barvinok/options.h>
   2 #include <barvinok/util.h>
   3 #include "bernoulli.h"
   4 #include "euler.h"
   5 #include "laurent.h"
   6 #include "summate.h"
   7 #include "section_array.h"
   8
   9 extern evalue *evalue_outer_floor(evalue *e);
  10 extern int evalue_replace_floor(evalue *e, const evalue *floor, int var);
  11 extern void evalue_drop_floor(evalue *e, const evalue *floor);
  12
  13 #define ALLOC(type) (type*)malloc(sizeof(type))
  14 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  15
  16 /* Apply the variable transformation specified by T and CP on
  17  * the polynomial e.  T expresses the old variables in terms
  18  * of the new variables (and optionally also the new parameters),
  19  * while CP expresses the old parameters in terms of the new
  20  * parameters.
  21  */
  22 static void transform_polynomial(evalue *E, Matrix *T, Matrix *CP,
  23                                  unsigned nvar, unsigned nparam,
  24                                  unsigned new_nvar, unsigned new_nparam)
  25 {
  26     int j;
  27     evalue **subs;
  28
  29     subs = ALLOCN(evalue *, nvar+nparam);
  30
  31     for (j = 0; j < nvar; ++j) {
  32         if (T)
  33             subs[j] = affine2evalue(T->p[j], T->p[T->NbRows-1][T->NbColumns-1],
  34                                     T->NbColumns-1);
  35         else
  36             subs[j] = evalue_var(j);
  37     }
  38     for (j = 0; j < nparam; ++j) {
  39         if (CP)
  40             subs[nvar+j] = affine2evalue(CP->p[j], CP->p[nparam][new_nparam],
  41                                          new_nparam);
  42         else
  43             subs[nvar+j] = evalue_var(j);
  44         evalue_shift_variables(subs[nvar+j], 0, new_nvar);
  45     }
  46
  47     evalue_substitute(E, subs);
  48     reduce_evalue(E);
  49
  50     for (j = 0; j < nvar+nparam; ++j)
  51         evalue_free(subs[j]);
  52     free(subs);
  53 }
  54
  55 static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E,
  56                                  unsigned nvar,
  57                                  struct evalue_section_array *sections,
  58                                  struct barvinok_options *options)
  59 {
  60     unsigned dim = P->Dimension;
  61     unsigned new_dim, new_nparam;
  62     Matrix *T = NULL, *CP = NULL;
  63     evalue *sum;
  64
  65     if (emptyQ(P))
  66         return evalue_zero();
  67
  68     assert(P->NbEq > 0);
  69
  70     remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays);
  71
  72     if (emptyQ(P)) {
  73         Polyhedron_Free(P);
  74         return evalue_zero();
  75     }
  76
  77     new_nparam = CP ? CP->NbColumns-1 : dim - nvar;
  78     new_dim = T ? T->NbColumns-1 : nvar + new_nparam;
  79
  80     /* We can avoid these substitutions if E is a constant */
  81     E = evalue_dup(E);
  82     transform_polynomial(E, T, CP, nvar, dim-nvar,
  83                          new_dim-new_nparam, new_nparam);
  84
  85     if (new_dim-new_nparam > 0) {
  86         sum = barvinok_sum_over_polytope(P, E, new_dim-new_nparam,
  87                                          sections, options);
  88         evalue_free(E);
  89         Polyhedron_Free(P);
  90     } else {
  91         sum = ALLOC(evalue);
  92         value_init(sum->d);
  93         sum->x.p = new_enode(partition, 2, new_dim);
  94         EVALUE_SET_DOMAIN(sum->x.p->arr[0], P);
  95         value_clear(sum->x.p->arr[1].d);
  96         sum->x.p->arr[1] = *E;
  97         free(E);
  98     }
  99
 100     if (CP) {
 101         evalue_backsubstitute(sum, CP, options->MaxRays);
 102         Matrix_Free(CP);
 103     }
 104
 105     if (T)
 106         Matrix_Free(T);
 107
 108     return sum;
 109 }
 110
 111 static evalue *sum_base(Polyhedron *P, evalue *E, unsigned nvar,
 112                                      struct barvinok_options *options)
 113 {
 114     if (options->summation == BV_SUM_EULER)
 115         return euler_summate(P, E, nvar, options);
 116     else if (options->summation == BV_SUM_LAURENT)
 117         return laurent_summate(P, E, nvar, options);
 118     assert(0);
 119 }
 120
 121 /* Count the number of non-zero terms in e when viewed as a polynomial
 122  * in only the first nvar variables.  "count" is the number counted
 123  * so far.
 124  */
 125 static int evalue_count_terms(const evalue *e, unsigned nvar, int count)
 126 {
 127     int i;
 128
 129     if (EVALUE_IS_ZERO(*e))
 130         return count;
 131
 132     if (value_zero_p(e->d))
 133         assert(e->x.p->type == polynomial);
 134     if (value_notzero_p(e->d) || e->x.p->pos >= nvar+1)
 135         return count+1;
 136
 137     for (i = 0; i < e->x.p->size; ++i)
 138         count = evalue_count_terms(&e->x.p->arr[i], nvar, count);
 139
 140     return count;
 141 }
 142
 143 /* Create placeholder structure for unzipping.
 144  * A "polynomial" is created with size terms in variable pos,
 145  * with each term having itself as coefficient.
 146  */
 147 static evalue *create_placeholder(int size, int pos)
 148 {
 149     int i, j;
 150     evalue *E = ALLOC(evalue);
 151     value_init(E->d);
 152     E->x.p = new_enode(polynomial, size, pos+1);
 153     for (i = 0; i < size; ++i) {
 154         E->x.p->arr[i].x.p = new_enode(polynomial, i+1, pos+1);
 155         for (j = 0; j < i; ++j)
 156             evalue_set_si(&E->x.p->arr[i].x.p->arr[j], 0, 1);
 157         evalue_set_si(&E->x.p->arr[i].x.p->arr[i], 1, 1);
 158     }
 159     return E;
 160 }
 161
 162 /* Interchange each non-zero term in e (when viewed as a polynomial
 163  * in only the first nvar variables) with a placeholder in ph (created
 164  * by create_placeholder), resulting in two polynomials in the
 165  * placeholder variable such that for each non-zero term in e
 166  * there is a power of the placeholder variable such that the factors
 167  * in the first nvar variables form the coefficient of that power in
 168  * the first polynomial (e) and the factors in the remaining variables
 169  * form the coefficient of that power in the second polynomial (ph).
 170  */
 171 static int evalue_unzip_terms(evalue *e, evalue *ph, unsigned nvar, int count)
 172 {
 173     int i;
 174
 175     if (EVALUE_IS_ZERO(*e))
 176         return count;
 177
 178     if (value_zero_p(e->d))
 179         assert(e->x.p->type == polynomial);
 180     if (value_notzero_p(e->d) || e->x.p->pos >= nvar+1) {
 181         evalue t = *e;
 182         *e = ph->x.p->arr[count];
 183         ph->x.p->arr[count] = t;
 184         return count+1;
 185     }
 186
 187     for (i = 0; i < e->x.p->size; ++i)
 188         count = evalue_unzip_terms(&e->x.p->arr[i], ph, nvar, count);
 189
 190     return count;
 191 }
 192
 193 /* Remove n variables at pos (0-based) from the polyhedron P.
 194  * Each of these variables is assumed to be completely free,
 195  * i.e., there is a line in the polyhedron corresponding to
 196  * each of these variables.
 197  */
 198 static Polyhedron *Polyhedron_Remove_Columns(Polyhedron *P, unsigned pos,
 199                                         unsigned n)
 200 {
 201     int i, j;
 202     unsigned NbConstraints = 0;
 203     unsigned NbRays = 0;
 204     Polyhedron *Q;
 205
 206     if (!P)
 207         return;
 208     if (n == 0)
 209         return;
 210
 211     assert(pos <= P->Dimension);
 212
 213     if (POL_HAS(P, POL_INEQUALITIES))
 214         NbConstraints = P->NbConstraints;
 215     if (POL_HAS(P, POL_POINTS))
 216         NbRays = P->NbRays - n;
 217
 218     Q = Polyhedron_Alloc(P->Dimension - n, NbConstraints, NbRays);
 219     if (POL_HAS(P, POL_INEQUALITIES)) {
 220         Q->NbEq = P->NbEq;
 221         for (i = 0; i < P->NbConstraints; ++i) {
 222             Vector_Copy(P->Constraint[i], Q->Constraint[i], 1+pos);
 223             Vector_Copy(P->Constraint[i]+1+pos+n, Q->Constraint[i]+1+pos,
 224                         Q->Dimension-pos+1);
 225         }
 226     }
 227     if (POL_HAS(P, POL_POINTS)) {
 228         Q->NbBid = P->NbBid - n;
 229         for (i = 0; i < n; ++i)
 230             value_set_si(Q->Ray[i][1+pos+i], 1);
 231         for (i = 0, j = 0; i < P->NbRays; ++i) {
 232             int line = First_Non_Zero(P->Ray[i], 1+P->Dimension+1);
 233             assert(line != -1);
 234             if (line-1 >= pos && line-1 < pos+n) {
 235                 ++j;
 236                 assert(j <= n);
 237                 continue;
 238             }
 239             assert(i-j < Q->NbRays);
 240             Vector_Copy(P->Ray[i], Q->Ray[i-j], 1+pos);
 241             Vector_Copy(P->Ray[i]+1+pos+n, Q->Ray[i-j]+1+pos,
 242                         Q->Dimension-pos+1);
 243         }
 244     }
 245     POL_SET(Q, POL_VALID);
 246     if (POL_HAS(P, POL_INEQUALITIES))
 247         POL_SET(Q, POL_INEQUALITIES);
 248     if (POL_HAS(P, POL_POINTS))
 249         POL_SET(Q, POL_POINTS);
 250     if (POL_HAS(P, POL_VERTICES))
 251         POL_SET(Q, POL_VERTICES);
 252     return Q;
 253 }
 254
 255 /* Remove n variables at pos (0-based) from the union of polyhedra P.
 256  * Each of these variables is assumed to be completely free,
 257  * i.e., there is a line in the polyhedron corresponding to
 258  * each of these variables.
 259  */
 260 static Polyhedron *Domain_Remove_Columns(Polyhedron *P, unsigned pos,
 261                                         unsigned n)
 262 {
 263     Polyhedron *R;
 264     Polyhedron **next = &R;
 265
 266     for (; P; P = P->next) {
 267         *next = Polyhedron_Remove_Columns(P, pos, n);
 268         next = &(*next)->next;
 269     }
 270     return R;
 271 }
 272
 273 /* Drop n parameters starting at first from partition evalue e */
 274 static void drop_parameters(evalue *e, int first, int n)
 275 {
 276     int i;
 277
 278     if (EVALUE_IS_ZERO(*e))
 279         return;
 280
 281     assert(value_zero_p(e->d) && e->x.p->type == partition);
 282     for (i = 0; i < e->x.p->size/2; ++i) {
 283         Polyhedron *P = EVALUE_DOMAIN(e->x.p->arr[2*i]);
 284         Polyhedron *Q = Domain_Remove_Columns(P, first, n);
 285         EVALUE_SET_DOMAIN(e->x.p->arr[2*i], Q);
 286         Domain_Free(P);
 287         evalue_shift_variables(&e->x.p->arr[2*i+1], first, -n);
 288     }
 289     e->x.p->pos -= n;
 290 }
 291
 292 static void extract_term_into(const evalue *src, int var, int exp, evalue *dst)
 293 {
 294     int i;
 295
 296     if (value_notzero_p(src->d) ||
 297             src->x.p->type != polynomial ||
 298             src->x.p->pos > var+1) {
 299         if (exp == 0)
 300             evalue_copy(dst, src);
 301         else
 302             evalue_set_si(dst, 0, 1);
 303         return;
 304     }
 305
 306     if (src->x.p->pos == var+1) {
 307         if (src->x.p->size > exp)
 308             evalue_copy(dst, &src->x.p->arr[exp]);
 309         else
 310             evalue_set_si(dst, 0, 1);
 311         return;
 312     }
 313
 314     dst->x.p = new_enode(polynomial, src->x.p->size, src->x.p->pos);
 315     for (i = 0; i < src->x.p->size; ++i)
 316         extract_term_into(&src->x.p->arr[i], var, exp,
 317                             &dst->x.p->arr[i]);
 318 }
 319
 320 /* Extract the coefficient of var^exp.
 321  */
 322 static evalue *extract_term(const evalue *e, int var, int exp)
 323 {
 324     int i;
 325     evalue *res;
 326
 327     if (EVALUE_IS_ZERO(*e))
 328         return evalue_zero();
 329
 330     assert(value_zero_p(e->d) && e->x.p->type == partition);
 331     res = ALLOC(evalue);
 332     value_init(res->d);
 333     res->x.p = new_enode(partition, e->x.p->size, e->x.p->pos);
 334     for (i = 0; i < e->x.p->size/2; ++i) {
 335         EVALUE_SET_DOMAIN(res->x.p->arr[2*i],
 336                           Domain_Copy(EVALUE_DOMAIN(e->x.p->arr[2*i])));
 337         extract_term_into(&e->x.p->arr[2*i+1], var, exp,
 338                             &res->x.p->arr[2*i+1]);
 339         reduce_evalue(&res->x.p->arr[2*i+1]);
 340     }
 341     return res;
 342 }
 343
 344 /* Insert n free variables at pos (0-based) in the polyhedron P.
 345  */
 346 static Polyhedron *Polyhedron_Insert_Columns(Polyhedron *P, unsigned pos,
 347                                         unsigned n)
 348 {
 349     int i;
 350     unsigned NbConstraints = 0;
 351     unsigned NbRays = 0;
 352     Polyhedron *Q;
 353
 354     if (!P)
 355         return;
 356     if (n == 0)
 357         return;
 358
 359     assert(pos <= P->Dimension);
 360
 361     if (POL_HAS(P, POL_INEQUALITIES))
 362         NbConstraints = P->NbConstraints;
 363     if (POL_HAS(P, POL_POINTS))
 364         NbRays = P->NbRays + n;
 365
 366     Q = Polyhedron_Alloc(P->Dimension+n, NbConstraints, NbRays);
 367     if (POL_HAS(P, POL_INEQUALITIES)) {
 368         Q->NbEq = P->NbEq;
 369         for (i = 0; i < P->NbConstraints; ++i) {
 370             Vector_Copy(P->Constraint[i], Q->Constraint[i], 1+pos);
 371             Vector_Copy(P->Constraint[i]+1+pos, Q->Constraint[i]+1+pos+n,
 372                         P->Dimension-pos+1);
 373         }
 374     }
 375     if (POL_HAS(P, POL_POINTS)) {
 376         Q->NbBid = P->NbBid + n;
 377         for (i = 0; i < n; ++i)
 378             value_set_si(Q->Ray[i][1+pos+i], 1);
 379         for (i = 0; i < P->NbRays; ++i) {
 380             Vector_Copy(P->Constraint[i], Q->Constraint[n+i], 1+pos);
 381             Vector_Copy(P->Constraint[i]+1+pos, Q->Constraint[n+i]+1+pos+n,
 382                         P->Dimension-pos+1);
 383         }
 384     }
 385     POL_SET(Q, POL_VALID);
 386     if (POL_HAS(P, POL_INEQUALITIES))
 387         POL_SET(Q, POL_INEQUALITIES);
 388     if (POL_HAS(P, POL_POINTS))
 389         POL_SET(Q, POL_POINTS);
 390     if (POL_HAS(P, POL_VERTICES))
 391         POL_SET(Q, POL_VERTICES);
 392     return Q;
 393 }
 394
 395 /* Perform summation of e over a list of 1 or more factors F, with context C.
 396  * nvar is the total number of variables in the remaining factors.
 397  * extra is the number of placeholder parameters introduced in e,
 398  * but not (yet) in F or C.
 399  *
 400  * If there is only one factor left, F is intersected with the
 401  * context C, the placeholder variables are added, and then
 402  * e is summed over the resulting parametric polytope.
 403  *
 404  * If there is more than one factor left, we create two polynomials
 405  * in a new placeholder variable (which is placed after the regular
 406  * parameters, but before any previously introduced placeholder
 407  * variables) that has the factors of the variables in the first
 408  * factor of F and the factor of the remaining variables of
 409  * each term as its coefficients.
 410  * These two polynomials are then summed over their domains
 411  * and afterwards the results are combined and the placeholder
 412  * variable is removed again.
 413  */
 414 static evalue *sum_factors(Polyhedron *F, Polyhedron *C, evalue *e,
 415                                      unsigned nvar, unsigned extra,
 416                                      struct barvinok_options *options)
 417 {
 418     Polyhedron *P = F;
 419     unsigned nparam = C->Dimension;
 420     unsigned F_var = F->Dimension - C->Dimension;
 421     int i, n;
 422     evalue *s1, *s2, *s;
 423     evalue *ph;
 424
 425     if (!F->next) {
 426         Polyhedron *CA = align_context(C, nvar+nparam, options->MaxRays);
 427         Polyhedron *P = DomainIntersection(F, CA, options->MaxRays);
 428         Polyhedron *Q = Polyhedron_Insert_Columns(P, nvar+nparam, extra);
 429         Polyhedron_Free(CA);
 430         Polyhedron_Free(F);
 431         Polyhedron_Free(P);
 432         evalue *sum = sum_base(Q, e, nvar, options);
 433         Polyhedron_Free(Q);
 434         return sum;
 435     }
 436
 437     n = evalue_count_terms(e, F_var, 0);
 438     ph = create_placeholder(n, nvar+nparam);
 439     evalue_shift_variables(e, nvar+nparam, 1);
 440     evalue_unzip_terms(e, ph, F_var, 0);
 441     evalue_shift_variables(e, nvar, -(nvar-F_var));
 442     evalue_reorder_terms(ph);
 443     evalue_shift_variables(ph, 0, -F_var);
 444
 445     s2 = sum_factors(F->next, C, ph, nvar-F_var, extra+1, options);
 446     evalue_free(ph);
 447     F->next = NULL;
 448     s1 = sum_factors(F, C, e, F_var, extra+1, options);
 449
 450     if (n == 1) {
 451         /* remove placeholder "polynomial" */
 452         reduce_evalue(s1);
 453         emul(s1, s2);
 454         evalue_free(s1);
 455         drop_parameters(s2, nparam, 1);
 456         return s2;
 457     }
 458
 459     s = evalue_zero();
 460     for (i = 0; i < n; ++i) {
 461         evalue *t1, *t2;
 462         t1 = extract_term(s1, nparam, i);
 463         t2 = extract_term(s2, nparam, i);
 464         emul(t1, t2);
 465         eadd(t2, s);
 466         evalue_free(t1);
 467         evalue_free(t2);
 468     }
 469     evalue_free(s1);
 470     evalue_free(s2);
 471
 472     drop_parameters(s, nparam, 1);
 473     return s;
 474 }
 475
 476 /* Perform summation over a product of factors F, obtained using
 477  * variable transformation T from the original problem specification.
 478  *
 479  * We first perform the corresponding transformation on the polynomial E,
 480  * compute the common context over all factors and then perform
 481  * the actual summation over the factors.
 482  */
 483 static evalue *sum_product(Polyhedron *F, evalue *E, Matrix *T, unsigned nparam,
 484                                      struct barvinok_options *options)
 485 {
 486     int i;
 487     Matrix *T2;
 488     unsigned nvar = T->NbRows;
 489     Polyhedron *C;
 490     evalue *sum;
 491
 492     assert(nvar == T->NbColumns);
 493     T2 = Matrix_Alloc(nvar+1, nvar+1);
 494     for (i = 0; i < nvar; ++i)
 495         Vector_Copy(T->p[i], T2->p[i], nvar);
 496     value_set_si(T2->p[nvar][nvar], 1);
 497
 498     transform_polynomial(E, T2, NULL, nvar, nparam, nvar, nparam);
 499
 500     C = Factor_Context(F, nparam, options->MaxRays);
 501     if (F->Dimension == nparam) {
 502         Polyhedron *T = F;
 503         F = F->next;
 504         T->next = NULL;
 505         Polyhedron_Free(T);
 506     }
 507     sum = sum_factors(F, C, E, nvar, 0, options);
 508
 509     Polyhedron_Free(C);
 510     Matrix_Free(T2);
 511     Matrix_Free(T);
 512     return sum;
 513 }
 514
 515 /* Add two constraints corresponding to floor = floor(e/d),
 516  *
 517  *       e - d t       >= 0
 518  *      -e + d t + d-1 >= 0
 519  *
 520  * e is assumed to be an affine expression.
 521  */
 522 Polyhedron *add_floor_var(Polyhedron *P, unsigned nvar, const evalue *floor,
 523                                      struct barvinok_options *options)
 524 {
 525     int i;
 526     unsigned dim = P->Dimension+1;
 527     Matrix *M = Matrix_Alloc(P->NbConstraints+2, 2+dim);
 528     Polyhedron *CP;
 529     Value *d = &M->p[0][1+nvar];
 530     evalue_extract_affine(floor, M->p[0]+1, M->p[0]+1+dim, d);
 531     value_oppose(*d, *d);
 532     value_set_si(M->p[0][0], 1);
 533     value_set_si(M->p[1][0], 1);
 534     Vector_Oppose(M->p[0]+1, M->p[1]+1, M->NbColumns-1);
 535     value_subtract(M->p[1][1+dim], M->p[1][1+dim], *d);
 536     value_decrement(M->p[1][1+dim], M->p[1][1+dim]);
 537
 538     for (i = 0; i < P->NbConstraints; ++i) {
 539         Vector_Copy(P->Constraint[i], M->p[i+2], 1+nvar);
 540         Vector_Copy(P->Constraint[i]+1+nvar, M->p[i+2]+1+nvar+1, dim-nvar-1+1);
 541     }
 542
 543     CP = Constraints2Polyhedron(M, options->MaxRays);
 544     Matrix_Free(M);
 545     return CP;
 546 }
 547
 548 static evalue *evalue_add(evalue *a, evalue *b)
 549 {
 550     if (!a)
 551         return b;
 552     if (!b)
 553         return a;
 554     eadd(a, b);
 555     evalue_free(a);
 556     return b;
 557 }
 558
 559 /* Compute sum of a step-polynomial over a polytope by grouping
 560  * terms containing the same floor-expressions and introducing
 561  * new variables for each such expression.
 562  * In particular, while there is any floor-expression left,
 563  * the step-polynomial is split into a polynomial containing
 564  * the expression, which is then converted to a new variable,
 565  * and a polynomial not containing the expression.
 566  */
 567 static evalue *sum_step_polynomial(Polyhedron *P, evalue *E, unsigned nvar,
 568                                      struct barvinok_options *options)
 569 {
 570     evalue *floor;
 571     evalue *cur = E;
 572     evalue *sum = NULL;
 573     evalue *t;
 574
 575     while ((floor = evalue_outer_floor(cur))) {
 576         Polyhedron *CP;
 577         evalue *converted;
 578         evalue *converted_floor;
 579
 580         /* Ignore floors that do not depend on variables. */
 581         if (value_notzero_p(floor->d) || floor->x.p->pos >= nvar+1)
 582             break;
 583
 584         converted = evalue_dup(cur);
 585         converted_floor = evalue_dup(floor);
 586         evalue_shift_variables(converted, nvar, 1);
 587         evalue_shift_variables(converted_floor, nvar, 1);
 588         evalue_replace_floor(converted, converted_floor, nvar);
 589         CP = add_floor_var(P, nvar, converted_floor, options);
 590         evalue_free(converted_floor);
 591         t = sum_step_polynomial(CP, converted, nvar+1, options);
 592         evalue_free(converted);
 593         Polyhedron_Free(CP);
 594         sum = evalue_add(t, sum);
 595
 596         if (cur == E)
 597             cur = evalue_dup(cur);
 598         evalue_drop_floor(cur, floor);
 599         evalue_free(floor);
 600     }
 601     if (floor) {
 602         evalue_floor2frac(cur);
 603         evalue_free(floor);
 604     }
 605
 606     if (EVALUE_IS_ZERO(*cur))
 607         t = NULL;
 608     else {
 609         Matrix *T;
 610         unsigned nparam = P->Dimension - nvar;
 611         Polyhedron *F = Polyhedron_Factor(P, nparam, &T, options->MaxRays);
 612         if (!F)
 613             t = sum_base(P, cur, nvar, options);
 614         else {
 615             if (cur == E)
 616                 cur = evalue_dup(cur);
 617             t = sum_product(F, cur, T, nparam, options);
 618         }
 619     }
 620
 621     if (E != cur)
 622         evalue_free(cur);
 623
 624     return evalue_add(t, sum);
 625 }
 626
 627 evalue *barvinok_sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
 628                                      struct evalue_section_array *sections,
 629                                      struct barvinok_options *options)
 630 {
 631     if (P->NbEq)
 632         return sum_over_polytope_with_equalities(P, E, nvar, sections, options);
 633
 634     if (options->summation == BV_SUM_BERNOULLI)
 635         return bernoulli_summate(P, E, nvar, sections, options);
 636     else if (options->summation == BV_SUM_BOX)
 637         return box_summate(P, E, nvar, options->MaxRays);
 638
 639     evalue_frac2floor2(E, 0);
 640
 641     return sum_step_polynomial(P, E, nvar, options);
 642 }
 643
 644 evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options)
 645 {
 646     int i;
 647     struct evalue_section_array sections;
 648     evalue *sum;
 649
 650     assert(nvar >= 0);
 651     if (nvar == 0 || EVALUE_IS_ZERO(*e))
 652         return evalue_dup(e);
 653
 654     assert(value_zero_p(e->d));
 655     assert(e->x.p->type == partition);
 656
 657     evalue_section_array_init(&sections);
 658     sum = evalue_zero();
 659
 660     for (i = 0; i < e->x.p->size/2; ++i) {
 661         Polyhedron *D;
 662         for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
 663             Polyhedron *next = D->next;
 664             evalue *tmp;
 665             D->next = NULL;
 666
 667             tmp = barvinok_sum_over_polytope(D, &e->x.p->arr[2*i+1], nvar,
 668                                              &sections, options);
 669             assert(tmp);
 670             eadd(tmp, sum);
 671             evalue_free(tmp);
 672
 673             D->next = next;
 674         }
 675     }
 676
 677     free(sections.s);
 678
 679     reduce_evalue(sum);
 680     return sum;
 681 }
 682
 683 evalue *evalue_sum(evalue *E, int nvar, unsigned MaxRays)
 684 {
 685     evalue *sum;
 686     struct barvinok_options *options = barvinok_options_new_with_defaults();
 687     options->MaxRays = MaxRays;
 688     sum = barvinok_summate(E, nvar, options);
 689     barvinok_options_free(options);
 690     return sum;
 691 }
 692
 693 evalue *esum(evalue *e, int nvar)
 694 {
 695     evalue *sum;
 696     struct barvinok_options *options = barvinok_options_new_with_defaults();
 697     sum = barvinok_summate(e, nvar, options);
 698     barvinok_options_free(options);
 699     return sum;
 700 }
 701
 702 /* Turn unweighted counting problem into "weighted" counting problem
 703  * with weight equal to 1 and call barvinok_summate on this weighted problem.
 704  */
 705 evalue *barvinok_summate_unweighted(Polyhedron *P, Polyhedron *C,
 706                                     struct barvinok_options *options)
 707 {
 708     Polyhedron *CA, *D;
 709     evalue e;
 710     evalue *sum;
 711
 712     if (emptyQ(P) || emptyQ(C))
 713         return evalue_zero();
 714
 715     CA = align_context(C, P->Dimension, options->MaxRays);
 716     D = DomainIntersection(P, CA, options->MaxRays);
 717     Domain_Free(CA);
 718
 719     if (emptyQ(D)) {
 720         Domain_Free(D);
 721         return evalue_zero();
 722     }
 723
 724     value_init(e.d);
 725     e.x.p = new_enode(partition, 2, P->Dimension);
 726     EVALUE_SET_DOMAIN(e.x.p->arr[0], D);
 727     evalue_set_si(&e.x.p->arr[1], 1, 1);
 728     sum = barvinok_summate(&e, P->Dimension - C->Dimension, options);
 729     free_evalue_refs(&e);
 730     return sum;
 731 }