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