bernoulli.c

   1 #include <assert.h>
   2 #include <barvinok/barvinok.h>
   3 #include <barvinok/options.h>
   4 #include <barvinok/util.h>
   5 #include "bernoulli.h"
   6 #include "lattice_point.h"
   7 #include "section_array.h"
   8
   9 #define ALLOC(type) (type*)malloc(sizeof(type))
  10 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  11
  12 static struct bernoulli_coef bernoulli_coef;
  13 static struct poly_list bernoulli;
  14 static struct poly_list faulhaber;
  15
  16 struct bernoulli_coef *bernoulli_coef_compute(int n)
  17 {
  18     int i, j;
  19     Value factor, tmp;
  20
  21     if (n < bernoulli_coef.n)
  22         return &bernoulli_coef;
  23
  24     if (n >= bernoulli_coef.size) {
  25         int size = 3*(n + 5)/2;
  26         Vector *b;
  27
  28         b = Vector_Alloc(size);
  29         if (bernoulli_coef.n) {
  30             Vector_Copy(bernoulli_coef.num->p, b->p, bernoulli_coef.n);
  31             Vector_Free(bernoulli_coef.num);
  32         }
  33         bernoulli_coef.num = b;
  34         b = Vector_Alloc(size);
  35         if (bernoulli_coef.n) {
  36             Vector_Copy(bernoulli_coef.den->p, b->p, bernoulli_coef.n);
  37             Vector_Free(bernoulli_coef.den);
  38         }
  39         bernoulli_coef.den = b;
  40         b = Vector_Alloc(size);
  41         if (bernoulli_coef.n) {
  42             Vector_Copy(bernoulli_coef.lcm->p, b->p, bernoulli_coef.n);
  43             Vector_Free(bernoulli_coef.lcm);
  44         }
  45         bernoulli_coef.lcm = b;
  46
  47         bernoulli_coef.size = size;
  48     }
  49     value_init(factor);
  50     value_init(tmp);
  51     for (i = bernoulli_coef.n; i <= n; ++i) {
  52         if (i == 0) {
  53             value_set_si(bernoulli_coef.num->p[0], 1);
  54             value_set_si(bernoulli_coef.den->p[0], 1);
  55             value_set_si(bernoulli_coef.lcm->p[0], 1);
  56             continue;
  57         }
  58         value_set_si(bernoulli_coef.num->p[i], 0);
  59         value_set_si(factor, -(i+1));
  60         for (j = i-1; j >= 0; --j) {
  61             mpz_mul_ui(factor, factor, j+1);
  62             mpz_divexact_ui(factor, factor, i+1-j);
  63             value_division(tmp, bernoulli_coef.lcm->p[i-1],
  64                            bernoulli_coef.den->p[j]);
  65             value_multiply(tmp, tmp, bernoulli_coef.num->p[j]);
  66             value_multiply(tmp, tmp, factor);
  67             value_addto(bernoulli_coef.num->p[i], bernoulli_coef.num->p[i], tmp);
  68         }
  69         mpz_mul_ui(bernoulli_coef.den->p[i], bernoulli_coef.lcm->p[i-1], i+1);
  70         value_gcd(tmp, bernoulli_coef.num->p[i], bernoulli_coef.den->p[i]);
  71         if (value_notone_p(tmp)) {
  72             value_division(bernoulli_coef.num->p[i],
  73                             bernoulli_coef.num->p[i], tmp);
  74             value_division(bernoulli_coef.den->p[i],
  75                             bernoulli_coef.den->p[i], tmp);
  76         }
  77         value_lcm(bernoulli_coef.lcm->p[i],
  78                   bernoulli_coef.lcm->p[i-1], bernoulli_coef.den->p[i]);
  79     }
  80     bernoulli_coef.n = n+1;
  81     value_clear(factor);
  82     value_clear(tmp);
  83
  84     return &bernoulli_coef;
  85 }
  86
  87 /*
  88  * Compute either Bernoulli B_n or Faulhaber F_n polynomials.
  89  *
  90  * B_n =         sum_{k=0}^n {  n  \choose k } b_k x^{n-k}
  91  * F_n = 1/(n+1) sum_{k=0}^n { n+1 \choose k } b_k x^{n+1-k}
  92  */
  93 static struct poly_list *bernoulli_faulhaber_compute(int n, struct poly_list *pl,
  94                                                      int faulhaber)
  95 {
  96     int i, j;
  97     Value factor;
  98     struct bernoulli_coef *bc;
  99
 100     if (n < pl->n)
 101         return pl;
 102
 103     if (n >= pl->size) {
 104         int size = 3*(n + 5)/2;
 105         Vector **poly;
 106
 107         poly = ALLOCN(Vector *, size);
 108         for (i = 0; i < pl->n; ++i)
 109             poly[i] = pl->poly[i];
 110         free(pl->poly);
 111         pl->poly = poly;
 112
 113         pl->size = size;
 114     }
 115
 116     bc = bernoulli_coef_compute(n);
 117
 118     value_init(factor);
 119     for (i = pl->n; i <= n; ++i) {
 120         pl->poly[i] = Vector_Alloc(i+faulhaber+2);
 121         value_assign(pl->poly[i]->p[i+faulhaber], bc->lcm->p[i]);
 122         if (faulhaber)
 123             mpz_mul_ui(pl->poly[i]->p[i+2], bc->lcm->p[i], i+1);
 124         else
 125             value_assign(pl->poly[i]->p[i+1], bc->lcm->p[i]);
 126         value_set_si(factor, 1);
 127         for (j = 1; j <= i; ++j) {
 128             mpz_mul_ui(factor, factor, i+faulhaber+1-j);
 129             mpz_divexact_ui(factor, factor, j);
 130             value_division(pl->poly[i]->p[i+faulhaber-j],
 131                            bc->lcm->p[i], bc->den->p[j]);
 132             value_multiply(pl->poly[i]->p[i+faulhaber-j],
 133                            pl->poly[i]->p[i+faulhaber-j], bc->num->p[j]);
 134             value_multiply(pl->poly[i]->p[i+faulhaber-j],
 135                            pl->poly[i]->p[i+faulhaber-j], factor);
 136         }
 137         Vector_Normalize(pl->poly[i]->p, i+faulhaber+2);
 138     }
 139     value_clear(factor);
 140     pl->n = n+1;
 141
 142     return pl;
 143 }
 144
 145 struct poly_list *bernoulli_compute(int n)
 146 {
 147     return bernoulli_faulhaber_compute(n, &bernoulli, 0);
 148 }
 149
 150 struct poly_list *faulhaber_compute(int n)
 151 {
 152     return bernoulli_faulhaber_compute(n, &faulhaber, 1);
 153 }
 154
 155 static evalue *shifted_copy(const evalue *src)
 156 {
 157     evalue *e = ALLOC(evalue);
 158     value_init(e->d);
 159     evalue_copy(e, src);
 160     evalue_shift_variables(e, -1);
 161     return e;
 162 }
 163
 164 /* Computes the argument for the Faulhaber polynomial.
 165  * If we are computing a polynomial approximation (exact == 0),
 166  * then this is simply arg/denom.
 167  * Otherwise, if neg == 0, then we are dealing with an upper bound
 168  * and we need to compute floor(arg/denom) = arg/denom - { arg/denom }
 169  * If neg == 1, then we are dealing with a lower bound
 170  * and we need to compute ceil(arg/denom) = arg/denom + { -arg/denom }
 171  *
 172  * Modifies arg (if exact == 1).
 173  */
 174 static evalue *power_sums_base(Vector *arg, Value denom, int neg, int exact)
 175 {
 176     evalue *fract;
 177     evalue *base = affine2evalue(arg->p, denom, arg->Size-1);
 178
 179     if (!exact || value_one_p(denom))
 180         return base;
 181
 182     if (neg)
 183         Vector_Oppose(arg->p, arg->p, arg->Size);
 184
 185     fract = fractional_part(arg->p, denom, arg->Size-1, NULL);
 186     if (!neg) {
 187         value_set_si(arg->p[0], -1);
 188         evalue_mul(fract, arg->p[0]);
 189     }
 190     eadd(fract, base);
 191     evalue_free(fract);
 192
 193     return base;
 194 }
 195
 196 static evalue *power_sums(struct poly_list *faulhaber, const evalue *poly,
 197                           Vector *arg, Value denom, int neg, int alt_neg,
 198                           int exact)
 199 {
 200     int i;
 201     evalue *base = power_sums_base(arg, denom, neg, exact);
 202     evalue *sum = evalue_zero();
 203
 204     for (i = 1; i < poly->x.p->size; ++i) {
 205         evalue *term = evalue_polynomial(faulhaber->poly[i], base);
 206         evalue *factor = shifted_copy(&poly->x.p->arr[i]);
 207         emul(factor, term);
 208         if (alt_neg && (i % 2))
 209             evalue_negate(term);
 210         eadd(term, sum);
 211         evalue_free(factor);
 212         evalue_free(term);
 213     }
 214     if (neg)
 215         evalue_negate(sum);
 216     evalue_free(base);
 217
 218     return sum;
 219 }
 220
 221 /* Given a constraint (cst_affine) a x + b y + c >= 0, compate a constraint (c)
 222  * +- (b y + c) +- a >=,> 0
 223  * ^            ^      ^
 224  * |            |      strict
 225  * sign_affine  sign_cst
 226  */
 227 static void bound_constraint(Value *c, unsigned dim,
 228                              Value *cst_affine,
 229                              int sign_affine, int sign_cst, int strict)
 230 {
 231     if (sign_affine == -1)
 232         Vector_Oppose(cst_affine+1, c, dim+1);
 233     else
 234         Vector_Copy(cst_affine+1, c, dim+1);
 235
 236     if (sign_cst == -1)
 237         value_subtract(c[dim], c[dim], cst_affine[0]);
 238     else if (sign_cst == 1)
 239         value_addto(c[dim], c[dim], cst_affine[0]);
 240
 241     if (strict)
 242         value_decrement(c[dim], c[dim]);
 243 }
 244
 245 struct Bernoulli_data {
 246     struct barvinok_options *options;
 247     struct evalue_section_array *sections;
 248     const evalue *e;
 249 };
 250
 251 static evalue *compute_poly_u(evalue *poly_u, Value *upper, Vector *row,
 252                               unsigned dim, Value tmp,
 253                               struct poly_list *faulhaber,
 254                               struct Bernoulli_data *data)
 255 {
 256     int exact = data->options->approximation_method == BV_APPROX_NONE;
 257     if (poly_u)
 258         return poly_u;
 259     Vector_Copy(upper+2, row->p, dim+1);
 260     value_oppose(tmp, upper[1]);
 261     value_addto(row->p[dim], row->p[dim], tmp);
 262     return power_sums(faulhaber, data->e, row, tmp, 0, 0, exact);
 263 }
 264
 265 static evalue *compute_poly_l(evalue *poly_l, Value *lower, Vector *row,
 266                               unsigned dim,
 267                               struct poly_list *faulhaber,
 268                               struct Bernoulli_data *data)
 269 {
 270     int exact = data->options->approximation_method == BV_APPROX_NONE;
 271     if (poly_l)
 272         return poly_l;
 273     Vector_Copy(lower+2, row->p, dim+1);
 274     value_addto(row->p[dim], row->p[dim], lower[1]);
 275     return power_sums(faulhaber, data->e, row, lower[1], 0, 1, exact);
 276 }
 277
 278 /* Compute sum of constant term.
 279  */
 280 static evalue *linear_term(const evalue *cst, Value *lower, Value *upper,
 281                            Vector *row, Value tmp, int exact)
 282 {
 283     evalue *linear;
 284     unsigned dim = row->Size - 1;
 285
 286     if (EVALUE_IS_ZERO(*cst))
 287         return evalue_zero();
 288
 289     if (!exact) {
 290         value_absolute(tmp, upper[1]);
 291         /* upper - lower */
 292         Vector_Combine(lower+2, upper+2, row->p, tmp, lower[1], dim+1);
 293         value_multiply(tmp, tmp, lower[1]);
 294         /* upper - lower + 1 */
 295         value_addto(row->p[dim], row->p[dim], tmp);
 296
 297         linear = affine2evalue(row->p, tmp, dim);
 298     } else {
 299         evalue *l;
 300
 301         value_absolute(tmp, upper[1]);
 302         Vector_Copy(upper+2, row->p, dim+1);
 303         value_addto(row->p[dim], row->p[dim], tmp);
 304         /* floor(upper+1) */
 305         linear = power_sums_base(row, tmp, 0, 1);
 306
 307         Vector_Copy(lower+2, row->p, dim+1);
 308         /* floor(-lower) */
 309         l = power_sums_base(row, lower[1], 0, 1);
 310
 311         /* floor(upper+1) + floor(-lower) = floor(upper) - ceil(lower) + 1 */
 312         eadd(l, linear);
 313         evalue_free(l);
 314     }
 315
 316     emul(cst, linear);
 317     return linear;
 318 }
 319
 320 static void Bernoulli_init(unsigned n, void *cb_data)
 321 {
 322     struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
 323     int exact = data->options->approximation_method == BV_APPROX_NONE;
 324     int cases = exact ? 3 : 5;
 325
 326     evalue_section_array_ensure(data->sections, cases * n);
 327 }
 328
 329 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data);
 330
 331 /*
 332  * This function requires that either the lower or the upper bound
 333  * represented by the constraints "lower" and "upper" is an integer
 334  * affine expression.
 335  * An affine substitution is performed to make this bound exactly
 336  * zero, ensuring that in the recursive call to Bernoulli_cb,
 337  * only one of the "cases" will apply.
 338  */
 339 static void transform_to_single_case(Matrix *M, Value *lower, Value *upper,
 340                                      struct Bernoulli_data *data)
 341 {
 342     unsigned dim = M->NbColumns-2;
 343     Vector *shadow;
 344     Value a, b;
 345     evalue **subs;
 346     const evalue *e = data->e;
 347     evalue *t;
 348     int i;
 349
 350     value_init(a);
 351     value_init(b);
 352     subs = ALLOCN(evalue *, dim+1);
 353     for (i = 0; i < dim; ++i)
 354         subs[1+i] = evalue_var(1+i);
 355     shadow = Vector_Alloc(2 * (2+dim+1));
 356     if (value_one_p(lower[1])) {
 357         /* Replace i by i + l' when b = 1 */
 358         value_set_si(shadow->p[0], 1);
 359         Vector_Oppose(lower+2, shadow->p+1, dim+1);
 360         subs[0] = affine2evalue(shadow->p, shadow->p[0], dim+1);
 361         /* new lower
 362          *              i >= 0
 363          * new upper
 364          *              (-a i + u') + a (-l') >= 0
 365          */
 366         value_assign(shadow->p[2+dim+1+1], upper[1]);
 367         value_oppose(a, upper[1]);
 368         value_set_si(b, 1);
 369         Vector_Combine(upper+2, lower+2, shadow->p+2+dim+1+2,
 370                        b, a, dim+1);
 371         upper = shadow->p+2+dim+1;
 372         lower = shadow->p;
 373         value_set_si(lower[1], 1);
 374         Vector_Set(lower+2, 0, dim+1);
 375     } else {
 376         /* Replace i by i + u' when a = 1 */
 377         value_set_si(shadow->p[0], 1);
 378         Vector_Copy(upper+2, shadow->p+1, dim+1);
 379         subs[0] = affine2evalue(shadow->p, shadow->p[0], dim+1);
 380         /* new lower
 381          *              (b i - l') + b u' >= 0
 382          * new upper
 383          *              -i >= 0
 384          */
 385         value_assign(shadow->p[1], lower[1]);
 386         value_set_si(a, 1);
 387         value_assign(b, lower[1]);
 388         Vector_Combine(upper+2, lower+2, shadow->p+2,
 389                        b, a, dim+1);
 390         upper = shadow->p+2+dim+1;
 391         lower = shadow->p;
 392         value_set_si(upper[1], -1);
 393         Vector_Set(upper+2, 0, dim+1);
 394     }
 395     value_clear(a);
 396     value_clear(b);
 397
 398     t = evalue_dup(data->e);
 399     evalue_substitute(t, subs);
 400     reduce_evalue(t);
 401     data->e = t;
 402     for (i = 0; i < dim+1; ++i)
 403         evalue_free(subs[i]);
 404     free(subs);
 405
 406     Bernoulli_cb(M, lower, upper, data);
 407
 408     evalue_free(t);
 409     data->e = e;
 410     Vector_Free(shadow);
 411 }
 412
 413 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
 414 {
 415     struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
 416     struct evalue_section_array *sections = data->sections;
 417     Matrix *M2;
 418     Polyhedron *T;
 419     const evalue *factor = NULL;
 420     evalue *linear = NULL;
 421     int constant = 0;
 422     Value tmp;
 423     unsigned dim = M->NbColumns-2;
 424     Vector *row;
 425     int exact = data->options->approximation_method == BV_APPROX_NONE;
 426     int cases = exact ? 3 : 5;
 427
 428     assert(lower);
 429     assert(upper);
 430     evalue_section_array_ensure(sections, sections->ns + cases);
 431
 432     M2 = Matrix_Copy(M);
 433     T = Constraints2Polyhedron(M2, data->options->MaxRays);
 434     Matrix_Free(M2);
 435
 436     POL_ENSURE_VERTICES(T);
 437     if (emptyQ(T)) {
 438         Polyhedron_Free(T);
 439         return;
 440     }
 441
 442     constant = value_notzero_p(data->e->d) ||
 443                 data->e->x.p->type != polynomial ||
 444                 data->e->x.p->pos != 1;
 445     if (!constant && (value_one_p(lower[1]) || value_mone_p(upper[1]))) {
 446         int single_case;
 447         int lower_cst, upper_cst;
 448
 449         lower_cst = First_Non_Zero(lower+2, dim) == -1;
 450         upper_cst = First_Non_Zero(upper+2, dim) == -1;
 451         single_case =
 452             (lower_cst && value_negz_p(lower[2+dim])) ||
 453             (upper_cst && value_negz_p(upper[2+dim])) ||
 454             (lower_cst && upper_cst &&
 455              value_posz_p(lower[2+dim]) && value_posz_p(upper[2+dim]));
 456         if (!single_case) {
 457             transform_to_single_case(M, lower, upper, data);
 458             Polyhedron_Free(T);
 459             return;
 460         }
 461     }
 462
 463     assert(lower != upper);
 464
 465     row = Vector_Alloc(dim+1);
 466     value_init(tmp);
 467     if (value_notzero_p(data->e->d)) {
 468         factor = data->e;
 469         constant = 1;
 470     } else {
 471         if (data->e->x.p->type == polynomial && data->e->x.p->pos == 1)
 472             factor = shifted_copy(&data->e->x.p->arr[0]);
 473         else {
 474             factor = shifted_copy(data->e);
 475             constant = 1;
 476         }
 477     }
 478     linear = linear_term(factor, lower, upper, row, tmp, exact);
 479
 480     if (constant) {
 481         evalue_section_array_add(sections, T, linear);
 482         data->options->stats->bernoulli_sums++;
 483     } else {
 484         evalue *poly_u = NULL, *poly_l = NULL;
 485         Polyhedron *D;
 486         struct poly_list *faulhaber;
 487         assert(data->e->x.p->type == polynomial);
 488         assert(data->e->x.p->pos == 1);
 489         faulhaber = faulhaber_compute(data->e->x.p->size-1);
 490         /* lower is the constraint
 491          *                       b i - l' >= 0          i >= l'/b = l
 492          * upper is the constraint
 493          *                      -a i + u' >= 0          i <= u'/a = u
 494          */
 495         M2 = Matrix_Alloc(3, 2+T->Dimension);
 496         value_set_si(M2->p[0][0], 1);
 497         value_set_si(M2->p[1][0], 1);
 498         value_set_si(M2->p[2][0], 1);
 499         /* Case 1:
 500          *              1 <= l          l' - b >= 0
 501          *              0 < l           l' - 1 >= 0     if exact
 502          */
 503         if (exact)
 504             bound_constraint(M2->p[0]+1, T->Dimension, lower+1, -1, 0, 1);
 505         else
 506             bound_constraint(M2->p[0]+1, T->Dimension, lower+1, -1, -1, 0);
 507         D = AddConstraints(M2->p_Init, 1, T, data->options->MaxRays);
 508         POL_ENSURE_VERTICES(D);
 509         if (emptyQ2(D))
 510             Polyhedron_Free(D);
 511         else {
 512             evalue *extra;
 513             poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
 514                                         faulhaber, data);
 515             Vector_Oppose(lower+2, row->p, dim+1);
 516             extra = power_sums(faulhaber, data->e, row, lower[1], 1, 0, exact);
 517             eadd(poly_u, extra);
 518             eadd(linear, extra);
 519
 520             evalue_section_array_add(sections, D, extra);
 521             data->options->stats->bernoulli_sums++;
 522         }
 523
 524         /* Case 2:
 525          *              1 <= -u         -u' - a >= 0
 526          *              0 < -u          -u' - 1 >= 0    if exact
 527          */
 528         if (exact)
 529             bound_constraint(M2->p[0]+1, T->Dimension, upper+1, -1, 0, 1);
 530         else
 531             bound_constraint(M2->p[0]+1, T->Dimension, upper+1, -1, 1, 0);
 532         D = AddConstraints(M2->p_Init, 1, T, data->options->MaxRays);
 533         POL_ENSURE_VERTICES(D);
 534         if (emptyQ2(D))
 535             Polyhedron_Free(D);
 536         else {
 537             evalue *extra;
 538             poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
 539             Vector_Oppose(upper+2, row->p, dim+1);
 540             value_oppose(tmp, upper[1]);
 541             extra = power_sums(faulhaber, data->e, row, tmp, 1, 1, exact);
 542             eadd(poly_l, extra);
 543             eadd(linear, extra);
 544
 545             evalue_section_array_add(sections, D, extra);
 546             data->options->stats->bernoulli_sums++;
 547         }
 548
 549         /* Case 3:
 550          *              u >= 0          u' >= 0
 551          *              -l >= 0         -l' >= 0
 552          */
 553         bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, 0, 0);
 554         bound_constraint(M2->p[1]+1, T->Dimension, lower+1, 1, 0, 0);
 555         D = AddConstraints(M2->p_Init, 2, T, data->options->MaxRays);
 556         POL_ENSURE_VERTICES(D);
 557         if (emptyQ2(D))
 558             Polyhedron_Free(D);
 559         else {
 560             evalue *e;
 561             poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
 562             poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
 563                                         faulhaber, data);
 564
 565             e = ALLOC(evalue);
 566             value_init(e->d);
 567             evalue_copy(e, poly_u);
 568             eadd(poly_l, e);
 569             eadd(linear, e);
 570             evalue_section_array_add(sections, D, e);
 571             data->options->stats->bernoulli_sums++;
 572         }
 573
 574         if (!exact) {
 575             /* Case 4:
 576              *          l < 1           -l' + b - 1 >= 0
 577              *          0 < l           l' - 1 >= 0
 578              *          u >= 1          u' - a >= 0
 579              */
 580             bound_constraint(M2->p[0]+1, T->Dimension, lower+1, 1, 1, 1);
 581             bound_constraint(M2->p[1]+1, T->Dimension, lower+1, -1, 0, 1);
 582             bound_constraint(M2->p[2]+1, T->Dimension, upper+1, 1, 1, 0);
 583             D = AddConstraints(M2->p_Init, 3, T, data->options->MaxRays);
 584             POL_ENSURE_VERTICES(D);
 585             if (emptyQ2(D))
 586                 Polyhedron_Free(D);
 587             else {
 588                 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
 589                                             faulhaber, data);
 590                 eadd(linear, poly_u);
 591                 evalue_section_array_add(sections, D, poly_u);
 592                 poly_u = NULL;
 593                 data->options->stats->bernoulli_sums++;
 594             }
 595
 596             /* Case 5:
 597              *          -u < 1          u' + a - 1 >= 0
 598              *          0 < -u          -u' - 1 >= 0
 599              *          l <= -1         -l' - b >= 0
 600              */
 601             bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, -1, 1);
 602             bound_constraint(M2->p[1]+1, T->Dimension, upper+1, -1, 0, 1);
 603             bound_constraint(M2->p[2]+1, T->Dimension, lower+1, 1, -1, 0);
 604             D = AddConstraints(M2->p_Init, 3, T, data->options->MaxRays);
 605             POL_ENSURE_VERTICES(D);
 606             if (emptyQ2(D))
 607                 Polyhedron_Free(D);
 608             else {
 609                 poly_l = compute_poly_l(poly_l, lower, row, dim,
 610                                         faulhaber, data);
 611                 eadd(linear, poly_l);
 612                 evalue_section_array_add(sections, D, poly_l);
 613                 poly_l = NULL;
 614                 data->options->stats->bernoulli_sums++;
 615             }
 616         }
 617
 618         Matrix_Free(M2);
 619         Polyhedron_Free(T);
 620         if (poly_l)
 621             evalue_free(poly_l);
 622         if (poly_u)
 623             evalue_free(poly_u);
 624         evalue_free(linear);
 625     }
 626     if (factor != data->e)
 627         evalue_free((evalue *)factor);
 628     value_clear(tmp);
 629     Vector_Free(row);
 630 }
 631
 632 /*
 633  * Move the variable at position pos to the front by exchanging
 634  * that variable with the first variable.
 635  */
 636 static void more_var_to_front(Polyhedron **P_p , evalue **E_p, int pos)
 637 {
 638     Polyhedron *P = *P_p;
 639     evalue *E = *E_p;
 640     unsigned dim = P->Dimension;
 641
 642     assert(pos != 0);
 643
 644     P = Polyhedron_Copy(P);
 645     Polyhedron_ExchangeColumns(P, 1, 1+pos);
 646     *P_p = P;
 647
 648     if (value_zero_p(E->d)) {
 649         int j;
 650         evalue **subs;
 651
 652         subs = ALLOCN(evalue *, dim);
 653         for (j = 0; j < dim; ++j)
 654             subs[j] = evalue_var(j);
 655         E = subs[0];
 656         subs[0] = subs[pos];
 657         subs[pos] = E;
 658         E = evalue_dup(*E_p);
 659         evalue_substitute(E, subs);
 660         for (j = 0; j < dim; ++j)
 661             evalue_free(subs[j]);
 662         free(subs);
 663         *E_p = E;
 664     }
 665 }
 666
 667 static int type_offset(enode *p)
 668 {
 669    return p->type == fractional ? 1 :
 670           p->type == flooring ? 1 : 0;
 671 }
 672
 673 static void adjust_periods(evalue *e, unsigned nvar, Vector *periods)
 674 {
 675     for (; value_zero_p(e->d); e = &e->x.p->arr[0]) {
 676         int pos;
 677         assert(e->x.p->type == polynomial);
 678         assert(e->x.p->size == 2);
 679         assert(value_notzero_p(e->x.p->arr[1].d));
 680
 681         pos = e->x.p->pos - 1;
 682         if (pos >= nvar)
 683             break;
 684
 685         value_lcm(periods->p[pos], periods->p[pos], e->x.p->arr[1].d);
 686     }
 687 }
 688
 689 static void fractional_periods_r(evalue *e, unsigned nvar, Vector *periods)
 690 {
 691     int i;
 692
 693     if (value_notzero_p(e->d))
 694         return;
 695
 696     assert(e->x.p->type == polynomial || e->x.p->type == fractional);
 697
 698     if (e->x.p->type == fractional)
 699         adjust_periods(&e->x.p->arr[0], nvar, periods);
 700
 701     for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
 702         fractional_periods_r(&e->x.p->arr[i], nvar, periods);
 703 }
 704
 705 /*
 706  * For each of the nvar variables, compute the lcm of the
 707  * denominators of the coefficients of that variable in
 708  * any of the fractional parts.
 709  */
 710 static Vector *fractional_periods(evalue *e, unsigned nvar)
 711 {
 712     int i;
 713     Vector *periods = Vector_Alloc(nvar);
 714
 715     for (i = 0; i < nvar; ++i)
 716         value_set_si(periods->p[i], 1);
 717
 718     fractional_periods_r(e, nvar, periods);
 719
 720     return periods;
 721 }
 722
 723 /* Move "best" variable to sum over into the first position,
 724  * possibly changing *P_p and *E_p.
 725  *
 726  * If there are any fractional parts (period != NULL), then we
 727  * first look for a variable with the smallest lcm of denominators
 728  * inside a fractional part.  This denominator is assigned to period
 729  * and corresponds to the number of "splinters" we need to construct
 730  * at this level.
 731  *
 732  * Of those with this denominator (all if period == NULL or if there
 733  * are no fractional parts), we select the variable with the smallest
 734  * maximal coefficient, as this coefficient will become a denominator
 735  * in new fractional parts (for an exact computation), which may
 736  * lead to splintering in the next step.
 737  */
 738 static void move_best_to_front(Polyhedron **P_p, evalue **E_p, unsigned nvar,
 739                                Value *period)
 740 {
 741     Polyhedron *P = *P_p;
 742     evalue *E = *E_p;
 743     int i, j, best_i = -1;
 744     Vector *periods;
 745     Value best, max;
 746
 747     assert(nvar >= 1);
 748
 749     if (period) {
 750         periods = fractional_periods(E, nvar);
 751         value_assign(*period, periods->p[0]);
 752         for (i = 1; i < nvar; ++i)
 753             if (value_lt(periods->p[i], *period))
 754                 value_assign(*period, periods->p[i]);
 755     }
 756
 757     value_init(best);
 758     value_init(max);
 759
 760     for (i = 0; i < nvar; ++i) {
 761         if (period && value_ne(*period, periods->p[i]))
 762             continue;
 763
 764         value_set_si(max, 0);
 765
 766         for (j = 0; j < P->NbConstraints; ++j) {
 767             if (value_zero_p(P->Constraint[j][1+i]))
 768                 continue;
 769             if (First_Non_Zero(P->Constraint[j]+1, i) == -1 &&
 770                 First_Non_Zero(P->Constraint[j]+1+i+1, nvar-i-1) == -1)
 771                 continue;
 772             if (value_abs_gt(P->Constraint[j][1+i], max))
 773                 value_absolute(max, P->Constraint[j][1+i]);
 774         }
 775
 776         if (best_i == -1 || value_lt(max, best)) {
 777             value_assign(best, max);
 778             best_i = i;
 779         }
 780     }
 781
 782     value_clear(best);
 783     value_clear(max);
 784
 785     if (period)
 786         Vector_Free(periods);
 787
 788     if (best_i != 0)
 789         more_var_to_front(P_p, E_p, best_i);
 790
 791     return;
 792 }
 793
 794 static evalue *sum_over_polytope_base(Polyhedron *P, evalue *E, unsigned nvar,
 795                                       struct evalue_section_array *sections,
 796                                       struct barvinok_options *options)
 797 {
 798     evalue *res;
 799     struct Bernoulli_data data;
 800
 801     assert(P->NbEq == 0);
 802
 803     sections->ns = 0;
 804     data.options = options;
 805     data.sections = sections;
 806     data.e = E;
 807
 808     for_each_lower_upper_bound(P, Bernoulli_init, Bernoulli_cb, &data);
 809
 810     res = evalue_from_section_array(sections->s, sections->ns);
 811
 812     if (nvar > 1) {
 813         evalue *tmp = barvinok_summate(res, nvar-1, options);
 814         evalue_free(res);
 815         res = tmp;
 816     }
 817
 818     return res;
 819 }
 820
 821 static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E,
 822                                  unsigned nvar,
 823                                  struct evalue_section_array *sections,
 824                                  struct barvinok_options *options)
 825 {
 826     unsigned dim = P->Dimension;
 827     unsigned new_dim, new_nparam;
 828     Matrix *T = NULL, *CP = NULL;
 829     evalue **subs;
 830     evalue *sum;
 831     int j;
 832
 833     if (emptyQ(P))
 834         return evalue_zero();
 835
 836     assert(P->NbEq > 0);
 837
 838     remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays);
 839
 840     if (emptyQ(P)) {
 841         Polyhedron_Free(P);
 842         return evalue_zero();
 843     }
 844
 845     new_nparam = CP ? CP->NbColumns-1 : dim - nvar;
 846     new_dim = T ? T->NbColumns-1 : nvar + new_nparam;
 847
 848     /* We can avoid these substitutions if E is a constant */
 849     subs = ALLOCN(evalue *, dim);
 850     for (j = 0; j < nvar; ++j) {
 851         if (T)
 852             subs[j] = affine2evalue(T->p[j], T->p[nvar+new_nparam][new_dim],
 853                                     new_dim);
 854         else
 855             subs[j] = evalue_var(j);
 856     }
 857     for (j = 0; j < dim-nvar; ++j) {
 858         if (CP)
 859             subs[nvar+j] = affine2evalue(CP->p[j], CP->p[dim-nvar][new_nparam],
 860                                          new_nparam);
 861         else
 862             subs[nvar+j] = evalue_var(j);
 863         evalue_shift_variables(subs[nvar+j], new_dim-new_nparam);
 864     }
 865
 866     E = evalue_dup(E);
 867     evalue_substitute(E, subs);
 868     reduce_evalue(E);
 869
 870     for (j = 0; j < dim; ++j)
 871         evalue_free(subs[j]);
 872     free(subs);
 873
 874     if (new_dim-new_nparam > 0) {
 875         sum = bernoulli_summate(P, E, new_dim-new_nparam, sections, options);
 876         evalue_free(E);
 877         Polyhedron_Free(P);
 878     } else {
 879         sum = ALLOC(evalue);
 880         value_init(sum->d);
 881         sum->x.p = new_enode(partition, 2, new_dim);
 882         EVALUE_SET_DOMAIN(sum->x.p->arr[0], P);
 883         value_clear(sum->x.p->arr[1].d);
 884         sum->x.p->arr[1] = *E;
 885         free(E);
 886     }
 887
 888     if (CP) {
 889         evalue_backsubstitute(sum, CP, options->MaxRays);
 890         Matrix_Free(CP);
 891     }
 892
 893     if (T)
 894         Matrix_Free(T);
 895
 896     return sum;
 897 }
 898
 899 /* Splinter P into period slices along the first variable x = period y + c,
 900  * 0 <= c < perdiod, * ensuring no fractional part contains the first variable,
 901  * and sum over all slices.
 902  */
 903 static evalue *sum_over_polytope_slices(Polyhedron *P, evalue *E,
 904                                         unsigned nvar,
 905                                         Value period,
 906                                         struct evalue_section_array *sections,
 907                                         struct barvinok_options *options)
 908 {
 909     evalue *sum = evalue_zero();
 910     evalue **subs;
 911     unsigned dim = P->Dimension;
 912     Matrix *T;
 913     Value i;
 914     Value one;
 915     int j;
 916
 917     value_init(i);
 918     value_init(one);
 919     value_set_si(one, 1);
 920
 921     subs = ALLOCN(evalue *, dim);
 922
 923     T = Matrix_Alloc(dim+1, dim+1);
 924     value_assign(T->p[0][0], period);
 925     for (j = 1; j < dim+1; ++j)
 926         value_set_si(T->p[j][j], 1);
 927
 928     for (j = 0; j < dim; ++j)
 929         subs[j] = evalue_var(j);
 930     evalue_mul(subs[0], period);
 931
 932     for (value_set_si(i, 0); value_lt(i, period); value_increment(i, i)) {
 933         evalue *tmp;
 934         Polyhedron *S = DomainPreimage(P, T, options->MaxRays);
 935         evalue *e = evalue_dup(E);
 936         evalue_substitute(e, subs);
 937         reduce_evalue(e);
 938
 939         if (S->NbEq)
 940             tmp = sum_over_polytope_with_equalities(S, e, nvar, sections, options);
 941         else
 942             tmp = sum_over_polytope_base(S, e, nvar, sections, options);
 943
 944         assert(tmp);
 945         eadd(tmp, sum);
 946         evalue_free(tmp);
 947
 948         Domain_Free(S);
 949         evalue_free(e);
 950
 951         value_increment(T->p[0][dim], T->p[0][dim]);
 952         evalue_add_constant(subs[0], one);
 953     }
 954
 955     value_clear(i);
 956     value_clear(one);
 957     Matrix_Free(T);
 958     for (j = 0; j < dim; ++j)
 959         evalue_free(subs[j]);
 960     free(subs);
 961
 962     reduce_evalue(sum);
 963     return sum;
 964 }
 965
 966 evalue *bernoulli_summate(Polyhedron *P, evalue *E, unsigned nvar,
 967                                  struct evalue_section_array *sections,
 968                                  struct barvinok_options *options)
 969 {
 970     Polyhedron *P_orig = P;
 971     evalue *E_orig = E;
 972     Value period;
 973     evalue *sum;
 974     int exact = options->approximation_method == BV_APPROX_NONE;
 975
 976     if (P->NbEq)
 977         return sum_over_polytope_with_equalities(P, E, nvar, sections, options);
 978
 979     value_init(period);
 980
 981     move_best_to_front(&P, &E, nvar, exact ? &period : NULL);
 982     if (exact && value_notone_p(period))
 983         sum = sum_over_polytope_slices(P, E, nvar, period, sections, options);
 984     else
 985         sum = sum_over_polytope_base(P, E, nvar, sections, options);
 986
 987     if (P != P_orig)
 988         Polyhedron_Free(P);
 989     if (E != E_orig)
 990         evalue_free(E);
 991
 992     value_clear(period);
 993
 994     return sum;
 995 }