bernoulli.c

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