?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Bernoulli_sum_evalue: make sure no empty partitions are created
[barvinok.git] / bernoulli.c
blob1293a2f5fd19126207aad1c38268c0f5cc34ea80
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 /* Given a constraint (cst_affine) a x + b y + c >= 0, compate a constraint (c)
200 * +- (b y + c) +- a >=,> 0
201 * ^ ^ ^
202 * | | strict
203 * sign_affine sign_cst
204 */
205 static void bound_constraint(Value *c, unsigned dim,
206 Value *cst_affine,
207 int sign_affine, int sign_cst, int strict)
208 {
209 if (sign_affine == -1)
210 Vector_Oppose(cst_affine+1, c, dim+1);
211 else
212 Vector_Copy(cst_affine+1, c, dim+1);
213
214 if (sign_cst == -1)
215 value_subtract(c[dim], c[dim], cst_affine[0]);
216 else if (sign_cst == 1)
217 value_addto(c[dim], c[dim], cst_affine[0]);
218
219 if (strict)
220 value_decrement(c[dim], c[dim]);
221 }
222
223 struct Bernoulli_data {
224 unsigned MaxRays;
225 struct evalue_section *s;
226 int size;
227 int ns;
228 evalue *e;
229 };
230
231 static evalue *compute_poly_u(evalue *poly_u, Value *upper, Vector *row,
232 unsigned dim, Value tmp,
233 struct poly_list *faulhaber,
234 struct Bernoulli_data *data)
235 {
236 if (poly_u)
237 return poly_u;
238 Vector_Copy(upper+2, row->p, dim+1);
239 value_oppose(tmp, upper[1]);
240 value_addto(row->p[dim], row->p[dim], tmp);
241 return power_sums(faulhaber, data->e, row, tmp, 0, 0);
242 }
243
244 static evalue *compute_poly_l(evalue *poly_l, Value *lower, Vector *row,
245 unsigned dim,
246 struct poly_list *faulhaber,
247 struct Bernoulli_data *data)
248 {
249 if (poly_l)
250 return poly_l;
251 Vector_Copy(lower+2, row->p, dim+1);
252 value_addto(row->p[dim], row->p[dim], lower[1]);
253 return power_sums(faulhaber, data->e, row, lower[1], 0, 1);
254 }
255
256 static void Bernoulli_init(unsigned n, void *cb_data)
257 {
258 struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
259 int cases = 5;
260
261 if (cases * n <= data->size)
262 return;
263
264 data->size = cases * (n + 16);
265 data->s = REALLOCN(data->s, struct evalue_section, data->size);
266 }
267
268 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
269 {
270 struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
271 Matrix *M2;
272 Polyhedron *T;
273 evalue *factor = NULL;
274 evalue *linear = NULL;
275 int constant = 0;
276 Value tmp;
277 unsigned dim = M->NbColumns-2;
278 Vector *row;
279 int cases = 5;
280
281 assert(lower);
282 assert(upper);
283 assert(data->ns + cases <= data->size);
284
285 M2 = Matrix_Copy(M);
286 T = Constraints2Polyhedron(M2, data->MaxRays);
287 Matrix_Free(M2);
288
289 POL_ENSURE_VERTICES(T);
290 if (emptyQ(T)) {
291 Polyhedron_Free(T);
292 return;
293 }
294
295 assert(lower != upper);
296
297 row = Vector_Alloc(dim+1);
298 value_init(tmp);
299 if (value_notzero_p(data->e->d)) {
300 factor = data->e;
301 constant = 1;
302 } else {
303 assert(data->e->x.p->type == polynomial);
304 if (data->e->x.p->pos > 1) {
305 factor = shifted_copy(data->e);
306 constant = 1;
307 } else
308 factor = shifted_copy(&data->e->x.p->arr[0]);
309 }
310 if (!EVALUE_IS_ZERO(*factor)) {
311 value_absolute(tmp, upper[1]);
312 /* upper - lower */
313 Vector_Combine(lower+2, upper+2, row->p, tmp, lower[1], dim+1);
314 value_multiply(tmp, tmp, lower[1]);
315 /* upper - lower + 1 */
316 value_addto(row->p[dim], row->p[dim], tmp);
317
318 linear = affine2evalue(row->p, tmp, dim);
319 emul(factor, linear);
320 } else
321 linear = evalue_zero();
322
323 if (constant) {
324 data->s[data->ns].E = linear;
325 data->s[data->ns].D = T;
326 ++data->ns;
327 } else {
328 evalue *poly_u = NULL, *poly_l = NULL;
329 Polyhedron *D;
330 struct poly_list *faulhaber;
331 assert(data->e->x.p->type == polynomial);
332 assert(data->e->x.p->pos == 1);
333 faulhaber = faulhaber_compute(data->e->x.p->size-1);
334 /* lower is the constraint
335 * b i - l' >= 0 i >= l'/b = l
336 * upper is the constraint
337 * -a i + u' >= 0 i <= u'/a = u
338 */
339 M2 = Matrix_Alloc(2, 2+T->Dimension);
340 value_set_si(M2->p[0][0], 1);
341 value_set_si(M2->p[1][0], 1);
342 /* Case 1:
343 * 1 <= l l' - b >= 0
344 */
345 bound_constraint(M2->p[0]+1, T->Dimension, lower+1, -1, -1, 0);
346 D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
347 POL_ENSURE_VERTICES(D);
348 if (emptyQ2(D))
349 Polyhedron_Free(D);
350 else {
351 evalue *extra;
352 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
353 faulhaber, data);
354 Vector_Oppose(lower+2, row->p, dim+1);
355 extra = power_sums(faulhaber, data->e, row, lower[1], 1, 0);
356 eadd(poly_u, extra);
357 eadd(linear, extra);
358
359 data->s[data->ns].E = extra;
360 data->s[data->ns].D = D;
361 ++data->ns;
362 }
363
364 /* Case 2:
365 * 1 <= -u -u' - a >= 0
366 */
367 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, -1, 1, 0);
368 D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
369 POL_ENSURE_VERTICES(D);
370 if (emptyQ2(D))
371 Polyhedron_Free(D);
372 else {
373 evalue *extra;
374 poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
375 Vector_Oppose(upper+2, row->p, dim+1);
376 value_oppose(tmp, upper[1]);
377 extra = power_sums(faulhaber, data->e, row, tmp, 1, 1);
378 eadd(poly_l, extra);
379 eadd(linear, extra);
380
381 data->s[data->ns].E = extra;
382 data->s[data->ns].D = D;
383 ++data->ns;
384 }
385
386 /* Case 3:
387 * u >= 0 u' >= 0
388 * -l >= 0 -l' >= 0
389 */
390 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, 0, 0);
391 bound_constraint(M2->p[1]+1, T->Dimension, lower+1, 1, 0, 0);
392 D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
393 POL_ENSURE_VERTICES(D);
394 if (emptyQ2(D))
395 Polyhedron_Free(D);
396 else {
397 poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
398 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
399 faulhaber, data);
400
401 data->s[data->ns].E = ALLOC(evalue);
402 value_init(data->s[data->ns].E->d);
403 evalue_copy(data->s[data->ns].E, poly_u);
404 eadd(poly_l, data->s[data->ns].E);
405 eadd(linear, data->s[data->ns].E);
406 data->s[data->ns].D = D;
407 ++data->ns;
408 }
409
410 /* Case 4:
411 * l < 1 -l' + b - 1 >= 0
412 * 0 < l l' - 1 >= 0
413 */
414 bound_constraint(M2->p[0]+1, T->Dimension, lower+1, 1, 1, 1);
415 bound_constraint(M2->p[1]+1, T->Dimension, lower+1, -1, 0, 1);
416 D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
417 POL_ENSURE_VERTICES(D);
418 if (emptyQ2(D))
419 Polyhedron_Free(D);
420 else {
421 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
422 faulhaber, data);
423
424 eadd(linear, poly_u);
425 data->s[data->ns].E = poly_u;
426 poly_u = NULL;
427 data->s[data->ns].D = D;
428 ++data->ns;
429 }
430
431 /* Case 5:
432 * -u < 1 u' + a - 1 >= 0
433 * 0 < -u -u' - 1 >= 0
434 */
435 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, -1, 1);
436 bound_constraint(M2->p[1]+1, T->Dimension, upper+1, -1, 0, 1);
437 D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
438 POL_ENSURE_VERTICES(D);
439 if (emptyQ2(D))
440 Polyhedron_Free(D);
441 else {
442 poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
443
444 eadd(linear, poly_l);
445 data->s[data->ns].E = poly_l;
446 poly_l = NULL;
447 data->s[data->ns].D = D;
448 ++data->ns;
449 }
450
451 Matrix_Free(M2);
452 Polyhedron_Free(T);
453 if (poly_l)
454 evalue_free(poly_l);
455 if (poly_u)
456 evalue_free(poly_u);
457 evalue_free(linear);
458 }
459 if (factor != data->e)
460 evalue_free(factor);
461 value_clear(tmp);
462 Vector_Free(row);
463 }
464
465 /* Looks for variable with integer bounds, i.e., with coefficients 0, 1 or -1.
466 * Returns 1 if such a variable is found and puts it in the first position,
467 * possibly changing *P_p and *E_p.
468 */
469 static int find_integer_bounds(Polyhedron **P_p, evalue **E_p, unsigned nvar)
470 {
471 Polyhedron *P = *P_p;
472 evalue *E = *E_p;
473 unsigned dim = P->Dimension;
474 int i, j;
475
476 for (i = 0; i < nvar; ++i) {
477 for (j = 0; j < P->NbConstraints; ++j) {
478 if (value_zero_p(P->Constraint[j][1+i]))
479 continue;
480 if (value_one_p(P->Constraint[j][1+i]))
481 continue;
482 if (value_mone_p(P->Constraint[j][1+i]))
483 continue;
484 break;
485 }
486 if (j == P->NbConstraints)
487 break;
488 }
489 if (i == nvar)
490 return 0;
491 if (i == 0)
492 return 1;
493 P = Polyhedron_Copy(P);
494 Polyhedron_ExchangeColumns(P, 1, 1+i);
495 *P_p = P;
496
497 if (value_zero_p(E->d)) {
498 evalue **subs;
499 subs = ALLOCN(evalue *, dim);
500 for (j = 0; j < dim; ++j)
501 subs[j] = evalue_var(j);
502 E = subs[0];
503 subs[0] = subs[i];
504 subs[i] = E;
505 E = evalue_dup(*E_p);
506 evalue_substitute(E, subs);
507 for (j = 0; j < dim; ++j)
508 evalue_free(subs[j]);
509 free(subs);
510 *E_p = E;
511 }
512
513 return 1;
514 }
515
516 static evalue *sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
517 struct Bernoulli_data *data,
518 struct barvinok_options *options)
519 {
520 unsigned dim = P->Dimension - 1;
521 evalue *res;
522
523 if (value_zero_p(P->Constraint[0][0]) &&
524 value_notzero_p(P->Constraint[0][1])) {
525 res = ALLOC(evalue);
526 value_init(res->d);
527 value_set_si(res->d, 0);
528 res->x.p = new_enode(partition, 2, dim);
529 EVALUE_SET_DOMAIN(res->x.p->arr[0], Polyhedron_Project(P, dim));
530 evalue_copy(&res->x.p->arr[1], E);
531 reduce_evalue_in_domain(&res->x.p->arr[1], P);
532 shift(&res->x.p->arr[1]);
533 } else {
534 data->ns = 0;
535 data->e = E;
536
537 for_each_lower_upper_bound(P, Bernoulli_init, Bernoulli_cb, data);
538
539 res = evalue_from_section_array(data->s, data->ns);
540 }
541
542 if (nvar > 1) {
543 evalue *tmp = Bernoulli_sum_evalue(res, nvar-1, options);
544 evalue_free(res);
545 res = tmp;
546 }
547
548 return res;
549 }
550
551 evalue *Bernoulli_sum_evalue(evalue *e, unsigned nvar,
552 struct barvinok_options *options)
553 {
554 struct Bernoulli_data data;
555 int i, j;
556 evalue *sum = evalue_zero();
557
558 if (EVALUE_IS_ZERO(*e))
559 return sum;
560
561 if (nvar == 0) {
562 eadd(e, sum);
563 return sum;
564 }
565
566 assert(value_zero_p(e->d));
567 assert(e->x.p->type == partition);
568
569 data.size = 16;
570 data.s = ALLOCN(struct evalue_section, data.size);
571 data.MaxRays = options->MaxRays;
572
573 for (i = 0; i < e->x.p->size/2; ++i) {
574 Polyhedron *D;
575 for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
576 evalue *E = &e->x.p->arr[2*i+1];
577 Polyhedron *P = D;
578 Polyhedron *next = D->next;
579 evalue *tmp;
580 int integer_bounds;
581
582 P->next = NULL;
583
584 integer_bounds = find_integer_bounds(&P, &E, nvar);
585 if (options->approximation_method == BV_APPROX_NONE &&
586 !integer_bounds) {
587 evalue_free(sum);
588 sum = NULL;
589 } else {
590 evalue *tmp = sum_over_polytope(P, E, nvar, &data, options);
591 eadd(tmp, sum);
592 evalue_free(tmp);
593 }
594
595 if (P != D)
596 Polyhedron_Free(P);
597 if (E != &e->x.p->arr[2*i+1])
598 evalue_free(E);
599
600 D->next = next;;
601
602 if (!sum)
603 break;
604 }
605
606 if (!sum)
607 break;
608 }
609
610 free(data.s);
611
612 if (sum)
613 reduce_evalue(sum);
614 return sum;
615 }
616
617 evalue *Bernoulli_sum(Polyhedron *P, Polyhedron *C,
618 struct barvinok_options *options)
619 {
620 Polyhedron *CA, *D;
621 evalue e;
622 evalue *sum;
623
624 if (emptyQ(P) || emptyQ(C))
625 return evalue_zero();
626
627 CA = align_context(C, P->Dimension, options->MaxRays);
628 D = DomainIntersection(P, CA, options->MaxRays);
629 Domain_Free(CA);
630
631 if (emptyQ(D)) {
632 Domain_Free(D);
633 return evalue_zero();
634 }
635
636 value_init(e.d);
637 e.x.p = new_enode(partition, 2, P->Dimension);
638 EVALUE_SET_DOMAIN(e.x.p->arr[0], D);
639 evalue_set_si(&e.x.p->arr[1], 1, 1);
640 sum = Bernoulli_sum_evalue(&e, P->Dimension - C->Dimension, options);
641 free_evalue_refs(&e);
642 return sum;
643 }
lattice point enumeration library