update AUTHORS
[barvinok.git] / bernoulli.c
blobb5cf1d0872f11099302c73b0a641aad02aa9b47e
1 #include <assert.h>
2 #include <barvinok/barvinok.h>
3 #include <barvinok/options.h>
4 #include <barvinok/util.h>
5 #include "bernoulli.h"
6 #include "lattice_point.h"
7 #include "section_array.h"
8 #include "summate.h"
10 #define ALLOC(type) (type*)malloc(sizeof(type))
11 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
13 static struct bernoulli_coef bernoulli_coef;
14 static struct poly_list bernoulli;
15 static struct poly_list faulhaber;
17 struct bernoulli_coef *bernoulli_coef_compute(int n)
19 int i, j;
20 Value factor, tmp;
22 if (n < bernoulli_coef.n)
23 return &bernoulli_coef;
25 if (n >= bernoulli_coef.size) {
26 int size = 3*(n + 5)/2;
27 Vector *b;
29 b = Vector_Alloc(size);
30 if (bernoulli_coef.n) {
31 Vector_Copy(bernoulli_coef.num->p, b->p, bernoulli_coef.n);
32 Vector_Free(bernoulli_coef.num);
34 bernoulli_coef.num = b;
35 b = Vector_Alloc(size);
36 if (bernoulli_coef.n) {
37 Vector_Copy(bernoulli_coef.den->p, b->p, bernoulli_coef.n);
38 Vector_Free(bernoulli_coef.den);
40 bernoulli_coef.den = b;
41 b = Vector_Alloc(size);
42 if (bernoulli_coef.n) {
43 Vector_Copy(bernoulli_coef.lcm->p, b->p, bernoulli_coef.n);
44 Vector_Free(bernoulli_coef.lcm);
46 bernoulli_coef.lcm = b;
48 bernoulli_coef.size = size;
50 value_init(factor);
51 value_init(tmp);
52 for (i = bernoulli_coef.n; i <= n; ++i) {
53 if (i == 0) {
54 value_set_si(bernoulli_coef.num->p[0], 1);
55 value_set_si(bernoulli_coef.den->p[0], 1);
56 value_set_si(bernoulli_coef.lcm->p[0], 1);
57 continue;
59 value_set_si(bernoulli_coef.num->p[i], 0);
60 value_set_si(factor, -(i+1));
61 for (j = i-1; j >= 0; --j) {
62 mpz_mul_ui(factor, factor, j+1);
63 mpz_divexact_ui(factor, factor, i+1-j);
64 value_division(tmp, bernoulli_coef.lcm->p[i-1],
65 bernoulli_coef.den->p[j]);
66 value_multiply(tmp, tmp, bernoulli_coef.num->p[j]);
67 value_multiply(tmp, tmp, factor);
68 value_addto(bernoulli_coef.num->p[i], bernoulli_coef.num->p[i], tmp);
70 mpz_mul_ui(bernoulli_coef.den->p[i], bernoulli_coef.lcm->p[i-1], i+1);
71 value_gcd(tmp, bernoulli_coef.num->p[i], bernoulli_coef.den->p[i]);
72 if (value_notone_p(tmp)) {
73 value_division(bernoulli_coef.num->p[i],
74 bernoulli_coef.num->p[i], tmp);
75 value_division(bernoulli_coef.den->p[i],
76 bernoulli_coef.den->p[i], tmp);
78 value_lcm(bernoulli_coef.lcm->p[i],
79 bernoulli_coef.lcm->p[i-1], bernoulli_coef.den->p[i]);
81 bernoulli_coef.n = n+1;
82 value_clear(factor);
83 value_clear(tmp);
85 return &bernoulli_coef;
89 * Compute either Bernoulli B_n or Faulhaber F_n polynomials.
91 * B_n = sum_{k=0}^n { n \choose k } b_k x^{n-k}
92 * F_n = 1/(n+1) sum_{k=0}^n { n+1 \choose k } b_k x^{n+1-k}
94 static struct poly_list *bernoulli_faulhaber_compute(int n, struct poly_list *pl,
95 int faulhaber)
97 int i, j;
98 Value factor;
99 struct bernoulli_coef *bc;
101 if (n < pl->n)
102 return pl;
104 if (n >= pl->size) {
105 int size = 3*(n + 5)/2;
106 Vector **poly;
108 poly = ALLOCN(Vector *, size);
109 for (i = 0; i < pl->n; ++i)
110 poly[i] = pl->poly[i];
111 free(pl->poly);
112 pl->poly = poly;
114 pl->size = size;
117 bc = bernoulli_coef_compute(n);
119 value_init(factor);
120 for (i = pl->n; i <= n; ++i) {
121 pl->poly[i] = Vector_Alloc(i+faulhaber+2);
122 value_assign(pl->poly[i]->p[i+faulhaber], bc->lcm->p[i]);
123 if (faulhaber)
124 mpz_mul_ui(pl->poly[i]->p[i+2], bc->lcm->p[i], i+1);
125 else
126 value_assign(pl->poly[i]->p[i+1], bc->lcm->p[i]);
127 value_set_si(factor, 1);
128 for (j = 1; j <= i; ++j) {
129 mpz_mul_ui(factor, factor, i+faulhaber+1-j);
130 mpz_divexact_ui(factor, factor, j);
131 value_division(pl->poly[i]->p[i+faulhaber-j],
132 bc->lcm->p[i], bc->den->p[j]);
133 value_multiply(pl->poly[i]->p[i+faulhaber-j],
134 pl->poly[i]->p[i+faulhaber-j], bc->num->p[j]);
135 value_multiply(pl->poly[i]->p[i+faulhaber-j],
136 pl->poly[i]->p[i+faulhaber-j], factor);
138 Vector_Normalize(pl->poly[i]->p, i+faulhaber+2);
140 value_clear(factor);
141 pl->n = n+1;
143 return pl;
146 struct poly_list *bernoulli_compute(int n)
148 return bernoulli_faulhaber_compute(n, &bernoulli, 0);
151 struct poly_list *faulhaber_compute(int n)
153 return bernoulli_faulhaber_compute(n, &faulhaber, 1);
156 static evalue *shifted_copy(const evalue *src)
158 evalue *e = ALLOC(evalue);
159 value_init(e->d);
160 evalue_copy(e, src);
161 evalue_shift_variables(e, 0, -1);
162 return e;
165 /* Computes the argument for the Faulhaber polynomial.
166 * If we are computing a polynomial approximation (exact == 0),
167 * then this is simply arg/denom.
168 * Otherwise, if neg == 0, then we are dealing with an upper bound
169 * and we need to compute floor(arg/denom) = arg/denom - { arg/denom }
170 * If neg == 1, then we are dealing with a lower bound
171 * and we need to compute ceil(arg/denom) = arg/denom + { -arg/denom }
173 * Modifies arg (if exact == 1).
175 static evalue *power_sums_base(Vector *arg, Value denom, int neg, int exact)
177 evalue *fract;
178 evalue *base = affine2evalue(arg->p, denom, arg->Size-1);
180 if (!exact || value_one_p(denom))
181 return base;
183 if (neg)
184 Vector_Oppose(arg->p, arg->p, arg->Size);
186 fract = fractional_part(arg->p, denom, arg->Size-1, NULL);
187 if (!neg) {
188 value_set_si(arg->p[0], -1);
189 evalue_mul(fract, arg->p[0]);
191 eadd(fract, base);
192 evalue_free(fract);
194 return base;
197 static evalue *power_sums(struct poly_list *faulhaber, const evalue *poly,
198 Vector *arg, Value denom, int neg, int alt_neg,
199 int exact)
201 int i;
202 evalue *base = power_sums_base(arg, denom, neg, exact);
203 evalue *sum = evalue_zero();
205 for (i = 1; i < poly->x.p->size; ++i) {
206 evalue *term = evalue_polynomial(faulhaber->poly[i], base);
207 evalue *factor = shifted_copy(&poly->x.p->arr[i]);
208 emul(factor, term);
209 if (alt_neg && (i % 2))
210 evalue_negate(term);
211 eadd(term, sum);
212 evalue_free(factor);
213 evalue_free(term);
215 if (neg)
216 evalue_negate(sum);
217 evalue_free(base);
219 return sum;
222 /* Given a constraint (cst_affine) a x + b y + c >= 0, compute a constraint (c)
223 * +- (b y + c) +- a >=,> 0
224 * ^ ^ ^
225 * | | strict
226 * sign_affine sign_cst
228 static void bound_constraint(Value *c, unsigned dim,
229 Value *cst_affine,
230 int sign_affine, int sign_cst, int strict)
232 if (sign_affine == -1)
233 Vector_Oppose(cst_affine+1, c, dim+1);
234 else
235 Vector_Copy(cst_affine+1, c, dim+1);
237 if (sign_cst == -1)
238 value_subtract(c[dim], c[dim], cst_affine[0]);
239 else if (sign_cst == 1)
240 value_addto(c[dim], c[dim], cst_affine[0]);
242 if (strict)
243 value_decrement(c[dim], c[dim]);
246 struct Bernoulli_data {
247 struct barvinok_options *options;
248 struct evalue_section_array *sections;
249 const evalue *e;
252 static evalue *compute_poly_u(evalue *poly_u, Value *upper, Vector *row,
253 unsigned dim, Value tmp,
254 struct poly_list *faulhaber,
255 struct Bernoulli_data *data)
257 int exact = data->options->approx->method == BV_APPROX_NONE;
258 if (poly_u)
259 return poly_u;
260 Vector_Copy(upper+2, row->p, dim+1);
261 value_oppose(tmp, upper[1]);
262 value_addto(row->p[dim], row->p[dim], tmp);
263 return power_sums(faulhaber, data->e, row, tmp, 0, 0, exact);
266 static evalue *compute_poly_l(evalue *poly_l, Value *lower, Vector *row,
267 unsigned dim,
268 struct poly_list *faulhaber,
269 struct Bernoulli_data *data)
271 int exact = data->options->approx->method == BV_APPROX_NONE;
272 if (poly_l)
273 return poly_l;
274 Vector_Copy(lower+2, row->p, dim+1);
275 value_addto(row->p[dim], row->p[dim], lower[1]);
276 return power_sums(faulhaber, data->e, row, lower[1], 0, 1, exact);
279 /* Compute sum of constant term.
281 static evalue *linear_term(const evalue *cst, Value *lower, Value *upper,
282 Vector *row, Value tmp, int exact)
284 evalue *linear;
285 unsigned dim = row->Size - 1;
287 if (EVALUE_IS_ZERO(*cst))
288 return evalue_zero();
290 if (!exact) {
291 value_absolute(tmp, upper[1]);
292 /* upper - lower */
293 Vector_Combine(lower+2, upper+2, row->p, tmp, lower[1], dim+1);
294 value_multiply(tmp, tmp, lower[1]);
295 /* upper - lower + 1 */
296 value_addto(row->p[dim], row->p[dim], tmp);
298 linear = affine2evalue(row->p, tmp, dim);
299 } else {
300 evalue *l;
302 value_absolute(tmp, upper[1]);
303 Vector_Copy(upper+2, row->p, dim+1);
304 value_addto(row->p[dim], row->p[dim], tmp);
305 /* floor(upper+1) */
306 linear = power_sums_base(row, tmp, 0, 1);
308 Vector_Copy(lower+2, row->p, dim+1);
309 /* floor(-lower) */
310 l = power_sums_base(row, lower[1], 0, 1);
312 /* floor(upper+1) + floor(-lower) = floor(upper) - ceil(lower) + 1 */
313 eadd(l, linear);
314 evalue_free(l);
317 emul(cst, linear);
318 return linear;
321 static void Bernoulli_init(unsigned n, void *cb_data)
323 struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
324 int exact = data->options->approx->method == BV_APPROX_NONE;
325 int cases = exact ? 3 : 5;
327 evalue_section_array_ensure(data->sections, cases * n);
330 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data);
333 * This function requires that either the lower or the upper bound
334 * represented by the constraints "lower" and "upper" is an integer
335 * affine expression.
336 * An affine substitution is performed to make this bound exactly
337 * zero, ensuring that in the recursive call to Bernoulli_cb,
338 * only one of the "cases" will apply.
340 static void transform_to_single_case(Matrix *M, Value *lower, Value *upper,
341 struct Bernoulli_data *data)
343 unsigned dim = M->NbColumns-2;
344 Vector *shadow;
345 Value a, b;
346 evalue **subs;
347 const evalue *e = data->e;
348 evalue *t;
349 int i;
351 value_init(a);
352 value_init(b);
353 subs = ALLOCN(evalue *, dim+1);
354 for (i = 0; i < dim; ++i)
355 subs[1+i] = evalue_var(1+i);
356 shadow = Vector_Alloc(2 * (2+dim+1));
357 if (value_one_p(lower[1])) {
358 /* Replace i by i + l' when b = 1 */
359 value_set_si(shadow->p[0], 1);
360 Vector_Oppose(lower+2, shadow->p+1, dim+1);
361 subs[0] = affine2evalue(shadow->p, shadow->p[0], dim+1);
362 /* new lower
363 * i >= 0
364 * new upper
365 * (-a i + u') + a (-l') >= 0
367 value_assign(shadow->p[2+dim+1+1], upper[1]);
368 value_oppose(a, upper[1]);
369 value_set_si(b, 1);
370 Vector_Combine(upper+2, lower+2, shadow->p+2+dim+1+2,
371 b, a, dim+1);
372 upper = shadow->p+2+dim+1;
373 lower = shadow->p;
374 value_set_si(lower[1], 1);
375 Vector_Set(lower+2, 0, dim+1);
376 } else {
377 /* Replace i by i + u' when a = 1 */
378 value_set_si(shadow->p[0], 1);
379 Vector_Copy(upper+2, shadow->p+1, dim+1);
380 subs[0] = affine2evalue(shadow->p, shadow->p[0], dim+1);
381 /* new lower
382 * (b i - l') + b u' >= 0
383 * new upper
384 * -i >= 0
386 value_assign(shadow->p[1], lower[1]);
387 value_set_si(a, 1);
388 value_assign(b, lower[1]);
389 Vector_Combine(upper+2, lower+2, shadow->p+2,
390 b, a, dim+1);
391 upper = shadow->p+2+dim+1;
392 lower = shadow->p;
393 value_set_si(upper[1], -1);
394 Vector_Set(upper+2, 0, dim+1);
396 value_clear(a);
397 value_clear(b);
399 t = evalue_dup(data->e);
400 evalue_substitute(t, subs);
401 reduce_evalue(t);
402 data->e = t;
403 for (i = 0; i < dim+1; ++i)
404 evalue_free(subs[i]);
405 free(subs);
407 Bernoulli_cb(M, lower, upper, data);
409 evalue_free(t);
410 data->e = e;
411 Vector_Free(shadow);
414 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
416 struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
417 struct evalue_section_array *sections = data->sections;
418 Matrix *M2;
419 Polyhedron *T;
420 const evalue *factor = NULL;
421 evalue *linear = NULL;
422 int constant = 0;
423 Value tmp;
424 unsigned dim = M->NbColumns-2;
425 Vector *row;
426 int exact = data->options->approx->method == BV_APPROX_NONE;
427 int cases = exact ? 3 : 5;
429 assert(lower);
430 assert(upper);
431 evalue_section_array_ensure(sections, sections->ns + cases);
433 M2 = Matrix_Copy(M);
434 T = Constraints2Polyhedron(M2, data->options->MaxRays);
435 Matrix_Free(M2);
437 POL_ENSURE_VERTICES(T);
438 if (emptyQ(T)) {
439 Polyhedron_Free(T);
440 return;
443 constant = value_notzero_p(data->e->d) ||
444 data->e->x.p->type != polynomial ||
445 data->e->x.p->pos != 1;
446 if (!constant && (value_one_p(lower[1]) || value_mone_p(upper[1]))) {
447 int single_case;
448 int lower_cst, upper_cst;
450 lower_cst = First_Non_Zero(lower+2, dim) == -1;
451 upper_cst = First_Non_Zero(upper+2, dim) == -1;
452 single_case =
453 (lower_cst && value_negz_p(lower[2+dim])) ||
454 (upper_cst && value_negz_p(upper[2+dim])) ||
455 (lower_cst && upper_cst &&
456 value_posz_p(lower[2+dim]) && value_posz_p(upper[2+dim]));
457 if (!single_case) {
458 transform_to_single_case(M, lower, upper, data);
459 Polyhedron_Free(T);
460 return;
464 assert(lower != upper);
466 row = Vector_Alloc(dim+1);
467 value_init(tmp);
468 if (value_notzero_p(data->e->d)) {
469 factor = data->e;
470 constant = 1;
471 } else {
472 if (data->e->x.p->type == polynomial && data->e->x.p->pos == 1)
473 factor = shifted_copy(&data->e->x.p->arr[0]);
474 else {
475 factor = shifted_copy(data->e);
476 constant = 1;
479 linear = linear_term(factor, lower, upper, row, tmp, exact);
481 if (constant) {
482 evalue_section_array_add(sections, T, linear);
483 data->options->stats->bernoulli_sums++;
484 } else {
485 evalue *poly_u = NULL, *poly_l = NULL;
486 Polyhedron *D;
487 struct poly_list *faulhaber;
488 assert(data->e->x.p->type == polynomial);
489 assert(data->e->x.p->pos == 1);
490 faulhaber = faulhaber_compute(data->e->x.p->size-1);
491 /* lower is the constraint
492 * b i - l' >= 0 i >= l'/b = l
493 * upper is the constraint
494 * -a i + u' >= 0 i <= u'/a = u
496 M2 = Matrix_Alloc(3, 2+T->Dimension);
497 value_set_si(M2->p[0][0], 1);
498 value_set_si(M2->p[1][0], 1);
499 value_set_si(M2->p[2][0], 1);
500 /* Case 1:
501 * 1 <= l l' - b >= 0
502 * 0 < l l' - 1 >= 0 if exact
504 if (exact)
505 bound_constraint(M2->p[0]+1, T->Dimension, lower+1, -1, 0, 1);
506 else
507 bound_constraint(M2->p[0]+1, T->Dimension, lower+1, -1, -1, 0);
508 D = AddConstraints(M2->p_Init, 1, T, data->options->MaxRays);
509 POL_ENSURE_VERTICES(D);
510 if (emptyQ2(D))
511 Polyhedron_Free(D);
512 else {
513 evalue *extra;
514 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
515 faulhaber, data);
516 Vector_Oppose(lower+2, row->p, dim+1);
517 extra = power_sums(faulhaber, data->e, row, lower[1], 1, 0, exact);
518 eadd(poly_u, extra);
519 eadd(linear, extra);
521 evalue_section_array_add(sections, D, extra);
522 data->options->stats->bernoulli_sums++;
525 /* Case 2:
526 * 1 <= -u -u' - a >= 0
527 * 0 < -u -u' - 1 >= 0 if exact
529 if (exact)
530 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, -1, 0, 1);
531 else
532 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, -1, 1, 0);
533 D = AddConstraints(M2->p_Init, 1, T, data->options->MaxRays);
534 POL_ENSURE_VERTICES(D);
535 if (emptyQ2(D))
536 Polyhedron_Free(D);
537 else {
538 evalue *extra;
539 poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
540 Vector_Oppose(upper+2, row->p, dim+1);
541 value_oppose(tmp, upper[1]);
542 extra = power_sums(faulhaber, data->e, row, tmp, 1, 1, exact);
543 eadd(poly_l, extra);
544 eadd(linear, extra);
546 evalue_section_array_add(sections, D, extra);
547 data->options->stats->bernoulli_sums++;
550 /* Case 3:
551 * u >= 0 u' >= 0
552 * -l >= 0 -l' >= 0
554 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, 0, 0);
555 bound_constraint(M2->p[1]+1, T->Dimension, lower+1, 1, 0, 0);
556 D = AddConstraints(M2->p_Init, 2, T, data->options->MaxRays);
557 POL_ENSURE_VERTICES(D);
558 if (emptyQ2(D))
559 Polyhedron_Free(D);
560 else {
561 evalue *e;
562 poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
563 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
564 faulhaber, data);
566 e = ALLOC(evalue);
567 value_init(e->d);
568 evalue_copy(e, poly_u);
569 eadd(poly_l, e);
570 eadd(linear, e);
571 evalue_section_array_add(sections, D, e);
572 data->options->stats->bernoulli_sums++;
575 if (!exact) {
576 /* Case 4:
577 * l < 1 -l' + b - 1 >= 0
578 * 0 < l l' - 1 >= 0
579 * u >= 1 u' - a >= 0
581 bound_constraint(M2->p[0]+1, T->Dimension, lower+1, 1, 1, 1);
582 bound_constraint(M2->p[1]+1, T->Dimension, lower+1, -1, 0, 1);
583 bound_constraint(M2->p[2]+1, T->Dimension, upper+1, 1, 1, 0);
584 D = AddConstraints(M2->p_Init, 3, T, data->options->MaxRays);
585 POL_ENSURE_VERTICES(D);
586 if (emptyQ2(D))
587 Polyhedron_Free(D);
588 else {
589 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
590 faulhaber, data);
591 eadd(linear, poly_u);
592 evalue_section_array_add(sections, D, poly_u);
593 poly_u = NULL;
594 data->options->stats->bernoulli_sums++;
597 /* Case 5:
598 * -u < 1 u' + a - 1 >= 0
599 * 0 < -u -u' - 1 >= 0
600 * l <= -1 -l' - b >= 0
602 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, -1, 1);
603 bound_constraint(M2->p[1]+1, T->Dimension, upper+1, -1, 0, 1);
604 bound_constraint(M2->p[2]+1, T->Dimension, lower+1, 1, -1, 0);
605 D = AddConstraints(M2->p_Init, 3, T, data->options->MaxRays);
606 POL_ENSURE_VERTICES(D);
607 if (emptyQ2(D))
608 Polyhedron_Free(D);
609 else {
610 poly_l = compute_poly_l(poly_l, lower, row, dim,
611 faulhaber, data);
612 eadd(linear, poly_l);
613 evalue_section_array_add(sections, D, poly_l);
614 poly_l = NULL;
615 data->options->stats->bernoulli_sums++;
619 Matrix_Free(M2);
620 Polyhedron_Free(T);
621 if (poly_l)
622 evalue_free(poly_l);
623 if (poly_u)
624 evalue_free(poly_u);
625 evalue_free(linear);
627 if (factor != data->e)
628 evalue_free((evalue *)factor);
629 value_clear(tmp);
630 Vector_Free(row);
634 * Move the variable at position pos to the front by exchanging
635 * that variable with the first variable.
637 static void more_var_to_front(Polyhedron **P_p , evalue **E_p, int pos)
639 Polyhedron *P = *P_p;
640 evalue *E = *E_p;
641 unsigned dim = P->Dimension;
643 assert(pos != 0);
645 P = Polyhedron_Copy(P);
646 Polyhedron_ExchangeColumns(P, 1, 1+pos);
647 *P_p = P;
649 if (value_zero_p(E->d)) {
650 int j;
651 evalue **subs;
653 subs = ALLOCN(evalue *, dim);
654 for (j = 0; j < dim; ++j)
655 subs[j] = evalue_var(j);
656 E = subs[0];
657 subs[0] = subs[pos];
658 subs[pos] = E;
659 E = evalue_dup(*E_p);
660 evalue_substitute(E, subs);
661 for (j = 0; j < dim; ++j)
662 evalue_free(subs[j]);
663 free(subs);
664 *E_p = E;
668 static int type_offset(enode *p)
670 return p->type == fractional ? 1 :
671 p->type == flooring ? 1 : 0;
674 static void adjust_periods(evalue *e, unsigned nvar, Vector *periods)
676 for (; value_zero_p(e->d); e = &e->x.p->arr[0]) {
677 int pos;
678 assert(e->x.p->type == polynomial);
679 assert(e->x.p->size == 2);
680 assert(value_notzero_p(e->x.p->arr[1].d));
682 pos = e->x.p->pos - 1;
683 if (pos >= nvar)
684 break;
686 value_lcm(periods->p[pos], periods->p[pos], e->x.p->arr[1].d);
690 static void fractional_periods_r(evalue *e, unsigned nvar, Vector *periods)
692 int i;
694 if (value_notzero_p(e->d))
695 return;
697 assert(e->x.p->type == polynomial || e->x.p->type == fractional);
699 if (e->x.p->type == fractional)
700 adjust_periods(&e->x.p->arr[0], nvar, periods);
702 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
703 fractional_periods_r(&e->x.p->arr[i], nvar, periods);
707 * For each of the nvar variables, compute the lcm of the
708 * denominators of the coefficients of that variable in
709 * any of the fractional parts.
711 static Vector *fractional_periods(evalue *e, unsigned nvar)
713 int i;
714 Vector *periods = Vector_Alloc(nvar);
716 for (i = 0; i < nvar; ++i)
717 value_set_si(periods->p[i], 1);
719 fractional_periods_r(e, nvar, periods);
721 return periods;
724 /* Move "best" variable to sum over into the first position,
725 * possibly changing *P_p and *E_p.
727 * If there are any fractional parts (period != NULL), then we
728 * first look for a variable with the smallest lcm of denominators
729 * inside a fractional part. This denominator is assigned to period
730 * and corresponds to the number of "splinters" we need to construct
731 * at this level.
733 * Of those with this denominator (all if period == NULL or if there
734 * are no fractional parts), we select the variable with the smallest
735 * maximal coefficient, as this coefficient will become a denominator
736 * in new fractional parts (for an exact computation), which may
737 * lead to splintering in the next step.
739 static void move_best_to_front(Polyhedron **P_p, evalue **E_p, unsigned nvar,
740 Value *period)
742 Polyhedron *P = *P_p;
743 evalue *E = *E_p;
744 int i, j, best_i = -1;
745 Vector *periods;
746 Value best, max;
748 assert(nvar >= 1);
750 if (period) {
751 periods = fractional_periods(E, nvar);
752 value_assign(*period, periods->p[0]);
753 for (i = 1; i < nvar; ++i)
754 if (value_lt(periods->p[i], *period))
755 value_assign(*period, periods->p[i]);
758 value_init(best);
759 value_init(max);
761 for (i = 0; i < nvar; ++i) {
762 if (period && value_ne(*period, periods->p[i]))
763 continue;
765 value_set_si(max, 0);
767 for (j = 0; j < P->NbConstraints; ++j) {
768 if (value_zero_p(P->Constraint[j][1+i]))
769 continue;
770 if (First_Non_Zero(P->Constraint[j]+1, i) == -1 &&
771 First_Non_Zero(P->Constraint[j]+1+i+1, nvar-i-1) == -1)
772 continue;
773 if (value_abs_gt(P->Constraint[j][1+i], max))
774 value_absolute(max, P->Constraint[j][1+i]);
777 if (best_i == -1 || value_lt(max, best)) {
778 value_assign(best, max);
779 best_i = i;
783 value_clear(best);
784 value_clear(max);
786 if (period)
787 Vector_Free(periods);
789 if (best_i != 0)
790 more_var_to_front(P_p, E_p, best_i);
792 return;
795 static evalue *sum_over_polytope_base(Polyhedron *P, evalue *E, unsigned nvar,
796 struct evalue_section_array *sections,
797 struct barvinok_options *options)
799 evalue *res;
800 struct Bernoulli_data data;
802 assert(P->NbEq == 0);
804 sections->ns = 0;
805 data.options = options;
806 data.sections = sections;
807 data.e = E;
809 for_each_lower_upper_bound(P, Bernoulli_init, Bernoulli_cb, &data);
811 res = evalue_from_section_array(sections->s, sections->ns);
813 if (nvar > 1) {
814 evalue *tmp = barvinok_summate(res, nvar-1, options);
815 evalue_free(res);
816 res = tmp;
819 return res;
822 /* Splinter P into period slices along the first variable x = period y + c,
823 * 0 <= c < period, ensuring no fractional part contains the first variable,
824 * and sum over all slices.
826 static evalue *sum_over_polytope_slices(Polyhedron *P, evalue *E,
827 unsigned nvar,
828 Value period,
829 struct evalue_section_array *sections,
830 struct barvinok_options *options)
832 evalue *sum = evalue_zero();
833 evalue **subs;
834 unsigned dim = P->Dimension;
835 Matrix *T;
836 Value i;
837 Value one;
838 int j;
840 value_init(i);
841 value_init(one);
842 value_set_si(one, 1);
844 subs = ALLOCN(evalue *, dim);
846 T = Matrix_Alloc(dim+1, dim+1);
847 value_assign(T->p[0][0], period);
848 for (j = 1; j < dim+1; ++j)
849 value_set_si(T->p[j][j], 1);
851 for (j = 0; j < dim; ++j)
852 subs[j] = evalue_var(j);
853 evalue_mul(subs[0], period);
855 for (value_set_si(i, 0); value_lt(i, period); value_increment(i, i)) {
856 evalue *tmp;
857 Polyhedron *S = DomainPreimage(P, T, options->MaxRays);
858 evalue *e = evalue_dup(E);
859 evalue_substitute(e, subs);
860 reduce_evalue(e);
862 if (S->NbEq)
863 tmp = barvinok_sum_over_polytope(S, e, nvar, sections, options);
864 else
865 tmp = sum_over_polytope_base(S, e, nvar, sections, options);
867 assert(tmp);
868 eadd(tmp, sum);
869 evalue_free(tmp);
871 Domain_Free(S);
872 evalue_free(e);
874 value_increment(T->p[0][dim], T->p[0][dim]);
875 evalue_add_constant(subs[0], one);
878 value_clear(i);
879 value_clear(one);
880 Matrix_Free(T);
881 for (j = 0; j < dim; ++j)
882 evalue_free(subs[j]);
883 free(subs);
885 reduce_evalue(sum);
886 return sum;
889 evalue *bernoulli_summate(Polyhedron *P, evalue *E, unsigned nvar,
890 struct evalue_section_array *sections,
891 struct barvinok_options *options)
893 Polyhedron *P_orig = P;
894 evalue *E_orig = E;
895 Value period;
896 evalue *sum;
897 int exact = options->approx->method == BV_APPROX_NONE;
899 value_init(period);
901 move_best_to_front(&P, &E, nvar, exact ? &period : NULL);
902 if (exact && value_notone_p(period))
903 sum = sum_over_polytope_slices(P, E, nvar, period, sections, options);
904 else
905 sum = sum_over_polytope_base(P, E, nvar, sections, options);
907 if (P != P_orig)
908 Polyhedron_Free(P);
909 if (E != E_orig)
910 evalue_free(E);
912 value_clear(period);
914 return sum;