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         Gcd(bernoulli_coef.num->p[i], bernoulli_coef.den->p[i], &tmp);
  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-1], bernoulli_coef.den->p[i],
  75                   &bernoulli_coef.lcm->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         free_evalue_refs(factor);
 189         free_evalue_refs(term);
 190         free(factor);
 191         free(term);
 192     }
 193     if (neg)
 194         evalue_negate(sum);
 195     free_evalue_refs(base);
 196     free(base);
 197
 198     return sum;
 199 }
 200
 201 struct section { Polyhedron *D; evalue *E; };
 202
 203 struct Bernoulli_data {
 204     unsigned MaxRays;
 205     struct section *s;
 206     int size;
 207     int ns;
 208     evalue *e;
 209 };
 210
 211 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
 212 {
 213     struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
 214     Matrix *M2;
 215     Polyhedron *T;
 216     evalue *factor = NULL;
 217     evalue *linear = NULL;
 218     int constant = 0;
 219     Value tmp;
 220     unsigned dim = M->NbColumns-2;
 221     Vector *row;
 222
 223     assert(data->ns < data->size);
 224
 225     M2 = Matrix_Copy(M);
 226     T = Constraints2Polyhedron(M2, data->MaxRays);
 227     Matrix_Free(M2);
 228
 229     POL_ENSURE_VERTICES(T);
 230     if (emptyQ(T)) {
 231         Polyhedron_Free(T);
 232         return;
 233     }
 234
 235     assert(lower != upper);
 236
 237     row = Vector_Alloc(dim+1);
 238     value_init(tmp);
 239     if (value_notzero_p(data->e->d)) {
 240         factor = data->e;
 241         constant = 1;
 242     } else {
 243         assert(data->e->x.p->type == polynomial);
 244         if (data->e->x.p->pos > 1) {
 245             factor = shifted_copy(data->e);
 246             constant = 1;
 247         } else
 248             factor = shifted_copy(&data->e->x.p->arr[0]);
 249     }
 250     if (!EVALUE_IS_ZERO(*factor)) {
 251         value_absolute(tmp, upper[1]);
 252         /* upper - lower */
 253         Vector_Combine(lower+2, upper+2, row->p, tmp, lower[1], dim+1);
 254         value_multiply(tmp, tmp, lower[1]);
 255         /* upper - lower + 1 */
 256         value_addto(row->p[dim], row->p[dim], tmp);
 257
 258         linear = affine2evalue(row->p, tmp, dim);
 259         emul(factor, linear);
 260     } else
 261         linear = evalue_zero();
 262
 263     if (constant) {
 264         data->s[data->ns].E = linear;
 265         data->s[data->ns].D = T;
 266         ++data->ns;
 267     } else {
 268         evalue *poly_u = NULL, *poly_l = NULL;
 269         Polyhedron *D;
 270         struct poly_list *faulhaber;
 271         assert(data->e->x.p->type == polynomial);
 272         assert(data->e->x.p->pos == 1);
 273         faulhaber = faulhaber_compute(data->e->x.p->size-1);
 274         /* lower is the constraint
 275          *                       b i - l' >= 0          i >= l'/b = l
 276          * upper is the constraint
 277          *                      -a i + u' >= 0          i <= u'/a = u
 278          */
 279         M2 = Matrix_Alloc(3, 2+T->Dimension);
 280         value_set_si(M2->p[0][0], 1);
 281         value_set_si(M2->p[1][0], 1);
 282         value_set_si(M2->p[2][0], 1);
 283         /* Case 1:
 284          *              1 <= l          l' - b >= 0
 285          */
 286         Vector_Oppose(lower+2, M2->p[0]+1, T->Dimension+1);
 287         value_subtract(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension],
 288                       lower[1]);
 289         D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
 290         if (emptyQ2(D))
 291             Polyhedron_Free(D);
 292         else {
 293             evalue *extra;
 294             if (!poly_u) {
 295                 Vector_Copy(upper+2, row->p, dim+1);
 296                 value_oppose(tmp, upper[1]);
 297                 value_addto(row->p[dim], row->p[dim], tmp);
 298                 poly_u = power_sums(faulhaber, data->e, row, tmp, 0, 0);
 299             }
 300             Vector_Oppose(lower+2, row->p, dim+1);
 301             extra = power_sums(faulhaber, data->e, row, lower[1], 1, 0);
 302             eadd(poly_u, extra);
 303             eadd(linear, extra);
 304
 305             data->s[data->ns].E = extra;
 306             data->s[data->ns].D = D;
 307             ++data->ns;
 308         }
 309
 310         /* Case 2:
 311          *              1 <= -u         -u' - a >= 0
 312          */
 313         Vector_Oppose(upper+2, M2->p[0]+1, T->Dimension+1);
 314         value_addto(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension],
 315                       upper[1]);
 316         D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
 317         if (emptyQ2(D))
 318             Polyhedron_Free(D);
 319         else {
 320             evalue *extra;
 321             if (!poly_l) {
 322                 Vector_Copy(lower+2, row->p, dim+1);
 323                 value_addto(row->p[dim], row->p[dim], lower[1]);
 324                 poly_l = power_sums(faulhaber, data->e, row, lower[1], 0, 1);
 325             }
 326             Vector_Oppose(upper+2, row->p, dim+1);
 327             value_oppose(tmp, upper[1]);
 328             extra = power_sums(faulhaber, data->e, row, tmp, 1, 1);
 329             eadd(poly_l, extra);
 330             eadd(linear, extra);
 331
 332             data->s[data->ns].E = extra;
 333             data->s[data->ns].D = D;
 334             ++data->ns;
 335         }
 336
 337         /* Case 3:
 338          *              u >= 0          u' >= 0
 339          *              -l >= 0         -l' >= 0
 340          */
 341         Vector_Copy(upper+2, M2->p[0]+1, T->Dimension+1);
 342         Vector_Copy(lower+2, M2->p[1]+1, T->Dimension+1);
 343         D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
 344         if (emptyQ2(D))
 345             Polyhedron_Free(D);
 346         else {
 347             if (!poly_l) {
 348                 Vector_Copy(lower+2, row->p, dim+1);
 349                 value_addto(row->p[dim], row->p[dim], lower[1]);
 350                 poly_l = power_sums(faulhaber, data->e, row, lower[1], 0, 1);
 351             }
 352             if (!poly_u) {
 353                 Vector_Copy(upper+2, row->p, dim+1);
 354                 value_oppose(tmp, upper[1]);
 355                 value_addto(row->p[dim], row->p[dim], tmp);
 356                 poly_u = power_sums(faulhaber, data->e, row, tmp, 0, 0);
 357             }
 358
 359             data->s[data->ns].E = ALLOC(evalue);
 360             value_init(data->s[data->ns].E->d);
 361             evalue_copy(data->s[data->ns].E, poly_u);
 362             eadd(poly_l, data->s[data->ns].E);
 363             eadd(linear, data->s[data->ns].E);
 364             data->s[data->ns].D = D;
 365             ++data->ns;
 366         }
 367
 368         /* Case 4:
 369          *              l < 1           -l' + b - 1 >= 0
 370          *              0 < l           l' - 1 >= 0
 371          */
 372         Vector_Copy(lower+2, M2->p[0]+1, T->Dimension+1);
 373         value_addto(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension], lower[1]);
 374         value_decrement(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension]);
 375         Vector_Oppose(lower+2, M2->p[1]+1, T->Dimension+1);
 376         value_decrement(M2->p[1][1+T->Dimension], M2->p[1][1+T->Dimension]);
 377         D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
 378         if (emptyQ2(D))
 379             Polyhedron_Free(D);
 380         else {
 381             if (!poly_u) {
 382                 Vector_Copy(upper+2, row->p, dim+1);
 383                 value_oppose(tmp, upper[1]);
 384                 value_addto(row->p[dim], row->p[dim], tmp);
 385                 poly_u = power_sums(faulhaber, data->e, row, tmp, 0, 0);
 386             }
 387
 388             eadd(linear, poly_u);
 389             data->s[data->ns].E = poly_u;
 390             poly_u = NULL;
 391             data->s[data->ns].D = D;
 392             ++data->ns;
 393         }
 394
 395         /* Case 5:
 396          *              -u < 1          u' + a - 1 >= 0
 397          *              0 < -u          -u' - 1 >= 0
 398          *              l <= 0          -l' >= 0
 399          */
 400         Vector_Copy(upper+2, M2->p[0]+1, T->Dimension+1);
 401         value_subtract(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension],
 402                         upper[1]);
 403         value_decrement(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension]);
 404         Vector_Oppose(upper+2, M2->p[1]+1, T->Dimension+1);
 405         value_decrement(M2->p[1][1+T->Dimension], M2->p[1][1+T->Dimension]);
 406         Vector_Copy(lower+2, M2->p[2]+1, T->Dimension+1);
 407         D = AddConstraints(M2->p_Init, 3, T, data->MaxRays);
 408         if (emptyQ2(D))
 409             Polyhedron_Free(D);
 410         else {
 411             if (!poly_l) {
 412                 Vector_Copy(lower+2, row->p, dim+1);
 413                 value_addto(row->p[dim], row->p[dim], lower[1]);
 414                 poly_l = power_sums(faulhaber, data->e, row, lower[1], 0, 1);
 415             }
 416
 417             eadd(linear, poly_l);
 418             data->s[data->ns].E = poly_l;
 419             poly_l = NULL;
 420             data->s[data->ns].D = D;
 421             ++data->ns;
 422         }
 423
 424         Matrix_Free(M2);
 425         Polyhedron_Free(T);
 426         if (poly_l) {
 427             free_evalue_refs(poly_l);
 428             free(poly_l);
 429         }
 430         if (poly_u) {
 431             free_evalue_refs(poly_u);
 432             free(poly_u);
 433         }
 434         free_evalue_refs(linear);
 435         free(linear);
 436     }
 437     if (factor != data->e) {
 438         free_evalue_refs(factor);
 439         free(factor);
 440     }
 441     value_clear(tmp);
 442     Vector_Free(row);
 443 }
 444
 445 evalue *Bernoulli_sum_evalue(evalue *e, unsigned nvar,
 446                              struct barvinok_options *options)
 447 {
 448     struct Bernoulli_data data;
 449     int i, j;
 450     evalue *sum = evalue_zero();
 451
 452     if (EVALUE_IS_ZERO(*e))
 453         return sum;
 454
 455     if (nvar == 0) {
 456         eadd(e, sum);
 457         return sum;
 458     }
 459
 460     assert(value_zero_p(e->d));
 461     assert(e->x.p->type == partition);
 462
 463     data.size = e->x.p->size * 2 + 16;
 464     data.s = ALLOCN(struct section, data.size);
 465     data.MaxRays = options->MaxRays;
 466
 467     for (i = 0; i < e->x.p->size/2; ++i) {
 468         Polyhedron *D;
 469         for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
 470             unsigned dim = D->Dimension - 1;
 471             Polyhedron *next = D->next;
 472             evalue tmp;
 473             D->next = NULL;
 474
 475             value_init(tmp.d);
 476             value_set_si(tmp.d, 0);
 477
 478             if (value_zero_p(D->Constraint[0][0]) &&
 479                     value_notzero_p(D->Constraint[0][1])) {
 480                 tmp.x.p = new_enode(partition, 2, dim);
 481                 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Project(D, dim));
 482                 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
 483                 reduce_evalue_in_domain(&tmp.x.p->arr[1], D);
 484                 shift(&tmp.x.p->arr[1]);
 485             } else {
 486                 data.ns = 0;
 487                 data.e = &e->x.p->arr[2*i+1];
 488
 489                 for_each_lower_upper_bound(D, Bernoulli_cb, &data);
 490
 491                 if (data.ns == 0)
 492                     evalue_set_si(&tmp, 0, 1);
 493                 else {
 494                     tmp.x.p = new_enode(partition, 2*data.ns, dim);
 495                     for (j = 0; j < data.ns; ++j) {
 496                         EVALUE_SET_DOMAIN(tmp.x.p->arr[2*j], data.s[j].D);
 497                         value_clear(tmp.x.p->arr[2*j+1].d);
 498                         tmp.x.p->arr[2*j+1] = *data.s[j].E;
 499                         free(data.s[j].E);
 500                     }
 501                 }
 502             }
 503
 504             if (nvar > 1) {
 505                 evalue *res = Bernoulli_sum_evalue(&tmp, nvar-1, options);
 506                 eadd(res, sum);
 507                 free_evalue_refs(res);
 508                 free(res);
 509             } else
 510                 eadd(&tmp, sum);
 511
 512             free_evalue_refs(&tmp);
 513             D->next = next;;
 514         }
 515     }
 516
 517     free(data.s);
 518
 519     reduce_evalue(sum);
 520     return sum;
 521 }
 522
 523 evalue *Bernoulli_sum(Polyhedron *P, Polyhedron *C,
 524                         struct barvinok_options *options)
 525 {
 526     evalue e;
 527     evalue *sum;
 528
 529     value_init(e.d);
 530     e.x.p = new_enode(partition, 2, P->Dimension);
 531     EVALUE_SET_DOMAIN(e.x.p->arr[0], Polyhedron_Copy(P));
 532     evalue_set_si(&e.x.p->arr[1], 1, 1);
 533     sum = Bernoulli_sum_evalue(&e, P->Dimension - C->Dimension, options);
 534     free_evalue_refs(&e);
 535     return sum;
 536 }