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"
9
10 #define ALLOC(type) (type*)malloc(sizeof(type))
11 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
12
13 static struct bernoulli_coef bernoulli_coef;
14 static struct poly_list bernoulli;
15 static struct poly_list faulhaber;
16
bernoulli_coef_compute(int n)17 struct bernoulli_coef *bernoulli_coef_compute(int n)
18 {
19 int i, j;
20 Value factor, tmp;
21
22 if (n < bernoulli_coef.n)
23 return &bernoulli_coef;
24
25 if (n >= bernoulli_coef.size) {
26 int size = 3*(n + 5)/2;
27 Vector *b;
28
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);
33 }
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);
39 }
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);
45 }
46 bernoulli_coef.lcm = b;
47
48 bernoulli_coef.size = size;
49 }
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;
58 }
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);
69 }
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);
77 }
78 value_lcm(bernoulli_coef.lcm->p[i],
79 bernoulli_coef.lcm->p[i-1], bernoulli_coef.den->p[i]);
80 }
81 bernoulli_coef.n = n+1;
82 value_clear(factor);
83 value_clear(tmp);
84
85 return &bernoulli_coef;
86 }
87
88 /*
89 * Compute either Bernoulli B_n or Faulhaber F_n polynomials.
90 *
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}
93 */
bernoulli_faulhaber_compute(int n,struct poly_list * pl,int faulhaber)94 static struct poly_list *bernoulli_faulhaber_compute(int n, struct poly_list *pl,
95 int faulhaber)
96 {
97 int i, j;
98 Value factor;
99 struct bernoulli_coef *bc;
100
101 if (n < pl->n)
102 return pl;
103
104 if (n >= pl->size) {
105 int size = 3*(n + 5)/2;
106 Vector **poly;
107
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;
113
114 pl->size = size;
115 }
116
117 bc = bernoulli_coef_compute(n);
118
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);
137 }
138 Vector_Normalize(pl->poly[i]->p, i+faulhaber+2);
139 }
140 value_clear(factor);
141 pl->n = n+1;
142
143 return pl;
144 }
145
bernoulli_compute(int n)146 struct poly_list *bernoulli_compute(int n)
147 {
148 return bernoulli_faulhaber_compute(n, &bernoulli, 0);
149 }
150
faulhaber_compute(int n)151 struct poly_list *faulhaber_compute(int n)
152 {
153 return bernoulli_faulhaber_compute(n, &faulhaber, 1);
154 }
155
shifted_copy(const evalue * src)156 static evalue *shifted_copy(const evalue *src)
157 {
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;
163 }
164
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 }
172 *
173 * Modifies arg (if exact == 1).
174 */
power_sums_base(Vector * arg,Value denom,int neg,int exact)175 static evalue *power_sums_base(Vector *arg, Value denom, int neg, int exact)
176 {
177 evalue *fract;
178 evalue *base = affine2evalue(arg->p, denom, arg->Size-1);
179
180 if (!exact || value_one_p(denom))
181 return base;
182
183 if (neg)
184 Vector_Oppose(arg->p, arg->p, arg->Size);
185
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]);
190 }
191 eadd(fract, base);
192 evalue_free(fract);
193
194 return base;
195 }
196
power_sums(struct poly_list * faulhaber,const evalue * poly,Vector * arg,Value denom,int neg,int alt_neg,int exact)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)
200 {
201 int i;
202 evalue *base = power_sums_base(arg, denom, neg, exact);
203 evalue *sum = evalue_zero();
204
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);
214 }
215 if (neg)
216 evalue_negate(sum);
217 evalue_free(base);
218
219 return sum;
220 }
221
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
227 */
bound_constraint(Value * c,unsigned dim,Value * cst_affine,int sign_affine,int sign_cst,int strict)228 static void bound_constraint(Value *c, unsigned dim,
229 Value *cst_affine,
230 int sign_affine, int sign_cst, int strict)
231 {
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);
236
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]);
241
242 if (strict)
243 value_decrement(c[dim], c[dim]);
244 }
245
246 struct Bernoulli_data {
247 struct barvinok_options *options;
248 struct evalue_section_array *sections;
249 const evalue *e;
250 };
251
compute_poly_u(evalue * poly_u,Value * upper,Vector * row,unsigned dim,Value tmp,struct poly_list * faulhaber,struct Bernoulli_data * data)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)
256 {
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);
264 }
265
compute_poly_l(evalue * poly_l,Value * lower,Vector * row,unsigned dim,struct poly_list * faulhaber,struct Bernoulli_data * data)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)
270 {
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);
277 }
278
279 /* Compute sum of constant term.
280 */
linear_term(const evalue * cst,Value * lower,Value * upper,Vector * row,Value tmp,int exact)281 static evalue *linear_term(const evalue *cst, Value *lower, Value *upper,
282 Vector *row, Value tmp, int exact)
283 {
284 evalue *linear;
285 unsigned dim = row->Size - 1;
286
287 if (EVALUE_IS_ZERO(*cst))
288 return evalue_zero();
289
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);
297
298 linear = affine2evalue(row->p, tmp, dim);
299 } else {
300 evalue *l;
301
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);
307
308 Vector_Copy(lower+2, row->p, dim+1);
309 /* floor(-lower) */
310 l = power_sums_base(row, lower[1], 0, 1);
311
312 /* floor(upper+1) + floor(-lower) = floor(upper) - ceil(lower) + 1 */
313 eadd(l, linear);
314 evalue_free(l);
315 }
316
317 emul(cst, linear);
318 return linear;
319 }
320
Bernoulli_init(unsigned n,void * cb_data)321 static void Bernoulli_init(unsigned n, void *cb_data)
322 {
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;
326
327 evalue_section_array_ensure(data->sections, cases * n);
328 }
329
330 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data);
331
332 /*
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.
339 */
transform_to_single_case(Matrix * M,Value * lower,Value * upper,struct Bernoulli_data * data)340 static void transform_to_single_case(Matrix *M, Value *lower, Value *upper,
341 struct Bernoulli_data *data)
342 {
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;
350
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
366 */
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
385 */
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);
395 }
396 value_clear(a);
397 value_clear(b);
398
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);
406
407 Bernoulli_cb(M, lower, upper, data);
408
409 evalue_free(t);
410 data->e = e;
411 Vector_Free(shadow);
412 }
413
Bernoulli_cb(Matrix * M,Value * lower,Value * upper,void * cb_data)414 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
415 {
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;
428
429 assert(lower);
430 assert(upper);
431 evalue_section_array_ensure(sections, sections->ns + cases);
432
433 M2 = Matrix_Copy(M);
434 T = Constraints2Polyhedron(M2, data->options->MaxRays);
435 Matrix_Free(M2);
436
437 POL_ENSURE_VERTICES(T);
438 if (emptyQ(T)) {
439 Polyhedron_Free(T);
440 return;
441 }
442
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;
449
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;
461 }
462 }
463
464 assert(lower != upper);
465
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;
477 }
478 }
479 linear = linear_term(factor, lower, upper, row, tmp, exact);
480
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
495 */
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
503 */
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);
520
521 evalue_section_array_add(sections, D, extra);
522 data->options->stats->bernoulli_sums++;
523 }
524
525 /* Case 2:
526 * 1 <= -u -u' - a >= 0
527 * 0 < -u -u' - 1 >= 0 if exact
528 */
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);
545
546 evalue_section_array_add(sections, D, extra);
547 data->options->stats->bernoulli_sums++;
548 }
549
550 /* Case 3:
551 * u >= 0 u' >= 0
552 * -l >= 0 -l' >= 0
553 */
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);
565
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++;
573 }
574
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
580 */
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++;
595 }
596
597 /* Case 5:
598 * -u < 1 u' + a - 1 >= 0
599 * 0 < -u -u' - 1 >= 0
600 * l <= -1 -l' - b >= 0
601 */
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++;
616 }
617 }
618
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);
626 }
627 if (factor != data->e)
628 evalue_free((evalue *)factor);
629 value_clear(tmp);
630 Vector_Free(row);
631 }
632
633 /*
634 * Move the variable at position pos to the front by exchanging
635 * that variable with the first variable.
636 */
more_var_to_front(Polyhedron ** P_p,evalue ** E_p,int pos)637 static void more_var_to_front(Polyhedron **P_p , evalue **E_p, int pos)
638 {
639 Polyhedron *P = *P_p;
640 evalue *E = *E_p;
641 unsigned dim = P->Dimension;
642
643 assert(pos != 0);
644
645 P = Polyhedron_Copy(P);
646 Polyhedron_ExchangeColumns(P, 1, 1+pos);
647 *P_p = P;
648
649 if (value_zero_p(E->d)) {
650 int j;
651 evalue **subs;
652
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;
665 }
666 }
667
type_offset(enode * p)668 static int type_offset(enode *p)
669 {
670 return p->type == fractional ? 1 :
671 p->type == flooring ? 1 : 0;
672 }
673
adjust_periods(evalue * e,unsigned nvar,Vector * periods)674 static void adjust_periods(evalue *e, unsigned nvar, Vector *periods)
675 {
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));
681
682 pos = e->x.p->pos - 1;
683 if (pos >= nvar)
684 break;
685
686 value_lcm(periods->p[pos], periods->p[pos], e->x.p->arr[1].d);
687 }
688 }
689
fractional_periods_r(evalue * e,unsigned nvar,Vector * periods)690 static void fractional_periods_r(evalue *e, unsigned nvar, Vector *periods)
691 {
692 int i;
693
694 if (value_notzero_p(e->d))
695 return;
696
697 assert(e->x.p->type == polynomial || e->x.p->type == fractional);
698
699 if (e->x.p->type == fractional)
700 adjust_periods(&e->x.p->arr[0], nvar, periods);
701
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);
704 }
705
706 /*
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.
710 */
fractional_periods(evalue * e,unsigned nvar)711 static Vector *fractional_periods(evalue *e, unsigned nvar)
712 {
713 int i;
714 Vector *periods = Vector_Alloc(nvar);
715
716 for (i = 0; i < nvar; ++i)
717 value_set_si(periods->p[i], 1);
718
719 fractional_periods_r(e, nvar, periods);
720
721 return periods;
722 }
723
724 /* Move "best" variable to sum over into the first position,
725 * possibly changing *P_p and *E_p.
726 *
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.
732 *
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.
738 */
move_best_to_front(Polyhedron ** P_p,evalue ** E_p,unsigned nvar,Value * period)739 static void move_best_to_front(Polyhedron **P_p, evalue **E_p, unsigned nvar,
740 Value *period)
741 {
742 Polyhedron *P = *P_p;
743 evalue *E = *E_p;
744 int i, j, best_i = -1;
745 Vector *periods;
746 Value best, max;
747
748 assert(nvar >= 1);
749
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]);
756 }
757
758 value_init(best);
759 value_init(max);
760
761 for (i = 0; i < nvar; ++i) {
762 if (period && value_ne(*period, periods->p[i]))
763 continue;
764
765 value_set_si(max, 0);
766
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]);
775 }
776
777 if (best_i == -1 || value_lt(max, best)) {
778 value_assign(best, max);
779 best_i = i;
780 }
781 }
782
783 value_clear(best);
784 value_clear(max);
785
786 if (period)
787 Vector_Free(periods);
788
789 if (best_i != 0)
790 more_var_to_front(P_p, E_p, best_i);
791
792 return;
793 }
794
sum_over_polytope_base(Polyhedron * P,evalue * E,unsigned nvar,struct evalue_section_array * sections,struct barvinok_options * options)795 static evalue *sum_over_polytope_base(Polyhedron *P, evalue *E, unsigned nvar,
796 struct evalue_section_array *sections,
797 struct barvinok_options *options)
798 {
799 evalue *res;
800 struct Bernoulli_data data;
801
802 assert(P->NbEq == 0);
803
804 sections->ns = 0;
805 data.options = options;
806 data.sections = sections;
807 data.e = E;
808
809 for_each_lower_upper_bound(P, Bernoulli_init, Bernoulli_cb, &data);
810
811 res = evalue_from_section_array(sections->s, sections->ns);
812
813 if (nvar > 1) {
814 evalue *tmp = barvinok_summate(res, nvar-1, options);
815 evalue_free(res);
816 res = tmp;
817 }
818
819 return res;
820 }
821
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.
825 */
sum_over_polytope_slices(Polyhedron * P,evalue * E,unsigned nvar,Value period,struct evalue_section_array * sections,struct barvinok_options * options)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)
831 {
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;
839
840 value_init(i);
841 value_init(one);
842 value_set_si(one, 1);
843
844 subs = ALLOCN(evalue *, dim);
845
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);
850
851 for (j = 0; j < dim; ++j)
852 subs[j] = evalue_var(j);
853 evalue_mul(subs[0], period);
854
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);
861
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);
866
867 assert(tmp);
868 eadd(tmp, sum);
869 evalue_free(tmp);
870
871 Domain_Free(S);
872 evalue_free(e);
873
874 value_increment(T->p[0][dim], T->p[0][dim]);
875 evalue_add_constant(subs[0], one);
876 }
877
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);
884
885 reduce_evalue(sum);
886 return sum;
887 }
888
bernoulli_summate(Polyhedron * P,evalue * E,unsigned nvar,struct evalue_section_array * sections,struct barvinok_options * options)889 evalue *bernoulli_summate(Polyhedron *P, evalue *E, unsigned nvar,
890 struct evalue_section_array *sections,
891 struct barvinok_options *options)
892 {
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;
898
899 value_init(period);
900
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);
906
907 if (P != P_orig)
908 Polyhedron_Free(P);
909 if (E != E_orig)
910 evalue_free(E);
911
912 value_clear(period);
913
914 return sum;
915 }
916