bernoulli.c

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