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(data->ns < data->size);
 219
 220     M2 = Matrix_Copy(M);
 221     T = Constraints2Polyhedron(M2, data->MaxRays);
 222     Matrix_Free(M2);
 223
 224     POL_ENSURE_VERTICES(T);
 225     if (emptyQ(T)) {
 226         Polyhedron_Free(T);
 227         return;
 228     }
 229
 230     assert(lower != upper);
 231
 232     row = Vector_Alloc(dim+1);
 233     value_init(tmp);
 234     if (value_notzero_p(data->e->d)) {
 235         factor = data->e;
 236         constant = 1;
 237     } else {
 238         assert(data->e->x.p->type == polynomial);
 239         if (data->e->x.p->pos > 1) {
 240             factor = shifted_copy(data->e);
 241             constant = 1;
 242         } else
 243             factor = shifted_copy(&data->e->x.p->arr[0]);
 244     }
 245     if (!EVALUE_IS_ZERO(*factor)) {
 246         value_absolute(tmp, upper[1]);
 247         /* upper - lower */
 248         Vector_Combine(lower+2, upper+2, row->p, tmp, lower[1], dim+1);
 249         value_multiply(tmp, tmp, lower[1]);
 250         /* upper - lower + 1 */
 251         value_addto(row->p[dim], row->p[dim], tmp);
 252
 253         linear = affine2evalue(row->p, tmp, dim);
 254         emul(factor, linear);
 255     } else
 256         linear = evalue_zero();
 257
 258     if (constant) {
 259         data->s[data->ns].E = linear;
 260         data->s[data->ns].D = T;
 261         ++data->ns;
 262     } else {
 263         evalue *poly_u = NULL, *poly_l = NULL;
 264         Polyhedron *D;
 265         struct poly_list *faulhaber;
 266         assert(data->e->x.p->type == polynomial);
 267         assert(data->e->x.p->pos == 1);
 268         faulhaber = faulhaber_compute(data->e->x.p->size-1);
 269         /* lower is the constraint
 270          *                       b i - l' >= 0          i >= l'/b = l
 271          * upper is the constraint
 272          *                      -a i + u' >= 0          i <= u'/a = u
 273          */
 274         M2 = Matrix_Alloc(3, 2+T->Dimension);
 275         value_set_si(M2->p[0][0], 1);
 276         value_set_si(M2->p[1][0], 1);
 277         value_set_si(M2->p[2][0], 1);
 278         /* Case 1:
 279          *              1 <= l          l' - b >= 0
 280          */
 281         Vector_Oppose(lower+2, M2->p[0]+1, T->Dimension+1);
 282         value_subtract(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension],
 283                       lower[1]);
 284         D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
 285         if (emptyQ2(D))
 286             Polyhedron_Free(D);
 287         else {
 288             evalue *extra;
 289             if (!poly_u) {
 290                 Vector_Copy(upper+2, row->p, dim+1);
 291                 value_oppose(tmp, upper[1]);
 292                 value_addto(row->p[dim], row->p[dim], tmp);
 293                 poly_u = power_sums(faulhaber, data->e, row, tmp, 0, 0);
 294             }
 295             Vector_Oppose(lower+2, row->p, dim+1);
 296             extra = power_sums(faulhaber, data->e, row, lower[1], 1, 0);
 297             eadd(poly_u, extra);
 298             eadd(linear, extra);
 299
 300             data->s[data->ns].E = extra;
 301             data->s[data->ns].D = D;
 302             ++data->ns;
 303         }
 304
 305         /* Case 2:
 306          *              1 <= -u         -u' - a >= 0
 307          */
 308         Vector_Oppose(upper+2, M2->p[0]+1, T->Dimension+1);
 309         value_addto(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension],
 310                       upper[1]);
 311         D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
 312         if (emptyQ2(D))
 313             Polyhedron_Free(D);
 314         else {
 315             evalue *extra;
 316             if (!poly_l) {
 317                 Vector_Copy(lower+2, row->p, dim+1);
 318                 value_addto(row->p[dim], row->p[dim], lower[1]);
 319                 poly_l = power_sums(faulhaber, data->e, row, lower[1], 0, 1);
 320             }
 321             Vector_Oppose(upper+2, row->p, dim+1);
 322             value_oppose(tmp, upper[1]);
 323             extra = power_sums(faulhaber, data->e, row, tmp, 1, 1);
 324             eadd(poly_l, extra);
 325             eadd(linear, extra);
 326
 327             data->s[data->ns].E = extra;
 328             data->s[data->ns].D = D;
 329             ++data->ns;
 330         }
 331
 332         /* Case 3:
 333          *              u >= 0          u' >= 0
 334          *              -l >= 0         -l' >= 0
 335          */
 336         Vector_Copy(upper+2, M2->p[0]+1, T->Dimension+1);
 337         Vector_Copy(lower+2, M2->p[1]+1, T->Dimension+1);
 338         D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
 339         if (emptyQ2(D))
 340             Polyhedron_Free(D);
 341         else {
 342             if (!poly_l) {
 343                 Vector_Copy(lower+2, row->p, dim+1);
 344                 value_addto(row->p[dim], row->p[dim], lower[1]);
 345                 poly_l = power_sums(faulhaber, data->e, row, lower[1], 0, 1);
 346             }
 347             if (!poly_u) {
 348                 Vector_Copy(upper+2, row->p, dim+1);
 349                 value_oppose(tmp, upper[1]);
 350                 value_addto(row->p[dim], row->p[dim], tmp);
 351                 poly_u = power_sums(faulhaber, data->e, row, tmp, 0, 0);
 352             }
 353
 354             data->s[data->ns].E = ALLOC(evalue);
 355             value_init(data->s[data->ns].E->d);
 356             evalue_copy(data->s[data->ns].E, poly_u);
 357             eadd(poly_l, data->s[data->ns].E);
 358             eadd(linear, data->s[data->ns].E);
 359             data->s[data->ns].D = D;
 360             ++data->ns;
 361         }
 362
 363         /* Case 4:
 364          *              l < 1           -l' + b - 1 >= 0
 365          *              0 < l           l' - 1 >= 0
 366          */
 367         Vector_Copy(lower+2, M2->p[0]+1, T->Dimension+1);
 368         value_addto(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension], lower[1]);
 369         value_decrement(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension]);
 370         Vector_Oppose(lower+2, M2->p[1]+1, T->Dimension+1);
 371         value_decrement(M2->p[1][1+T->Dimension], M2->p[1][1+T->Dimension]);
 372         D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
 373         if (emptyQ2(D))
 374             Polyhedron_Free(D);
 375         else {
 376             if (!poly_u) {
 377                 Vector_Copy(upper+2, row->p, dim+1);
 378                 value_oppose(tmp, upper[1]);
 379                 value_addto(row->p[dim], row->p[dim], tmp);
 380                 poly_u = power_sums(faulhaber, data->e, row, tmp, 0, 0);
 381             }
 382
 383             eadd(linear, poly_u);
 384             data->s[data->ns].E = poly_u;
 385             poly_u = NULL;
 386             data->s[data->ns].D = D;
 387             ++data->ns;
 388         }
 389
 390         /* Case 5:
 391          *              -u < 1          u' + a - 1 >= 0
 392          *              0 < -u          -u' - 1 >= 0
 393          *              l <= 0          -l' >= 0
 394          */
 395         Vector_Copy(upper+2, M2->p[0]+1, T->Dimension+1);
 396         value_subtract(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension],
 397                         upper[1]);
 398         value_decrement(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension]);
 399         Vector_Oppose(upper+2, M2->p[1]+1, T->Dimension+1);
 400         value_decrement(M2->p[1][1+T->Dimension], M2->p[1][1+T->Dimension]);
 401         Vector_Copy(lower+2, M2->p[2]+1, T->Dimension+1);
 402         D = AddConstraints(M2->p_Init, 3, T, data->MaxRays);
 403         if (emptyQ2(D))
 404             Polyhedron_Free(D);
 405         else {
 406             if (!poly_l) {
 407                 Vector_Copy(lower+2, row->p, dim+1);
 408                 value_addto(row->p[dim], row->p[dim], lower[1]);
 409                 poly_l = power_sums(faulhaber, data->e, row, lower[1], 0, 1);
 410             }
 411
 412             eadd(linear, poly_l);
 413             data->s[data->ns].E = poly_l;
 414             poly_l = NULL;
 415             data->s[data->ns].D = D;
 416             ++data->ns;
 417         }
 418
 419         Matrix_Free(M2);
 420         Polyhedron_Free(T);
 421         if (poly_l)
 422             evalue_free(poly_l);
 423         if (poly_u)
 424             evalue_free(poly_u);
 425         evalue_free(linear);
 426     }
 427     if (factor != data->e)
 428         evalue_free(factor);
 429     value_clear(tmp);
 430     Vector_Free(row);
 431 }
 432
 433 /* Looks for variable with integer bounds, i.e., with coefficients 0, 1 or -1.
 434  * Returns 1 if such a variable is found and puts it in the first position,
 435  * possibly changing *P_p and *E_p.
 436  */
 437 static int find_integer_bounds(Polyhedron **P_p, evalue **E_p, unsigned nvar)
 438 {
 439     Polyhedron *P = *P_p;
 440     evalue *E = *E_p;
 441     unsigned dim = P->Dimension;
 442     int i, j;
 443
 444     for (i = 0; i < nvar; ++i) {
 445         for (j = 0; j < P->NbConstraints; ++j) {
 446             if (value_zero_p(P->Constraint[j][1+i]))
 447                 continue;
 448             if (value_one_p(P->Constraint[j][1+i]))
 449                 continue;
 450             if (value_mone_p(P->Constraint[j][1+i]))
 451                 continue;
 452             break;
 453         }
 454         if (j == P->NbConstraints)
 455             break;
 456     }
 457     if (i == nvar)
 458         return 0;
 459     if (i == 0)
 460         return 1;
 461     P = Polyhedron_Copy(P);
 462     Polyhedron_ExchangeColumns(P, 1, 1+i);
 463     *P_p = P;
 464
 465     if (value_zero_p(E->d)) {
 466         evalue **subs;
 467         subs = ALLOCN(evalue *, dim);
 468         for (j = 0; j < dim; ++j)
 469             subs[j] = evalue_var(j);
 470         E = subs[0];
 471         subs[0] = subs[i];
 472         subs[i] = E;
 473         E = evalue_dup(*E_p);
 474         evalue_substitute(E, subs);
 475         for (j = 0; j < dim; ++j)
 476             evalue_free(subs[j]);
 477         free(subs);
 478         *E_p = E;
 479     }
 480
 481     return 1;
 482 }
 483
 484 static evalue *sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
 485                                  struct Bernoulli_data *data,
 486                                  struct barvinok_options *options)
 487 {
 488     unsigned dim = P->Dimension - 1;
 489     evalue *res;
 490
 491     if (value_zero_p(P->Constraint[0][0]) &&
 492             value_notzero_p(P->Constraint[0][1])) {
 493         res = ALLOC(evalue);
 494         value_init(res->d);
 495         value_set_si(res->d, 0);
 496         res->x.p = new_enode(partition, 2, dim);
 497         EVALUE_SET_DOMAIN(res->x.p->arr[0], Polyhedron_Project(P, dim));
 498         evalue_copy(&res->x.p->arr[1], E);
 499         reduce_evalue_in_domain(&res->x.p->arr[1], P);
 500         shift(&res->x.p->arr[1]);
 501     } else {
 502         data->ns = 0;
 503         data->e = E;
 504
 505         for_each_lower_upper_bound(P, Bernoulli_cb, data);
 506
 507         res = evalue_from_section_array(data->s, data->ns);
 508     }
 509
 510     if (nvar > 1) {
 511         evalue *tmp = Bernoulli_sum_evalue(res, nvar-1, options);
 512         evalue_free(res);
 513         res = tmp;
 514     }
 515
 516     return res;
 517 }
 518
 519 evalue *Bernoulli_sum_evalue(evalue *e, unsigned nvar,
 520                              struct barvinok_options *options)
 521 {
 522     struct Bernoulli_data data;
 523     int i, j;
 524     evalue *sum = evalue_zero();
 525
 526     if (EVALUE_IS_ZERO(*e))
 527         return sum;
 528
 529     if (nvar == 0) {
 530         eadd(e, sum);
 531         return sum;
 532     }
 533
 534     assert(value_zero_p(e->d));
 535     assert(e->x.p->type == partition);
 536
 537     data.size = e->x.p->size * 2 + 16;
 538     data.s = ALLOCN(struct evalue_section, data.size);
 539     data.MaxRays = options->MaxRays;
 540
 541     for (i = 0; i < e->x.p->size/2; ++i) {
 542         Polyhedron *D;
 543         for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
 544             evalue *E = &e->x.p->arr[2*i+1];
 545             Polyhedron *P = D;
 546             Polyhedron *next = D->next;
 547             evalue *tmp;
 548             int integer_bounds;
 549
 550             P->next = NULL;
 551
 552             integer_bounds = find_integer_bounds(&P, &E, nvar);
 553             if (options->approximation_method == BV_APPROX_NONE &&
 554                 !integer_bounds) {
 555                 evalue_free(sum);
 556                 sum = NULL;
 557             } else {
 558                 evalue *tmp = sum_over_polytope(P, E, nvar, &data, options);
 559                 eadd(tmp, sum);
 560                 evalue_free(tmp);
 561             }
 562
 563             if (P != D)
 564                 Polyhedron_Free(P);
 565             if (E != &e->x.p->arr[2*i+1])
 566                 evalue_free(E);
 567
 568             D->next = next;;
 569
 570             if (!sum)
 571                 break;
 572         }
 573
 574         if (!sum)
 575             break;
 576     }
 577
 578     free(data.s);
 579
 580     if (sum)
 581         reduce_evalue(sum);
 582     return sum;
 583 }
 584
 585 evalue *Bernoulli_sum(Polyhedron *P, Polyhedron *C,
 586                         struct barvinok_options *options)
 587 {
 588     evalue e;
 589     evalue *sum;
 590
 591     value_init(e.d);
 592     e.x.p = new_enode(partition, 2, P->Dimension);
 593     EVALUE_SET_DOMAIN(e.x.p->arr[0], Polyhedron_Copy(P));
 594     evalue_set_si(&e.x.p->arr[1], 1, 1);
 595     sum = Bernoulli_sum_evalue(&e, P->Dimension - C->Dimension, options);
 596     free_evalue_refs(&e);
 597     return sum;
 598 }