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