bernoulli.c

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