bernoulli.c

   1 #include <assert.h>
   2 #include <barvinok/options.h>
   3 #include <barvinok/util.h>
   4 #include "bernoulli.h"
   5
   6 #define ALLOC(type) (type*)malloc(sizeof(type))
   7 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
   8 #define REALLOCN(ptr,type,n) (type*)realloc(ptr, (n) * sizeof(type))
   9
  10 static struct bernoulli_coef bernoulli_coef;
  11 static struct poly_list bernoulli;
  12 static struct poly_list faulhaber;
  13
  14 struct bernoulli_coef *bernoulli_coef_compute(int n)
  15 {
  16     int i, j;
  17     Value factor, tmp;
  18
  19     if (n < bernoulli_coef.n)
  20         return &bernoulli_coef;
  21
  22     if (n >= bernoulli_coef.size) {
  23         int size = 3*(n + 5)/2;
  24         Vector *b;
  25
  26         b = Vector_Alloc(size);
  27         if (bernoulli_coef.n) {
  28             Vector_Copy(bernoulli_coef.num->p, b->p, bernoulli_coef.n);
  29             Vector_Free(bernoulli_coef.num);
  30         }
  31         bernoulli_coef.num = b;
  32         b = Vector_Alloc(size);
  33         if (bernoulli_coef.n) {
  34             Vector_Copy(bernoulli_coef.den->p, b->p, bernoulli_coef.n);
  35             Vector_Free(bernoulli_coef.den);
  36         }
  37         bernoulli_coef.den = b;
  38         b = Vector_Alloc(size);
  39         if (bernoulli_coef.n) {
  40             Vector_Copy(bernoulli_coef.lcm->p, b->p, bernoulli_coef.n);
  41             Vector_Free(bernoulli_coef.lcm);
  42         }
  43         bernoulli_coef.lcm = b;
  44
  45         bernoulli_coef.size = size;
  46     }
  47     value_init(factor);
  48     value_init(tmp);
  49     for (i = bernoulli_coef.n; i <= n; ++i) {
  50         if (i == 0) {
  51             value_set_si(bernoulli_coef.num->p[0], 1);
  52             value_set_si(bernoulli_coef.den->p[0], 1);
  53             value_set_si(bernoulli_coef.lcm->p[0], 1);
  54             continue;
  55         }
  56         value_set_si(bernoulli_coef.num->p[i], 0);
  57         value_set_si(factor, -(i+1));
  58         for (j = i-1; j >= 0; --j) {
  59             mpz_mul_ui(factor, factor, j+1);
  60             mpz_divexact_ui(factor, factor, i+1-j);
  61             value_division(tmp, bernoulli_coef.lcm->p[i-1],
  62                            bernoulli_coef.den->p[j]);
  63             value_multiply(tmp, tmp, bernoulli_coef.num->p[j]);
  64             value_multiply(tmp, tmp, factor);
  65             value_addto(bernoulli_coef.num->p[i], bernoulli_coef.num->p[i], tmp);
  66         }
  67         mpz_mul_ui(bernoulli_coef.den->p[i], bernoulli_coef.lcm->p[i-1], i+1);
  68         value_gcd(tmp, bernoulli_coef.num->p[i], bernoulli_coef.den->p[i]);
  69         if (value_notone_p(tmp)) {
  70             value_division(bernoulli_coef.num->p[i],
  71                             bernoulli_coef.num->p[i], tmp);
  72             value_division(bernoulli_coef.den->p[i],
  73                             bernoulli_coef.den->p[i], tmp);
  74         }
  75         value_lcm(bernoulli_coef.lcm->p[i],
  76                   bernoulli_coef.lcm->p[i-1], bernoulli_coef.den->p[i]);
  77     }
  78     bernoulli_coef.n = n+1;
  79     value_clear(factor);
  80     value_clear(tmp);
  81
  82     return &bernoulli_coef;
  83 }
  84
  85 /*
  86  * Compute either Bernoulli B_n or Faulhaber F_n polynomials.
  87  *
  88  * B_n =         sum_{k=0}^n {  n  \choose k } b_k x^{n-k}
  89  * F_n = 1/(n+1) sum_{k=0}^n { n+1 \choose k } b_k x^{n+1-k}
  90  */
  91 static struct poly_list *bernoulli_faulhaber_compute(int n, struct poly_list *pl,
  92                                                      int faulhaber)
  93 {
  94     int i, j;
  95     Value factor;
  96     struct bernoulli_coef *bc;
  97
  98     if (n < pl->n)
  99         return pl;
 100
 101     if (n >= pl->size) {
 102         int size = 3*(n + 5)/2;
 103         Vector **poly;
 104
 105         poly = ALLOCN(Vector *, size);
 106         for (i = 0; i < pl->n; ++i)
 107             poly[i] = pl->poly[i];
 108         free(pl->poly);
 109         pl->poly = poly;
 110
 111         pl->size = size;
 112     }
 113
 114     bc = bernoulli_coef_compute(n);
 115
 116     value_init(factor);
 117     for (i = pl->n; i <= n; ++i) {
 118         pl->poly[i] = Vector_Alloc(i+faulhaber+2);
 119         value_assign(pl->poly[i]->p[i+faulhaber], bc->lcm->p[i]);
 120         if (faulhaber)
 121             mpz_mul_ui(pl->poly[i]->p[i+2], bc->lcm->p[i], i+1);
 122         else
 123             value_assign(pl->poly[i]->p[i+1], bc->lcm->p[i]);
 124         value_set_si(factor, 1);
 125         for (j = 1; j <= i; ++j) {
 126             mpz_mul_ui(factor, factor, i+faulhaber+1-j);
 127             mpz_divexact_ui(factor, factor, j);
 128             value_division(pl->poly[i]->p[i+faulhaber-j],
 129                            bc->lcm->p[i], bc->den->p[j]);
 130             value_multiply(pl->poly[i]->p[i+faulhaber-j],
 131                            pl->poly[i]->p[i+faulhaber-j], bc->num->p[j]);
 132             value_multiply(pl->poly[i]->p[i+faulhaber-j],
 133                            pl->poly[i]->p[i+faulhaber-j], factor);
 134         }
 135         Vector_Normalize(pl->poly[i]->p, i+faulhaber+2);
 136     }
 137     value_clear(factor);
 138     pl->n = n+1;
 139
 140     return pl;
 141 }
 142
 143 struct poly_list *bernoulli_compute(int n)
 144 {
 145     return bernoulli_faulhaber_compute(n, &bernoulli, 0);
 146 }
 147
 148 struct poly_list *faulhaber_compute(int n)
 149 {
 150     return bernoulli_faulhaber_compute(n, &faulhaber, 1);
 151 }
 152
 153 /* shift variables in polynomial one down */
 154 static void shift(evalue *e)
 155 {
 156     int i;
 157     if (value_notzero_p(e->d))
 158         return;
 159     assert(e->x.p->type == polynomial);
 160     assert(e->x.p->pos > 1);
 161     e->x.p->pos--;
 162     for (i = 0; i < e->x.p->size; ++i)
 163         shift(&e->x.p->arr[i]);
 164 }
 165
 166 /* shift variables in polynomial n up */
 167 static void unshift(evalue *e, unsigned n)
 168 {
 169     int i;
 170     if (value_notzero_p(e->d))
 171         return;
 172     assert(e->x.p->type == polynomial || e->x.p->type == fractional);
 173     if (e->x.p->type == polynomial)
 174         e->x.p->pos += n;
 175     for (i = 0; i < e->x.p->size; ++i)
 176         unshift(&e->x.p->arr[i], n);
 177 }
 178
 179 static evalue *shifted_copy(evalue *src)
 180 {
 181     evalue *e = ALLOC(evalue);
 182     value_init(e->d);
 183     evalue_copy(e, src);
 184     shift(e);
 185     return e;
 186 }
 187
 188 static evalue *power_sums(struct poly_list *faulhaber, evalue *poly,
 189                           Vector *arg, Value denom, int neg, int alt_neg)
 190 {
 191     int i;
 192     evalue *base = affine2evalue(arg->p, denom, arg->Size-1);
 193     evalue *sum = evalue_zero();
 194
 195     for (i = 1; i < poly->x.p->size; ++i) {
 196         evalue *term = evalue_polynomial(faulhaber->poly[i], base);
 197         evalue *factor = shifted_copy(&poly->x.p->arr[i]);
 198         emul(factor, term);
 199         if (alt_neg && (i % 2))
 200             evalue_negate(term);
 201         eadd(term, sum);
 202         evalue_free(factor);
 203         evalue_free(term);
 204     }
 205     if (neg)
 206         evalue_negate(sum);
 207     evalue_free(base);
 208
 209     return sum;
 210 }
 211
 212 /* Given a constraint (cst_affine) a x + b y + c >= 0, compate a constraint (c)
 213  * +- (b y + c) +- a >=,> 0
 214  * ^            ^      ^
 215  * |            |      strict
 216  * sign_affine  sign_cst
 217  */
 218 static void bound_constraint(Value *c, unsigned dim,
 219                              Value *cst_affine,
 220                              int sign_affine, int sign_cst, int strict)
 221 {
 222     if (sign_affine == -1)
 223         Vector_Oppose(cst_affine+1, c, dim+1);
 224     else
 225         Vector_Copy(cst_affine+1, c, dim+1);
 226
 227     if (sign_cst == -1)
 228         value_subtract(c[dim], c[dim], cst_affine[0]);
 229     else if (sign_cst == 1)
 230         value_addto(c[dim], c[dim], cst_affine[0]);
 231
 232     if (strict)
 233         value_decrement(c[dim], c[dim]);
 234 }
 235
 236 struct Bernoulli_data {
 237     unsigned MaxRays;
 238     struct evalue_section *s;
 239     int size;
 240     int ns;
 241     evalue *e;
 242 };
 243
 244 static evalue *compute_poly_u(evalue *poly_u, Value *upper, Vector *row,
 245                               unsigned dim, Value tmp,
 246                               struct poly_list *faulhaber,
 247                               struct Bernoulli_data *data)
 248 {
 249     if (poly_u)
 250         return poly_u;
 251     Vector_Copy(upper+2, row->p, dim+1);
 252     value_oppose(tmp, upper[1]);
 253     value_addto(row->p[dim], row->p[dim], tmp);
 254     return power_sums(faulhaber, data->e, row, tmp, 0, 0);
 255 }
 256
 257 static evalue *compute_poly_l(evalue *poly_l, Value *lower, Vector *row,
 258                               unsigned dim,
 259                               struct poly_list *faulhaber,
 260                               struct Bernoulli_data *data)
 261 {
 262     if (poly_l)
 263         return poly_l;
 264     Vector_Copy(lower+2, row->p, dim+1);
 265     value_addto(row->p[dim], row->p[dim], lower[1]);
 266     return power_sums(faulhaber, data->e, row, lower[1], 0, 1);
 267 }
 268
 269 static void Bernoulli_init(unsigned n, void *cb_data)
 270 {
 271     struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
 272     int cases = 5;
 273
 274     if (cases * n <= data->size)
 275         return;
 276
 277     data->size = cases * (n + 16);
 278     data->s = REALLOCN(data->s, struct evalue_section, data->size);
 279 }
 280
 281 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
 282 {
 283     struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
 284     Matrix *M2;
 285     Polyhedron *T;
 286     evalue *factor = NULL;
 287     evalue *linear = NULL;
 288     int constant = 0;
 289     Value tmp;
 290     unsigned dim = M->NbColumns-2;
 291     Vector *row;
 292     int cases = 5;
 293
 294     assert(lower);
 295     assert(upper);
 296     assert(data->ns + cases <= data->size);
 297
 298     M2 = Matrix_Copy(M);
 299     T = Constraints2Polyhedron(M2, data->MaxRays);
 300     Matrix_Free(M2);
 301
 302     POL_ENSURE_VERTICES(T);
 303     if (emptyQ(T)) {
 304         Polyhedron_Free(T);
 305         return;
 306     }
 307
 308     assert(lower != upper);
 309
 310     row = Vector_Alloc(dim+1);
 311     value_init(tmp);
 312     if (value_notzero_p(data->e->d)) {
 313         factor = data->e;
 314         constant = 1;
 315     } else {
 316         assert(data->e->x.p->type == polynomial);
 317         if (data->e->x.p->pos > 1) {
 318             factor = shifted_copy(data->e);
 319             constant = 1;
 320         } else
 321             factor = shifted_copy(&data->e->x.p->arr[0]);
 322     }
 323     if (!EVALUE_IS_ZERO(*factor)) {
 324         value_absolute(tmp, upper[1]);
 325         /* upper - lower */
 326         Vector_Combine(lower+2, upper+2, row->p, tmp, lower[1], dim+1);
 327         value_multiply(tmp, tmp, lower[1]);
 328         /* upper - lower + 1 */
 329         value_addto(row->p[dim], row->p[dim], tmp);
 330
 331         linear = affine2evalue(row->p, tmp, dim);
 332         emul(factor, linear);
 333     } else
 334         linear = evalue_zero();
 335
 336     if (constant) {
 337         data->s[data->ns].E = linear;
 338         data->s[data->ns].D = T;
 339         ++data->ns;
 340     } else {
 341         evalue *poly_u = NULL, *poly_l = NULL;
 342         Polyhedron *D;
 343         struct poly_list *faulhaber;
 344         assert(data->e->x.p->type == polynomial);
 345         assert(data->e->x.p->pos == 1);
 346         faulhaber = faulhaber_compute(data->e->x.p->size-1);
 347         /* lower is the constraint
 348          *                       b i - l' >= 0          i >= l'/b = l
 349          * upper is the constraint
 350          *                      -a i + u' >= 0          i <= u'/a = u
 351          */
 352         M2 = Matrix_Alloc(3, 2+T->Dimension);
 353         value_set_si(M2->p[0][0], 1);
 354         value_set_si(M2->p[1][0], 1);
 355         value_set_si(M2->p[2][0], 1);
 356         /* Case 1:
 357          *              1 <= l          l' - b >= 0
 358          */
 359         bound_constraint(M2->p[0]+1, T->Dimension, lower+1, -1, -1, 0);
 360         D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
 361         POL_ENSURE_VERTICES(D);
 362         if (emptyQ2(D))
 363             Polyhedron_Free(D);
 364         else {
 365             evalue *extra;
 366             poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
 367                                         faulhaber, data);
 368             Vector_Oppose(lower+2, row->p, dim+1);
 369             extra = power_sums(faulhaber, data->e, row, lower[1], 1, 0);
 370             eadd(poly_u, extra);
 371             eadd(linear, extra);
 372
 373             data->s[data->ns].E = extra;
 374             data->s[data->ns].D = D;
 375             ++data->ns;
 376         }
 377
 378         /* Case 2:
 379          *              1 <= -u         -u' - a >= 0
 380          */
 381         bound_constraint(M2->p[0]+1, T->Dimension, upper+1, -1, 1, 0);
 382         D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
 383         POL_ENSURE_VERTICES(D);
 384         if (emptyQ2(D))
 385             Polyhedron_Free(D);
 386         else {
 387             evalue *extra;
 388             poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
 389             Vector_Oppose(upper+2, row->p, dim+1);
 390             value_oppose(tmp, upper[1]);
 391             extra = power_sums(faulhaber, data->e, row, tmp, 1, 1);
 392             eadd(poly_l, extra);
 393             eadd(linear, extra);
 394
 395             data->s[data->ns].E = extra;
 396             data->s[data->ns].D = D;
 397             ++data->ns;
 398         }
 399
 400         /* Case 3:
 401          *              u >= 0          u' >= 0
 402          *              -l >= 0         -l' >= 0
 403          */
 404         bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, 0, 0);
 405         bound_constraint(M2->p[1]+1, T->Dimension, lower+1, 1, 0, 0);
 406         D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
 407         POL_ENSURE_VERTICES(D);
 408         if (emptyQ2(D))
 409             Polyhedron_Free(D);
 410         else {
 411             poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
 412             poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
 413                                         faulhaber, data);
 414
 415             data->s[data->ns].E = ALLOC(evalue);
 416             value_init(data->s[data->ns].E->d);
 417             evalue_copy(data->s[data->ns].E, poly_u);
 418             eadd(poly_l, data->s[data->ns].E);
 419             eadd(linear, data->s[data->ns].E);
 420             data->s[data->ns].D = D;
 421             ++data->ns;
 422         }
 423
 424         /* Case 4:
 425          *              l < 1           -l' + b - 1 >= 0
 426          *              0 < l           l' - 1 >= 0
 427          *              u >= 1          u' - a >= 0
 428          */
 429         bound_constraint(M2->p[0]+1, T->Dimension, lower+1, 1, 1, 1);
 430         bound_constraint(M2->p[1]+1, T->Dimension, lower+1, -1, 0, 1);
 431         bound_constraint(M2->p[2]+1, T->Dimension, upper+1, 1, 1, 0);
 432         D = AddConstraints(M2->p_Init, 3, T, data->MaxRays);
 433         POL_ENSURE_VERTICES(D);
 434         if (emptyQ2(D))
 435             Polyhedron_Free(D);
 436         else {
 437             poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
 438                                         faulhaber, data);
 439
 440             eadd(linear, poly_u);
 441             data->s[data->ns].E = poly_u;
 442             poly_u = NULL;
 443             data->s[data->ns].D = D;
 444             ++data->ns;
 445         }
 446
 447         /* Case 5:
 448          *              -u < 1          u' + a - 1 >= 0
 449          *              0 < -u          -u' - 1 >= 0
 450          *              l <= -1         -l' - b >= 0
 451          */
 452         bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, -1, 1);
 453         bound_constraint(M2->p[1]+1, T->Dimension, upper+1, -1, 0, 1);
 454         bound_constraint(M2->p[2]+1, T->Dimension, lower+1, 1, -1, 0);
 455         D = AddConstraints(M2->p_Init, 3, T, data->MaxRays);
 456         POL_ENSURE_VERTICES(D);
 457         if (emptyQ2(D))
 458             Polyhedron_Free(D);
 459         else {
 460             poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
 461
 462             eadd(linear, poly_l);
 463             data->s[data->ns].E = poly_l;
 464             poly_l = NULL;
 465             data->s[data->ns].D = D;
 466             ++data->ns;
 467         }
 468
 469         Matrix_Free(M2);
 470         Polyhedron_Free(T);
 471         if (poly_l)
 472             evalue_free(poly_l);
 473         if (poly_u)
 474             evalue_free(poly_u);
 475         evalue_free(linear);
 476     }
 477     if (factor != data->e)
 478         evalue_free(factor);
 479     value_clear(tmp);
 480     Vector_Free(row);
 481 }
 482
 483 /* Looks for variable with integer bounds, i.e., with coefficients 0, 1 or -1.
 484  * Returns 1 if such a variable is found and puts it in the first position,
 485  * possibly changing *P_p and *E_p.
 486  */
 487 static int find_integer_bounds(Polyhedron **P_p, evalue **E_p, unsigned nvar)
 488 {
 489     Polyhedron *P = *P_p;
 490     evalue *E = *E_p;
 491     unsigned dim = P->Dimension;
 492     int i, j;
 493
 494     for (i = 0; i < nvar; ++i) {
 495         for (j = 0; j < P->NbConstraints; ++j) {
 496             if (value_zero_p(P->Constraint[j][1+i]))
 497                 continue;
 498             if (value_one_p(P->Constraint[j][1+i]))
 499                 continue;
 500             if (value_mone_p(P->Constraint[j][1+i]))
 501                 continue;
 502             break;
 503         }
 504         if (j == P->NbConstraints)
 505             break;
 506     }
 507     if (i == nvar)
 508         return 0;
 509     if (i == 0)
 510         return 1;
 511     P = Polyhedron_Copy(P);
 512     Polyhedron_ExchangeColumns(P, 1, 1+i);
 513     *P_p = P;
 514
 515     if (value_zero_p(E->d)) {
 516         evalue **subs;
 517         subs = ALLOCN(evalue *, dim);
 518         for (j = 0; j < dim; ++j)
 519             subs[j] = evalue_var(j);
 520         E = subs[0];
 521         subs[0] = subs[i];
 522         subs[i] = E;
 523         E = evalue_dup(*E_p);
 524         evalue_substitute(E, subs);
 525         for (j = 0; j < dim; ++j)
 526             evalue_free(subs[j]);
 527         free(subs);
 528         *E_p = E;
 529     }
 530
 531     return 1;
 532 }
 533
 534 static evalue *sum_over_polytope_base(Polyhedron *P, evalue *E, unsigned nvar,
 535                                       struct Bernoulli_data *data,
 536                                       struct barvinok_options *options)
 537 {
 538     evalue *res;
 539
 540     assert(P->NbEq == 0);
 541
 542     data->ns = 0;
 543     data->e = E;
 544
 545     for_each_lower_upper_bound(P, Bernoulli_init, Bernoulli_cb, data);
 546
 547     res = evalue_from_section_array(data->s, data->ns);
 548
 549     if (nvar > 1) {
 550         evalue *tmp = Bernoulli_sum_evalue(res, nvar-1, options);
 551         evalue_free(res);
 552         res = tmp;
 553     }
 554
 555     return res;
 556 }
 557
 558 static evalue *sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
 559                                  struct Bernoulli_data *data,
 560                                  struct barvinok_options *options);
 561
 562 static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E,
 563                                  unsigned nvar,
 564                                  struct Bernoulli_data *data,
 565                                  struct barvinok_options *options)
 566 {
 567     unsigned dim = P->Dimension;
 568     unsigned new_dim, new_nparam;
 569     Matrix *T = NULL, *CP = NULL;
 570     evalue **subs;
 571     evalue *sum;
 572     int j;
 573
 574     assert(P->NbEq > 0);
 575
 576     remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays);
 577
 578     if (emptyQ(P)) {
 579         Polyhedron_Free(P);
 580         return evalue_zero();
 581     }
 582
 583     new_nparam = CP ? CP->NbColumns-1 : dim - nvar;
 584     new_dim = T ? T->NbColumns-1 : nvar + new_nparam;
 585
 586     /* We can avoid these substitutions if E is a constant */
 587     subs = ALLOCN(evalue *, dim);
 588     for (j = 0; j < nvar; ++j) {
 589         if (T)
 590             subs[j] = affine2evalue(T->p[j], T->p[nvar+new_nparam][new_dim],
 591                                     new_dim);
 592         else
 593             subs[j] = evalue_var(j);
 594     }
 595     for (j = 0; j < dim-nvar; ++j) {
 596         if (CP)
 597             subs[nvar+j] = affine2evalue(CP->p[j], CP->p[dim-nvar][new_nparam],
 598                                          new_nparam);
 599         else
 600             subs[nvar+j] = evalue_var(j);
 601         unshift(subs[nvar+j], new_dim-new_nparam);
 602     }
 603
 604     E = evalue_dup(E);
 605     evalue_substitute(E, subs);
 606     reduce_evalue(E);
 607
 608     for (j = 0; j < dim; ++j)
 609         evalue_free(subs[j]);
 610     free(subs);
 611
 612     if (new_dim-new_nparam > 0) {
 613         sum = sum_over_polytope(P, E, new_dim-new_nparam, data, options);
 614         evalue_free(E);
 615         Polyhedron_Free(P);
 616     } else {
 617         sum = ALLOC(evalue);
 618         value_init(sum->d);
 619         sum->x.p = new_enode(partition, 2, new_dim);
 620         EVALUE_SET_DOMAIN(sum->x.p->arr[0], P);
 621         value_clear(sum->x.p->arr[1].d);
 622         sum->x.p->arr[1] = *E;
 623         free(E);
 624     }
 625
 626     if (CP) {
 627         evalue_backsubstitute(sum, CP, options->MaxRays);
 628         Matrix_Free(CP);
 629     }
 630
 631     if (T)
 632         Matrix_Free(T);
 633
 634     return sum;
 635 }
 636
 637 static evalue *sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
 638                                  struct Bernoulli_data *data,
 639                                  struct barvinok_options *options)
 640 {
 641     if (P->NbEq)
 642         return sum_over_polytope_with_equalities(P, E, nvar, data, options);
 643     else
 644         return sum_over_polytope_base(P, E, nvar, data, options);
 645 }
 646
 647 evalue *Bernoulli_sum_evalue(evalue *e, unsigned nvar,
 648                              struct barvinok_options *options)
 649 {
 650     struct Bernoulli_data data;
 651     int i, j;
 652     evalue *sum = evalue_zero();
 653
 654     if (EVALUE_IS_ZERO(*e))
 655         return sum;
 656
 657     if (nvar == 0) {
 658         eadd(e, sum);
 659         return sum;
 660     }
 661
 662     assert(value_zero_p(e->d));
 663     assert(e->x.p->type == partition);
 664
 665     data.size = 16;
 666     data.s = ALLOCN(struct evalue_section, data.size);
 667     data.MaxRays = options->MaxRays;
 668
 669     for (i = 0; i < e->x.p->size/2; ++i) {
 670         Polyhedron *D;
 671         for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
 672             evalue *E = &e->x.p->arr[2*i+1];
 673             Polyhedron *P = D;
 674             Polyhedron *next = D->next;
 675             evalue *tmp;
 676             int integer_bounds;
 677
 678             P->next = NULL;
 679
 680             integer_bounds = find_integer_bounds(&P, &E, nvar);
 681             if (options->approximation_method == BV_APPROX_NONE &&
 682                 !integer_bounds) {
 683                 evalue_free(sum);
 684                 sum = NULL;
 685             } else {
 686                 evalue *tmp = sum_over_polytope(P, E, nvar, &data, options);
 687                 eadd(tmp, sum);
 688                 evalue_free(tmp);
 689             }
 690
 691             if (P != D)
 692                 Polyhedron_Free(P);
 693             if (E != &e->x.p->arr[2*i+1])
 694                 evalue_free(E);
 695
 696             D->next = next;;
 697
 698             if (!sum)
 699                 break;
 700         }
 701
 702         if (!sum)
 703             break;
 704     }
 705
 706     free(data.s);
 707
 708     if (sum)
 709         reduce_evalue(sum);
 710     return sum;
 711 }
 712
 713 evalue *Bernoulli_sum(Polyhedron *P, Polyhedron *C,
 714                         struct barvinok_options *options)
 715 {
 716     Polyhedron *CA, *D;
 717     evalue e;
 718     evalue *sum;
 719
 720     if (emptyQ(P) || emptyQ(C))
 721         return evalue_zero();
 722
 723     CA = align_context(C, P->Dimension, options->MaxRays);
 724     D = DomainIntersection(P, CA, options->MaxRays);
 725     Domain_Free(CA);
 726
 727     if (emptyQ(D)) {
 728         Domain_Free(D);
 729         return evalue_zero();
 730     }
 731
 732     value_init(e.d);
 733     e.x.p = new_enode(partition, 2, P->Dimension);
 734     EVALUE_SET_DOMAIN(e.x.p->arr[0], D);
 735     evalue_set_si(&e.x.p->arr[1], 1, 1);
 736     sum = Bernoulli_sum_evalue(&e, P->Dimension - C->Dimension, options);
 737     free_evalue_refs(&e);
 738     return sum;
 739 }