remove_parameter_equalities: stop "removing" equalities if polyhedron is empty
[barvinok.git] / bernoulli.c
blob677a4cda068b80d80eaa130add5dda2aa3b43ebb
1 #include <assert.h>
2 #include <barvinok/options.h>
3 #include <barvinok/util.h>
4 #include "bernoulli.h"
6 #define ALLOC(type) (type*)malloc(sizeof(type))
7 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
9 static struct bernoulli_coef bernoulli_coef;
10 static struct poly_list bernoulli;
11 static struct poly_list faulhaber;
13 struct bernoulli_coef *bernoulli_coef_compute(int n)
15 int i, j;
16 Value factor, tmp;
18 if (n < bernoulli_coef.n)
19 return &bernoulli_coef;
21 if (n >= bernoulli_coef.size) {
22 int size = 3*(n + 5)/2;
23 Vector *b;
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);
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);
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);
42 bernoulli_coef.lcm = b;
44 bernoulli_coef.size = size;
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;
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);
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);
74 value_lcm(bernoulli_coef.lcm->p[i],
75 bernoulli_coef.lcm->p[i-1], bernoulli_coef.den->p[i]);
77 bernoulli_coef.n = n+1;
78 value_clear(factor);
79 value_clear(tmp);
81 return &bernoulli_coef;
85 * Compute either Bernoulli B_n or Faulhaber F_n polynomials.
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}
90 static struct poly_list *bernoulli_faulhaber_compute(int n, struct poly_list *pl,
91 int faulhaber)
93 int i, j;
94 Value factor;
95 struct bernoulli_coef *bc;
97 if (n < pl->n)
98 return pl;
100 if (n >= pl->size) {
101 int size = 3*(n + 5)/2;
102 Vector **poly;
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;
110 pl->size = size;
113 bc = bernoulli_coef_compute(n);
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);
134 Vector_Normalize(pl->poly[i]->p, i+faulhaber+2);
136 value_clear(factor);
137 pl->n = n+1;
139 return pl;
142 struct poly_list *bernoulli_compute(int n)
144 return bernoulli_faulhaber_compute(n, &bernoulli, 0);
147 struct poly_list *faulhaber_compute(int n)
149 return bernoulli_faulhaber_compute(n, &faulhaber, 1);
152 /* shift variables in polynomial one down */
153 static void shift(evalue *e)
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]);
165 static evalue *shifted_copy(evalue *src)
167 evalue *e = ALLOC(evalue);
168 value_init(e->d);
169 evalue_copy(e, src);
170 shift(e);
171 return e;
174 static evalue *power_sums(struct poly_list *faulhaber, evalue *poly,
175 Vector *arg, Value denom, int neg, int alt_neg)
177 int i;
178 evalue *base = affine2evalue(arg->p, denom, arg->Size-1);
179 evalue *sum = evalue_zero();
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);
191 if (neg)
192 evalue_negate(sum);
193 evalue_free(base);
195 return sum;
198 struct Bernoulli_data {
199 unsigned MaxRays;
200 struct evalue_section *s;
201 int size;
202 int ns;
203 evalue *e;
206 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
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;
218 assert(data->ns < data->size);
220 M2 = Matrix_Copy(M);
221 T = Constraints2Polyhedron(M2, data->MaxRays);
222 Matrix_Free(M2);
224 POL_ENSURE_VERTICES(T);
225 if (emptyQ(T)) {
226 Polyhedron_Free(T);
227 return;
230 assert(lower != upper);
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]);
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);
253 linear = affine2evalue(row->p, tmp, dim);
254 emul(factor, linear);
255 } else
256 linear = evalue_zero();
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
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
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);
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);
300 data->s[data->ns].E = extra;
301 data->s[data->ns].D = D;
302 ++data->ns;
305 /* Case 2:
306 * 1 <= -u -u' - a >= 0
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);
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);
327 data->s[data->ns].E = extra;
328 data->s[data->ns].D = D;
329 ++data->ns;
332 /* Case 3:
333 * u >= 0 u' >= 0
334 * -l >= 0 -l' >= 0
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);
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);
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;
363 /* Case 4:
364 * l < 1 -l' + b - 1 >= 0
365 * 0 < l l' - 1 >= 0
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);
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;
390 /* Case 5:
391 * -u < 1 u' + a - 1 >= 0
392 * 0 < -u -u' - 1 >= 0
393 * l <= 0 -l' >= 0
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);
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;
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);
427 if (factor != data->e)
428 evalue_free(factor);
429 value_clear(tmp);
430 Vector_Free(row);
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.
437 static int find_integer_bounds(Polyhedron **P_p, evalue **E_p, unsigned nvar)
439 Polyhedron *P = *P_p;
440 evalue *E = *E_p;
441 unsigned dim = P->Dimension;
442 int i, j;
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;
454 if (j == P->NbConstraints)
455 break;
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;
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;
481 return 1;
484 static evalue *sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
485 struct Bernoulli_data *data,
486 struct barvinok_options *options)
488 unsigned dim = P->Dimension - 1;
489 evalue *res;
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;
505 for_each_lower_upper_bound(P, Bernoulli_cb, data);
507 res = evalue_from_section_array(data->s, data->ns);
510 if (nvar > 1) {
511 evalue *tmp = Bernoulli_sum_evalue(res, nvar-1, options);
512 evalue_free(res);
513 res = tmp;
516 return res;
519 evalue *Bernoulli_sum_evalue(evalue *e, unsigned nvar,
520 struct barvinok_options *options)
522 struct Bernoulli_data data;
523 int i, j;
524 evalue *sum = evalue_zero();
526 if (EVALUE_IS_ZERO(*e))
527 return sum;
529 if (nvar == 0) {
530 eadd(e, sum);
531 return sum;
534 assert(value_zero_p(e->d));
535 assert(e->x.p->type == partition);
537 data.size = e->x.p->size * 2 + 16;
538 data.s = ALLOCN(struct evalue_section, data.size);
539 data.MaxRays = options->MaxRays;
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;
550 P->next = NULL;
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);
563 if (P != D)
564 Polyhedron_Free(P);
565 if (E != &e->x.p->arr[2*i+1])
566 evalue_free(E);
568 D->next = next;;
570 if (!sum)
571 break;
574 if (!sum)
575 break;
578 free(data.s);
580 if (sum)
581 reduce_evalue(sum);
582 return sum;
585 evalue *Bernoulli_sum(Polyhedron *P, Polyhedron *C,
586 struct barvinok_options *options)
588 evalue e;
589 evalue *sum;
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;