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 bernoulli;
  10
  11 struct bernoulli *bernoulli_compute(int n)
  12 {
  13     int i, j;
  14     Value lcm, factor, tmp;
  15
  16     if (n < bernoulli.n)
  17         return &bernoulli;
  18
  19     if (n >= bernoulli.size) {
  20         int size = n + 5;
  21         Vector *b;
  22         Vector **sum;
  23
  24         b = Vector_Alloc(size);
  25         if (bernoulli.n) {
  26             Vector_Copy(bernoulli.b_num->p, b->p, bernoulli.n);
  27             Vector_Free(bernoulli.b_num);
  28         }
  29         bernoulli.b_num = b;
  30         b = Vector_Alloc(size);
  31         if (bernoulli.n) {
  32             Vector_Copy(bernoulli.b_den->p, b->p, bernoulli.n);
  33             Vector_Free(bernoulli.b_den);
  34         }
  35         bernoulli.b_den = b;
  36         sum = ALLOCN(Vector *, size);
  37         for (i = 0; i < bernoulli.n; ++i)
  38             sum[i] = bernoulli.sum[i];
  39         free(bernoulli.sum);
  40         bernoulli.sum = sum;
  41
  42         bernoulli.size = size;
  43     }
  44     value_init(lcm);
  45     value_init(factor);
  46     value_init(tmp);
  47     value_set_si(lcm, 1);
  48     for (i = 0; i < bernoulli.n; ++i)
  49         value_lcm(lcm, bernoulli.b_den->p[i], &lcm);
  50     for (i = bernoulli.n; i <= n; ++i) {
  51         if (i == 0) {
  52             value_set_si(bernoulli.b_num->p[0], 1);
  53             value_set_si(bernoulli.b_den->p[0], 1);
  54             continue;
  55         }
  56         value_set_si(bernoulli.b_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, lcm, bernoulli.b_den->p[j]);
  62             value_multiply(tmp, tmp, bernoulli.b_num->p[j]);
  63             value_multiply(tmp, tmp, factor);
  64             value_addto(bernoulli.b_num->p[i], bernoulli.b_num->p[i], tmp);
  65         }
  66         mpz_mul_ui(bernoulli.b_den->p[i], lcm, i+1);
  67         Gcd(bernoulli.b_num->p[i], bernoulli.b_den->p[i], &tmp);
  68         if (value_notone_p(tmp)) {
  69             value_division(bernoulli.b_num->p[i], bernoulli.b_num->p[i], tmp);
  70             value_division(bernoulli.b_den->p[i], bernoulli.b_den->p[i], tmp);
  71         }
  72         value_lcm(lcm, bernoulli.b_den->p[i], &lcm);
  73     }
  74     for (i = bernoulli.n; i <= n; ++i) {
  75         bernoulli.sum[i] = Vector_Alloc(i+3);
  76         value_assign(bernoulli.sum[i]->p[i+1], lcm);
  77         mpz_mul_ui(bernoulli.sum[i]->p[i+2], lcm, i+1);
  78         value_set_si(factor, 1);
  79         for (j = 1; j <= i; ++j) {
  80             mpz_mul_ui(factor, factor, i+2-j);
  81             mpz_divexact_ui(factor, factor, j);
  82             value_division(bernoulli.sum[i]->p[i+1-j], lcm, bernoulli.b_den->p[j]);
  83             value_multiply(bernoulli.sum[i]->p[i+1-j],
  84                            bernoulli.sum[i]->p[i+1-j], bernoulli.b_num->p[j]);
  85             value_multiply(bernoulli.sum[i]->p[i+1-j],
  86                            bernoulli.sum[i]->p[i+1-j], factor);
  87         }
  88         Vector_Normalize(bernoulli.sum[i]->p, i+3);
  89     }
  90     bernoulli.n = n+1;
  91     value_clear(lcm);
  92     value_clear(factor);
  93     value_clear(tmp);
  94
  95     return &bernoulli;
  96 }
  97
  98 /* shift variables in polynomial one down */
  99 static void shift(evalue *e)
 100 {
 101     int i;
 102     if (value_notzero_p(e->d))
 103         return;
 104     assert(e->x.p->type == polynomial);
 105     assert(e->x.p->pos > 1);
 106     e->x.p->pos--;
 107     for (i = 0; i < e->x.p->size; ++i)
 108         shift(&e->x.p->arr[i]);
 109 }
 110
 111 static evalue *shifted_copy(evalue *src)
 112 {
 113     evalue *e = ALLOC(evalue);
 114     value_init(e->d);
 115     evalue_copy(e, src);
 116     shift(e);
 117     return e;
 118 }
 119
 120 static evalue *power_sums(struct bernoulli *bernoulli, evalue *poly,
 121                           Vector *arg, Value denom, int neg, int alt_neg)
 122 {
 123     int i;
 124     evalue *base = affine2evalue(arg->p, denom, arg->Size-1);
 125     evalue *sum = evalue_zero();
 126
 127     for (i = 1; i < poly->x.p->size; ++i) {
 128         evalue *term = evalue_polynomial(bernoulli->sum[i], base);
 129         evalue *factor = shifted_copy(&poly->x.p->arr[i]);
 130         emul(factor, term);
 131         if (alt_neg && (i % 2))
 132             evalue_negate(term);
 133         eadd(term, sum);
 134         free_evalue_refs(factor);
 135         free_evalue_refs(term);
 136         free(factor);
 137         free(term);
 138     }
 139     if (neg)
 140         evalue_negate(sum);
 141     free_evalue_refs(base);
 142     free(base);
 143
 144     return sum;
 145 }
 146
 147 struct section { Polyhedron *D; evalue *E; };
 148
 149 struct Bernoulli_data {
 150     unsigned MaxRays;
 151     struct section *s;
 152     int size;
 153     int ns;
 154     evalue *e;
 155 };
 156
 157 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
 158 {
 159     struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
 160     Matrix *M2;
 161     Polyhedron *T;
 162     evalue *factor = NULL;
 163     evalue *linear = NULL;
 164     int constant = 0;
 165     Value tmp;
 166     unsigned dim = M->NbColumns-2;
 167     Vector *row;
 168
 169     assert(data->ns < data->size);
 170
 171     M2 = Matrix_Copy(M);
 172     T = Constraints2Polyhedron(M2, data->MaxRays);
 173     Matrix_Free(M2);
 174
 175     POL_ENSURE_VERTICES(T);
 176     if (emptyQ(T)) {
 177         Polyhedron_Free(T);
 178         return;
 179     }
 180
 181     assert(lower != upper);
 182
 183     row = Vector_Alloc(dim+1);
 184     value_init(tmp);
 185     if (value_notzero_p(data->e->d)) {
 186         factor = data->e;
 187         constant = 1;
 188     } else {
 189         assert(data->e->x.p->type == polynomial);
 190         if (data->e->x.p->pos > 1) {
 191             factor = shifted_copy(data->e);
 192             constant = 1;
 193         } else
 194             factor = shifted_copy(&data->e->x.p->arr[0]);
 195     }
 196     if (!EVALUE_IS_ZERO(*factor)) {
 197         value_absolute(tmp, upper[1]);
 198         /* upper - lower */
 199         Vector_Combine(lower+2, upper+2, row->p, tmp, lower[1], dim+1);
 200         value_multiply(tmp, tmp, lower[1]);
 201         /* upper - lower + 1 */
 202         value_addto(row->p[dim], row->p[dim], tmp);
 203
 204         linear = affine2evalue(row->p, tmp, dim);
 205         emul(factor, linear);
 206     } else
 207         linear = evalue_zero();
 208
 209     if (constant) {
 210         data->s[data->ns].E = linear;
 211         data->s[data->ns].D = T;
 212         ++data->ns;
 213     } else {
 214         evalue *poly_u = NULL, *poly_l = NULL;
 215         Polyhedron *D;
 216         struct bernoulli *bernoulli;
 217         assert(data->e->x.p->type == polynomial);
 218         assert(data->e->x.p->pos == 1);
 219         bernoulli = bernoulli_compute(data->e->x.p->size-1);
 220         /* lower is the constraint
 221          *                       b i - l' >= 0          i >= l'/b = l
 222          * upper is the constraint
 223          *                      -a i + u' >= 0          i <= u'/a = u
 224          */
 225         M2 = Matrix_Alloc(3, 2+T->Dimension);
 226         value_set_si(M2->p[0][0], 1);
 227         value_set_si(M2->p[1][0], 1);
 228         value_set_si(M2->p[2][0], 1);
 229         /* Case 1:
 230          *              1 <= l          l' - b >= 0
 231          */
 232         Vector_Oppose(lower+2, M2->p[0]+1, T->Dimension+1);
 233         value_subtract(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension],
 234                       lower[1]);
 235         D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
 236         if (emptyQ2(D))
 237             Polyhedron_Free(D);
 238         else {
 239             evalue *extra;
 240             if (!poly_u) {
 241                 Vector_Copy(upper+2, row->p, dim+1);
 242                 value_oppose(tmp, upper[1]);
 243                 value_addto(row->p[dim], row->p[dim], tmp);
 244                 poly_u = power_sums(bernoulli, data->e, row, tmp, 0, 0);
 245             }
 246             Vector_Oppose(lower+2, row->p, dim+1);
 247             extra = power_sums(bernoulli, data->e, row, lower[1], 1, 0);
 248             eadd(poly_u, extra);
 249             eadd(linear, extra);
 250
 251             data->s[data->ns].E = extra;
 252             data->s[data->ns].D = D;
 253             ++data->ns;
 254         }
 255
 256         /* Case 2:
 257          *              1 <= -u         -u' - a >= 0
 258          */
 259         Vector_Oppose(upper+2, M2->p[0]+1, T->Dimension+1);
 260         value_addto(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension],
 261                       upper[1]);
 262         D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
 263         if (emptyQ2(D))
 264             Polyhedron_Free(D);
 265         else {
 266             evalue *extra;
 267             if (!poly_l) {
 268                 Vector_Copy(lower+2, row->p, dim+1);
 269                 value_addto(row->p[dim], row->p[dim], lower[1]);
 270                 poly_l = power_sums(bernoulli, data->e, row, lower[1], 0, 1);
 271             }
 272             Vector_Oppose(upper+2, row->p, dim+1);
 273             value_oppose(tmp, upper[1]);
 274             extra = power_sums(bernoulli, data->e, row, tmp, 1, 1);
 275             eadd(poly_l, extra);
 276             eadd(linear, extra);
 277
 278             data->s[data->ns].E = extra;
 279             data->s[data->ns].D = D;
 280             ++data->ns;
 281         }
 282
 283         /* Case 3:
 284          *              u >= 0          u' >= 0
 285          *              -l >= 0         -l' >= 0
 286          */
 287         Vector_Copy(upper+2, M2->p[0]+1, T->Dimension+1);
 288         Vector_Copy(lower+2, M2->p[1]+1, T->Dimension+1);
 289         D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
 290         if (emptyQ2(D))
 291             Polyhedron_Free(D);
 292         else {
 293             if (!poly_l) {
 294                 Vector_Copy(lower+2, row->p, dim+1);
 295                 value_addto(row->p[dim], row->p[dim], lower[1]);
 296                 poly_l = power_sums(bernoulli, data->e, row, lower[1], 0, 1);
 297             }
 298             if (!poly_u) {
 299                 Vector_Copy(upper+2, row->p, dim+1);
 300                 value_oppose(tmp, upper[1]);
 301                 value_addto(row->p[dim], row->p[dim], tmp);
 302                 poly_u = power_sums(bernoulli, data->e, row, tmp, 0, 0);
 303             }
 304
 305             data->s[data->ns].E = ALLOC(evalue);
 306             value_init(data->s[data->ns].E->d);
 307             evalue_copy(data->s[data->ns].E, poly_u);
 308             eadd(poly_l, data->s[data->ns].E);
 309             eadd(linear, data->s[data->ns].E);
 310             data->s[data->ns].D = D;
 311             ++data->ns;
 312         }
 313
 314         /* Case 4:
 315          *              l < 1           -l' + b - 1 >= 0
 316          *              0 < l           l' - 1 >= 0
 317          */
 318         Vector_Copy(lower+2, M2->p[0]+1, T->Dimension+1);
 319         value_addto(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension], lower[1]);
 320         value_decrement(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension]);
 321         Vector_Oppose(lower+2, M2->p[1]+1, T->Dimension+1);
 322         value_decrement(M2->p[1][1+T->Dimension], M2->p[1][1+T->Dimension]);
 323         D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
 324         if (emptyQ2(D))
 325             Polyhedron_Free(D);
 326         else {
 327             if (!poly_u) {
 328                 Vector_Copy(upper+2, row->p, dim+1);
 329                 value_oppose(tmp, upper[1]);
 330                 value_addto(row->p[dim], row->p[dim], tmp);
 331                 poly_u = power_sums(bernoulli, data->e, row, tmp, 0, 0);
 332             }
 333
 334             eadd(linear, poly_u);
 335             data->s[data->ns].E = poly_u;
 336             poly_u = NULL;
 337             data->s[data->ns].D = D;
 338             ++data->ns;
 339         }
 340
 341         /* Case 5:
 342          *              -u < 1          u' + a - 1 >= 0
 343          *              0 < -u          -u' - 1 >= 0
 344          *              l <= 0          -l' >= 0
 345          */
 346         Vector_Copy(upper+2, M2->p[0]+1, T->Dimension+1);
 347         value_subtract(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension],
 348                         upper[1]);
 349         value_decrement(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension]);
 350         Vector_Oppose(upper+2, M2->p[1]+1, T->Dimension+1);
 351         value_decrement(M2->p[1][1+T->Dimension], M2->p[1][1+T->Dimension]);
 352         Vector_Copy(lower+2, M2->p[2]+1, T->Dimension+1);
 353         D = AddConstraints(M2->p_Init, 3, T, data->MaxRays);
 354         if (emptyQ2(D))
 355             Polyhedron_Free(D);
 356         else {
 357             if (!poly_l) {
 358                 Vector_Copy(lower+2, row->p, dim+1);
 359                 value_addto(row->p[dim], row->p[dim], lower[1]);
 360                 poly_l = power_sums(bernoulli, data->e, row, lower[1], 0, 1);
 361             }
 362
 363             eadd(linear, poly_l);
 364             data->s[data->ns].E = poly_l;
 365             poly_l = NULL;
 366             data->s[data->ns].D = D;
 367             ++data->ns;
 368         }
 369
 370         Matrix_Free(M2);
 371         Polyhedron_Free(T);
 372         if (poly_l) {
 373             free_evalue_refs(poly_l);
 374             free(poly_l);
 375         }
 376         if (poly_u) {
 377             free_evalue_refs(poly_u);
 378             free(poly_u);
 379         }
 380         free_evalue_refs(linear);
 381         free(linear);
 382     }
 383     if (factor != data->e) {
 384         free_evalue_refs(factor);
 385         free(factor);
 386     }
 387     value_clear(tmp);
 388     Vector_Free(row);
 389 }
 390
 391 evalue *Bernoulli_sum_evalue(evalue *e, unsigned nvar,
 392                              struct barvinok_options *options)
 393 {
 394     struct Bernoulli_data data;
 395     int i, j;
 396     evalue *sum = evalue_zero();
 397
 398     if (EVALUE_IS_ZERO(*e))
 399         return sum;
 400
 401     if (nvar == 0) {
 402         eadd(e, sum);
 403         return sum;
 404     }
 405
 406     assert(value_zero_p(e->d));
 407     assert(e->x.p->type == partition);
 408
 409     data.size = e->x.p->size * 2 + 16;
 410     data.s = ALLOCN(struct section, data.size);
 411     data.MaxRays = options->MaxRays;
 412
 413     for (i = 0; i < e->x.p->size/2; ++i) {
 414         Polyhedron *D;
 415         for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
 416             unsigned dim = D->Dimension - 1;
 417             Polyhedron *next = D->next;
 418             evalue tmp;
 419             D->next = NULL;
 420
 421             value_init(tmp.d);
 422             value_set_si(tmp.d, 0);
 423
 424             if (value_zero_p(D->Constraint[0][0]) &&
 425                     value_notzero_p(D->Constraint[0][1])) {
 426                 tmp.x.p = new_enode(partition, 2, dim);
 427                 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Project(D, dim));
 428                 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
 429                 reduce_evalue_in_domain(&tmp.x.p->arr[1], D);
 430                 shift(&tmp.x.p->arr[1]);
 431             } else {
 432                 data.ns = 0;
 433                 data.e = &e->x.p->arr[2*i+1];
 434
 435                 for_each_lower_upper_bound(D, Bernoulli_cb, &data);
 436
 437                 if (data.ns == 0)
 438                     evalue_set_si(&tmp, 0, 1);
 439                 else {
 440                     tmp.x.p = new_enode(partition, 2*data.ns, dim);
 441                     for (j = 0; j < data.ns; ++j) {
 442                         EVALUE_SET_DOMAIN(tmp.x.p->arr[2*j], data.s[j].D);
 443                         value_clear(tmp.x.p->arr[2*j+1].d);
 444                         tmp.x.p->arr[2*j+1] = *data.s[j].E;
 445                         free(data.s[j].E);
 446                     }
 447                 }
 448             }
 449
 450             if (nvar > 1) {
 451                 evalue *res = Bernoulli_sum_evalue(&tmp, nvar-1, options);
 452                 eadd(res, sum);
 453                 free_evalue_refs(res);
 454                 free(res);
 455             } else
 456                 eadd(&tmp, sum);
 457
 458             free_evalue_refs(&tmp);
 459             D->next = next;;
 460         }
 461     }
 462
 463     free(data.s);
 464
 465     reduce_evalue(sum);
 466     return sum;
 467 }
 468
 469 evalue *Bernoulli_sum(Polyhedron *P, Polyhedron *C,
 470                         struct barvinok_options *options)
 471 {
 472     evalue e;
 473     evalue *sum;
 474
 475     value_init(e.d);
 476     e.x.p = new_enode(partition, 2, P->Dimension);
 477     EVALUE_SET_DOMAIN(e.x.p->arr[0], Polyhedron_Copy(P));
 478     evalue_set_si(&e.x.p->arr[1], 1, 1);
 479     sum = Bernoulli_sum_evalue(&e, P->Dimension - C->Dimension, options);
 480     free_evalue_refs(&e);
 481     return sum;
 482 }