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