?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Bernoulli_sum_evalue: optionally handle fractional bounds exactly
[barvinok.git] / bernoulli.c
bloba1dac8b31778f3ba48a5998a72adf9ed28c887b1
1 #include <assert.h>
2 #include <barvinok/options.h>
3 #include <barvinok/util.h>
4 #include "bernoulli.h"
5 #include "lattice_point.h"
6
7 #define ALLOC(type) (type*)malloc(sizeof(type))
8 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
9 #define REALLOCN(ptr,type,n) (type*)realloc(ptr, (n) * sizeof(type))
10
11 static struct bernoulli_coef bernoulli_coef;
12 static struct poly_list bernoulli;
13 static struct poly_list faulhaber;
14
15 struct bernoulli_coef *bernoulli_coef_compute(int n)
16 {
17 int i, j;
18 Value factor, tmp;
19
20 if (n < bernoulli_coef.n)
21 return &bernoulli_coef;
22
23 if (n >= bernoulli_coef.size) {
24 int size = 3*(n + 5)/2;
25 Vector *b;
26
27 b = Vector_Alloc(size);
28 if (bernoulli_coef.n) {
29 Vector_Copy(bernoulli_coef.num->p, b->p, bernoulli_coef.n);
30 Vector_Free(bernoulli_coef.num);
31 }
32 bernoulli_coef.num = b;
33 b = Vector_Alloc(size);
34 if (bernoulli_coef.n) {
35 Vector_Copy(bernoulli_coef.den->p, b->p, bernoulli_coef.n);
36 Vector_Free(bernoulli_coef.den);
37 }
38 bernoulli_coef.den = b;
39 b = Vector_Alloc(size);
40 if (bernoulli_coef.n) {
41 Vector_Copy(bernoulli_coef.lcm->p, b->p, bernoulli_coef.n);
42 Vector_Free(bernoulli_coef.lcm);
43 }
44 bernoulli_coef.lcm = b;
45
46 bernoulli_coef.size = size;
47 }
48 value_init(factor);
49 value_init(tmp);
50 for (i = bernoulli_coef.n; i <= n; ++i) {
51 if (i == 0) {
52 value_set_si(bernoulli_coef.num->p[0], 1);
53 value_set_si(bernoulli_coef.den->p[0], 1);
54 value_set_si(bernoulli_coef.lcm->p[0], 1);
55 continue;
56 }
57 value_set_si(bernoulli_coef.num->p[i], 0);
58 value_set_si(factor, -(i+1));
59 for (j = i-1; j >= 0; --j) {
60 mpz_mul_ui(factor, factor, j+1);
61 mpz_divexact_ui(factor, factor, i+1-j);
62 value_division(tmp, bernoulli_coef.lcm->p[i-1],
63 bernoulli_coef.den->p[j]);
64 value_multiply(tmp, tmp, bernoulli_coef.num->p[j]);
65 value_multiply(tmp, tmp, factor);
66 value_addto(bernoulli_coef.num->p[i], bernoulli_coef.num->p[i], tmp);
67 }
68 mpz_mul_ui(bernoulli_coef.den->p[i], bernoulli_coef.lcm->p[i-1], i+1);
69 value_gcd(tmp, bernoulli_coef.num->p[i], bernoulli_coef.den->p[i]);
70 if (value_notone_p(tmp)) {
71 value_division(bernoulli_coef.num->p[i],
72 bernoulli_coef.num->p[i], tmp);
73 value_division(bernoulli_coef.den->p[i],
74 bernoulli_coef.den->p[i], tmp);
75 }
76 value_lcm(bernoulli_coef.lcm->p[i],
77 bernoulli_coef.lcm->p[i-1], bernoulli_coef.den->p[i]);
78 }
79 bernoulli_coef.n = n+1;
80 value_clear(factor);
81 value_clear(tmp);
82
83 return &bernoulli_coef;
84 }
85
86 /*
87 * Compute either Bernoulli B_n or Faulhaber F_n polynomials.
88 *
89 * B_n = sum_{k=0}^n { n \choose k } b_k x^{n-k}
90 * F_n = 1/(n+1) sum_{k=0}^n { n+1 \choose k } b_k x^{n+1-k}
91 */
92 static struct poly_list *bernoulli_faulhaber_compute(int n, struct poly_list *pl,
93 int faulhaber)
94 {
95 int i, j;
96 Value factor;
97 struct bernoulli_coef *bc;
98
99 if (n < pl->n)
100 return pl;
101
102 if (n >= pl->size) {
103 int size = 3*(n + 5)/2;
104 Vector **poly;
105
106 poly = ALLOCN(Vector *, size);
107 for (i = 0; i < pl->n; ++i)
108 poly[i] = pl->poly[i];
109 free(pl->poly);
110 pl->poly = poly;
111
112 pl->size = size;
113 }
114
115 bc = bernoulli_coef_compute(n);
116
117 value_init(factor);
118 for (i = pl->n; i <= n; ++i) {
119 pl->poly[i] = Vector_Alloc(i+faulhaber+2);
120 value_assign(pl->poly[i]->p[i+faulhaber], bc->lcm->p[i]);
121 if (faulhaber)
122 mpz_mul_ui(pl->poly[i]->p[i+2], bc->lcm->p[i], i+1);
123 else
124 value_assign(pl->poly[i]->p[i+1], bc->lcm->p[i]);
125 value_set_si(factor, 1);
126 for (j = 1; j <= i; ++j) {
127 mpz_mul_ui(factor, factor, i+faulhaber+1-j);
128 mpz_divexact_ui(factor, factor, j);
129 value_division(pl->poly[i]->p[i+faulhaber-j],
130 bc->lcm->p[i], bc->den->p[j]);
131 value_multiply(pl->poly[i]->p[i+faulhaber-j],
132 pl->poly[i]->p[i+faulhaber-j], bc->num->p[j]);
133 value_multiply(pl->poly[i]->p[i+faulhaber-j],
134 pl->poly[i]->p[i+faulhaber-j], factor);
135 }
136 Vector_Normalize(pl->poly[i]->p, i+faulhaber+2);
137 }
138 value_clear(factor);
139 pl->n = n+1;
140
141 return pl;
142 }
143
144 struct poly_list *bernoulli_compute(int n)
145 {
146 return bernoulli_faulhaber_compute(n, &bernoulli, 0);
147 }
148
149 struct poly_list *faulhaber_compute(int n)
150 {
151 return bernoulli_faulhaber_compute(n, &faulhaber, 1);
152 }
153
154 /* shift variables in polynomial one down */
155 static void shift(evalue *e)
156 {
157 int i;
158 if (value_notzero_p(e->d))
159 return;
160 assert(e->x.p->type == polynomial || e->x.p->type == fractional);
161 if (e->x.p->type == polynomial) {
162 assert(e->x.p->pos > 1);
163 e->x.p->pos--;
164 }
165 for (i = 0; i < e->x.p->size; ++i)
166 shift(&e->x.p->arr[i]);
167 }
168
169 /* shift variables in polynomial n up */
170 static void unshift(evalue *e, unsigned n)
171 {
172 int i;
173 if (value_notzero_p(e->d))
174 return;
175 assert(e->x.p->type == polynomial || e->x.p->type == fractional);
176 if (e->x.p->type == polynomial)
177 e->x.p->pos += n;
178 for (i = 0; i < e->x.p->size; ++i)
179 unshift(&e->x.p->arr[i], n);
180 }
181
182 static evalue *shifted_copy(evalue *src)
183 {
184 evalue *e = ALLOC(evalue);
185 value_init(e->d);
186 evalue_copy(e, src);
187 shift(e);
188 return e;
189 }
190
191 /* Computes the argument for the Faulhaber polynomial.
192 * If we are computing a polynomial approximation (exact == 0),
193 * then this is simply arg/denom.
194 * Otherwise, if neg == 0, then we are dealing with an upper bound
195 * and we need to compute floor(arg/denom) = arg/denom - { arg/denom }
196 * If neg == 1, then we are dealing with a lower bound
197 * and we need to compute ceil(arg/denom) = arg/denom + { -arg/denom }
198 *
199 * Modifies arg (if exact == 1).
200 */
201 static evalue *power_sums_base(Vector *arg, Value denom, int neg, int exact)
202 {
203 evalue *fract;
204 evalue *base = affine2evalue(arg->p, denom, arg->Size-1);
205
206 if (!exact || value_one_p(denom))
207 return base;
208
209 if (neg)
210 Vector_Oppose(arg->p, arg->p, arg->Size);
211
212 fract = fractional_part(arg->p, denom, arg->Size-1, NULL);
213 if (!neg) {
214 value_set_si(arg->p[0], -1);
215 evalue_mul(fract, arg->p[0]);
216 }
217 eadd(fract, base);
218 evalue_free(fract);
219
220 return base;
221 }
222
223 static evalue *power_sums(struct poly_list *faulhaber, evalue *poly,
224 Vector *arg, Value denom, int neg, int alt_neg,
225 int exact)
226 {
227 int i;
228 evalue *base = power_sums_base(arg, denom, neg, exact);
229 evalue *sum = evalue_zero();
230
231 for (i = 1; i < poly->x.p->size; ++i) {
232 evalue *term = evalue_polynomial(faulhaber->poly[i], base);
233 evalue *factor = shifted_copy(&poly->x.p->arr[i]);
234 emul(factor, term);
235 if (alt_neg && (i % 2))
236 evalue_negate(term);
237 eadd(term, sum);
238 evalue_free(factor);
239 evalue_free(term);
240 }
241 if (neg)
242 evalue_negate(sum);
243 evalue_free(base);
244
245 return sum;
246 }
247
248 /* Given a constraint (cst_affine) a x + b y + c >= 0, compate a constraint (c)
249 * +- (b y + c) +- a >=,> 0
250 * ^ ^ ^
251 * | | strict
252 * sign_affine sign_cst
253 */
254 static void bound_constraint(Value *c, unsigned dim,
255 Value *cst_affine,
256 int sign_affine, int sign_cst, int strict)
257 {
258 if (sign_affine == -1)
259 Vector_Oppose(cst_affine+1, c, dim+1);
260 else
261 Vector_Copy(cst_affine+1, c, dim+1);
262
263 if (sign_cst == -1)
264 value_subtract(c[dim], c[dim], cst_affine[0]);
265 else if (sign_cst == 1)
266 value_addto(c[dim], c[dim], cst_affine[0]);
267
268 if (strict)
269 value_decrement(c[dim], c[dim]);
270 }
271
272 struct Bernoulli_data {
273 struct barvinok_options *options;
274 struct evalue_section *s;
275 int size;
276 int ns;
277 evalue *e;
278 };
279
280 static evalue *compute_poly_u(evalue *poly_u, Value *upper, Vector *row,
281 unsigned dim, Value tmp,
282 struct poly_list *faulhaber,
283 struct Bernoulli_data *data)
284 {
285 int exact = data->options->approximation_method == BV_APPROX_NONE;
286 if (poly_u)
287 return poly_u;
288 Vector_Copy(upper+2, row->p, dim+1);
289 value_oppose(tmp, upper[1]);
290 value_addto(row->p[dim], row->p[dim], tmp);
291 return power_sums(faulhaber, data->e, row, tmp, 0, 0, exact);
292 }
293
294 static evalue *compute_poly_l(evalue *poly_l, Value *lower, Vector *row,
295 unsigned dim,
296 struct poly_list *faulhaber,
297 struct Bernoulli_data *data)
298 {
299 int exact = data->options->approximation_method == BV_APPROX_NONE;
300 if (poly_l)
301 return poly_l;
302 Vector_Copy(lower+2, row->p, dim+1);
303 value_addto(row->p[dim], row->p[dim], lower[1]);
304 return power_sums(faulhaber, data->e, row, lower[1], 0, 1, exact);
305 }
306
307 /* Compute sum of constant term.
308 */
309 static evalue *linear_term(evalue *cst, Value *lower, Value *upper,
310 Vector *row, Value tmp, int exact)
311 {
312 evalue *linear;
313 unsigned dim = row->Size - 1;
314
315 if (EVALUE_IS_ZERO(*cst))
316 return evalue_zero();
317
318 if (!exact) {
319 value_absolute(tmp, upper[1]);
320 /* upper - lower */
321 Vector_Combine(lower+2, upper+2, row->p, tmp, lower[1], dim+1);
322 value_multiply(tmp, tmp, lower[1]);
323 /* upper - lower + 1 */
324 value_addto(row->p[dim], row->p[dim], tmp);
325
326 linear = affine2evalue(row->p, tmp, dim);
327 } else {
328 evalue *l;
329
330 value_absolute(tmp, upper[1]);
331 Vector_Copy(upper+2, row->p, dim+1);
332 value_addto(row->p[dim], row->p[dim], tmp);
333 /* floor(upper+1) */
334 linear = power_sums_base(row, tmp, 0, 1);
335
336 Vector_Copy(lower+2, row->p, dim+1);
337 /* floor(-lower) */
338 l = power_sums_base(row, lower[1], 0, 1);
339
340 /* floor(upper+1) + floor(-lower) = floor(upper) - ceil(lower) + 1 */
341 eadd(l, linear);
342 evalue_free(l);
343 }
344
345 emul(cst, linear);
346 return linear;
347 }
348
349 static void Bernoulli_init(unsigned n, void *cb_data)
350 {
351 struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
352 int exact = data->options->approximation_method == BV_APPROX_NONE;
353 int cases = exact ? 3 : 5;
354
355 if (cases * n <= data->size)
356 return;
357
358 data->size = cases * (n + 16);
359 data->s = REALLOCN(data->s, struct evalue_section, data->size);
360 }
361
362 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
363 {
364 struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
365 Matrix *M2;
366 Polyhedron *T;
367 evalue *factor = NULL;
368 evalue *linear = NULL;
369 int constant = 0;
370 Value tmp;
371 unsigned dim = M->NbColumns-2;
372 Vector *row;
373 int exact = data->options->approximation_method == BV_APPROX_NONE;
374 int cases = exact ? 3 : 5;
375
376 assert(lower);
377 assert(upper);
378 assert(data->ns + cases <= data->size);
379
380 M2 = Matrix_Copy(M);
381 T = Constraints2Polyhedron(M2, data->options->MaxRays);
382 Matrix_Free(M2);
383
384 POL_ENSURE_VERTICES(T);
385 if (emptyQ(T)) {
386 Polyhedron_Free(T);
387 return;
388 }
389
390 assert(lower != upper);
391
392 row = Vector_Alloc(dim+1);
393 value_init(tmp);
394 if (value_notzero_p(data->e->d)) {
395 factor = data->e;
396 constant = 1;
397 } else {
398 if (data->e->x.p->type == polynomial && data->e->x.p->pos == 1)
399 factor = shifted_copy(&data->e->x.p->arr[0]);
400 else {
401 factor = shifted_copy(data->e);
402 constant = 1;
403 }
404 }
405 linear = linear_term(factor, lower, upper, row, tmp, exact);
406
407 if (constant) {
408 data->s[data->ns].E = linear;
409 data->s[data->ns].D = T;
410 ++data->ns;
411 } else {
412 evalue *poly_u = NULL, *poly_l = NULL;
413 Polyhedron *D;
414 struct poly_list *faulhaber;
415 assert(data->e->x.p->type == polynomial);
416 assert(data->e->x.p->pos == 1);
417 faulhaber = faulhaber_compute(data->e->x.p->size-1);
418 /* lower is the constraint
419 * b i - l' >= 0 i >= l'/b = l
420 * upper is the constraint
421 * -a i + u' >= 0 i <= u'/a = u
422 */
423 M2 = Matrix_Alloc(3, 2+T->Dimension);
424 value_set_si(M2->p[0][0], 1);
425 value_set_si(M2->p[1][0], 1);
426 value_set_si(M2->p[2][0], 1);
427 /* Case 1:
428 * 1 <= l l' - b >= 0
429 * 0 < l l' - 1 >= 0 if exact
430 */
431 if (exact)
432 bound_constraint(M2->p[0]+1, T->Dimension, lower+1, -1, 0, 1);
433 else
434 bound_constraint(M2->p[0]+1, T->Dimension, lower+1, -1, -1, 0);
435 D = AddConstraints(M2->p_Init, 1, T, data->options->MaxRays);
436 POL_ENSURE_VERTICES(D);
437 if (emptyQ2(D))
438 Polyhedron_Free(D);
439 else {
440 evalue *extra;
441 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
442 faulhaber, data);
443 Vector_Oppose(lower+2, row->p, dim+1);
444 extra = power_sums(faulhaber, data->e, row, lower[1], 1, 0, exact);
445 eadd(poly_u, extra);
446 eadd(linear, extra);
447
448 data->s[data->ns].E = extra;
449 data->s[data->ns].D = D;
450 ++data->ns;
451 }
452
453 /* Case 2:
454 * 1 <= -u -u' - a >= 0
455 * 0 < -u -u' - 1 >= 0 if exact
456 */
457 if (exact)
458 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, -1, 0, 1);
459 else
460 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, -1, 1, 0);
461 D = AddConstraints(M2->p_Init, 1, T, data->options->MaxRays);
462 POL_ENSURE_VERTICES(D);
463 if (emptyQ2(D))
464 Polyhedron_Free(D);
465 else {
466 evalue *extra;
467 poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
468 Vector_Oppose(upper+2, row->p, dim+1);
469 value_oppose(tmp, upper[1]);
470 extra = power_sums(faulhaber, data->e, row, tmp, 1, 1, exact);
471 eadd(poly_l, extra);
472 eadd(linear, extra);
473
474 data->s[data->ns].E = extra;
475 data->s[data->ns].D = D;
476 ++data->ns;
477 }
478
479 /* Case 3:
480 * u >= 0 u' >= 0
481 * -l >= 0 -l' >= 0
482 */
483 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, 0, 0);
484 bound_constraint(M2->p[1]+1, T->Dimension, lower+1, 1, 0, 0);
485 D = AddConstraints(M2->p_Init, 2, T, data->options->MaxRays);
486 POL_ENSURE_VERTICES(D);
487 if (emptyQ2(D))
488 Polyhedron_Free(D);
489 else {
490 poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
491 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
492 faulhaber, data);
493
494 data->s[data->ns].E = ALLOC(evalue);
495 value_init(data->s[data->ns].E->d);
496 evalue_copy(data->s[data->ns].E, poly_u);
497 eadd(poly_l, data->s[data->ns].E);
498 eadd(linear, data->s[data->ns].E);
499 data->s[data->ns].D = D;
500 ++data->ns;
501 }
502
503 if (!exact) {
504 /* Case 4:
505 * l < 1 -l' + b - 1 >= 0
506 * 0 < l l' - 1 >= 0
507 * u >= 1 u' - a >= 0
508 */
509 bound_constraint(M2->p[0]+1, T->Dimension, lower+1, 1, 1, 1);
510 bound_constraint(M2->p[1]+1, T->Dimension, lower+1, -1, 0, 1);
511 bound_constraint(M2->p[2]+1, T->Dimension, upper+1, 1, 1, 0);
512 D = AddConstraints(M2->p_Init, 3, T, data->options->MaxRays);
513 POL_ENSURE_VERTICES(D);
514 if (emptyQ2(D))
515 Polyhedron_Free(D);
516 else {
517 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
518 faulhaber, data);
519 eadd(linear, poly_u);
520 data->s[data->ns].E = poly_u;
521 poly_u = NULL;
522 data->s[data->ns].D = D;
523 ++data->ns;
524 }
525
526 /* Case 5:
527 * -u < 1 u' + a - 1 >= 0
528 * 0 < -u -u' - 1 >= 0
529 * l <= -1 -l' - b >= 0
530 */
531 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, -1, 1);
532 bound_constraint(M2->p[1]+1, T->Dimension, upper+1, -1, 0, 1);
533 bound_constraint(M2->p[2]+1, T->Dimension, lower+1, 1, -1, 0);
534 D = AddConstraints(M2->p_Init, 3, T, data->options->MaxRays);
535 POL_ENSURE_VERTICES(D);
536 if (emptyQ2(D))
537 Polyhedron_Free(D);
538 else {
539 poly_l = compute_poly_l(poly_l, lower, row, dim,
540 faulhaber, data);
541 eadd(linear, poly_l);
542 data->s[data->ns].E = poly_l;
543 poly_l = NULL;
544 data->s[data->ns].D = D;
545 ++data->ns;
546 }
547 }
548
549 Matrix_Free(M2);
550 Polyhedron_Free(T);
551 if (poly_l)
552 evalue_free(poly_l);
553 if (poly_u)
554 evalue_free(poly_u);
555 evalue_free(linear);
556 }
557 if (factor != data->e)
558 evalue_free(factor);
559 value_clear(tmp);
560 Vector_Free(row);
561 }
562
563 /*
564 * Move the variable at position pos to the front by exchanging
565 * that variable with the first variable.
566 */
567 static void more_var_to_front(Polyhedron **P_p , evalue **E_p, int pos)
568 {
569 Polyhedron *P = *P_p;
570 evalue *E = *E_p;
571 unsigned dim = P->Dimension;
572
573 assert(pos != 0);
574
575 P = Polyhedron_Copy(P);
576 Polyhedron_ExchangeColumns(P, 1, 1+pos);
577 *P_p = P;
578
579 if (value_zero_p(E->d)) {
580 int j;
581 evalue **subs;
582
583 subs = ALLOCN(evalue *, dim);
584 for (j = 0; j < dim; ++j)
585 subs[j] = evalue_var(j);
586 E = subs[0];
587 subs[0] = subs[pos];
588 subs[pos] = E;
589 E = evalue_dup(*E_p);
590 evalue_substitute(E, subs);
591 for (j = 0; j < dim; ++j)
592 evalue_free(subs[j]);
593 free(subs);
594 *E_p = E;
595 }
596 }
597
598 static int type_offset(enode *p)
599 {
600 return p->type == fractional ? 1 :
601 p->type == flooring ? 1 : 0;
602 }
603
604 static void adjust_periods(evalue *e, unsigned nvar, Vector *periods)
605 {
606 for (; value_zero_p(e->d); e = &e->x.p->arr[0]) {
607 int pos;
608 assert(e->x.p->type == polynomial);
609 assert(e->x.p->size == 2);
610 assert(value_notzero_p(e->x.p->arr[1].d));
611
612 pos = e->x.p->pos - 1;
613 if (pos >= nvar)
614 break;
615
616 value_lcm(periods->p[pos], periods->p[pos], e->x.p->arr[1].d);
617 }
618 }
619
620 static void fractional_periods_r(evalue *e, unsigned nvar, Vector *periods)
621 {
622 int i;
623
624 if (value_notzero_p(e->d))
625 return;
626
627 assert(e->x.p->type == polynomial || e->x.p->type == fractional);
628
629 if (e->x.p->type == fractional)
630 adjust_periods(&e->x.p->arr[0], nvar, periods);
631
632 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
633 fractional_periods_r(&e->x.p->arr[i], nvar, periods);
634 }
635
636 /*
637 * For each of the nvar variables, compute the lcm of the
638 * denominators of the coefficients of that variable in
639 * any of the fractional parts.
640 */
641 static Vector *fractional_periods(evalue *e, unsigned nvar)
642 {
643 int i;
644 Vector *periods = Vector_Alloc(nvar);
645
646 for (i = 0; i < nvar; ++i)
647 value_set_si(periods->p[i], 1);
648
649 fractional_periods_r(e, nvar, periods);
650
651 return periods;
652 }
653
654 /* Move "best" variable to sum over into the first position,
655 * possibly changing *P_p and *E_p.
656 *
657 * If there are any fractional parts (period != NULL), then we
658 * first look for a variable with the smallest lcm of denominators
659 * inside a fractional part. This denominator is assigned to period
660 * and corresponds to the number of "splinters" we need to construct
661 * at this level.
662 *
663 * Of those with this denominator (all if period == NULL or if there
664 * are no fractional parts), we select the variable with the smallest
665 * maximal coefficient, as this coefficient will become a denominator
666 * in new fractional parts (for an exact computation), which may
667 * lead to splintering in the next step.
668 */
669 static void move_best_to_front(Polyhedron **P_p, evalue **E_p, unsigned nvar,
670 Value *period)
671 {
672 Polyhedron *P = *P_p;
673 evalue *E = *E_p;
674 int i, j, best_i = -1;
675 Vector *periods;
676 Value best, max;
677
678 assert(nvar >= 1);
679
680 if (period) {
681 periods = fractional_periods(E, nvar);
682 value_assign(*period, periods->p[0]);
683 for (i = 1; i < nvar; ++i)
684 if (value_lt(periods->p[i], *period))
685 value_assign(*period, periods->p[i]);
686 }
687
688 value_init(best);
689 value_init(max);
690
691 for (i = 0; i < nvar; ++i) {
692 if (period && value_ne(*period, periods->p[i]))
693 continue;
694
695 value_set_si(max, 0);
696
697 for (j = 0; j < P->NbConstraints; ++j) {
698 if (value_zero_p(P->Constraint[j][1+i]))
699 continue;
700 if (First_Non_Zero(P->Constraint[j]+1, i) == -1 &&
701 First_Non_Zero(P->Constraint[j]+1+i+1, nvar-i-1) == -1)
702 continue;
703 if (value_abs_gt(P->Constraint[j][1+i], max))
704 value_absolute(max, P->Constraint[j][1+i]);
705 }
706
707 if (best_i == -1 || value_lt(max, best)) {
708 value_assign(best, max);
709 best_i = i;
710 }
711 }
712
713 value_clear(best);
714 value_clear(max);
715
716 if (period)
717 Vector_Free(periods);
718
719 if (best_i != 0)
720 more_var_to_front(P_p, E_p, best_i);
721
722 return;
723 }
724
725 static evalue *sum_over_polytope_base(Polyhedron *P, evalue *E, unsigned nvar,
726 struct Bernoulli_data *data,
727 struct barvinok_options *options)
728 {
729 evalue *res;
730
731 assert(P->NbEq == 0);
732
733 data->ns = 0;
734 data->e = E;
735
736 for_each_lower_upper_bound(P, Bernoulli_init, Bernoulli_cb, data);
737
738 res = evalue_from_section_array(data->s, data->ns);
739
740 if (nvar > 1) {
741 evalue *tmp = Bernoulli_sum_evalue(res, nvar-1, options);
742 evalue_free(res);
743 res = tmp;
744 }
745
746 return res;
747 }
748
749 static evalue *sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
750 struct Bernoulli_data *data,
751 struct barvinok_options *options);
752
753 static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E,
754 unsigned nvar,
755 struct Bernoulli_data *data,
756 struct barvinok_options *options)
757 {
758 unsigned dim = P->Dimension;
759 unsigned new_dim, new_nparam;
760 Matrix *T = NULL, *CP = NULL;
761 evalue **subs;
762 evalue *sum;
763 int j;
764
765 if (emptyQ(P))
766 return evalue_zero();
767
768 assert(P->NbEq > 0);
769
770 remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays);
771
772 if (emptyQ(P)) {
773 Polyhedron_Free(P);
774 return evalue_zero();
775 }
776
777 new_nparam = CP ? CP->NbColumns-1 : dim - nvar;
778 new_dim = T ? T->NbColumns-1 : nvar + new_nparam;
779
780 /* We can avoid these substitutions if E is a constant */
781 subs = ALLOCN(evalue *, dim);
782 for (j = 0; j < nvar; ++j) {
783 if (T)
784 subs[j] = affine2evalue(T->p[j], T->p[nvar+new_nparam][new_dim],
785 new_dim);
786 else
787 subs[j] = evalue_var(j);
788 }
789 for (j = 0; j < dim-nvar; ++j) {
790 if (CP)
791 subs[nvar+j] = affine2evalue(CP->p[j], CP->p[dim-nvar][new_nparam],
792 new_nparam);
793 else
794 subs[nvar+j] = evalue_var(j);
795 unshift(subs[nvar+j], new_dim-new_nparam);
796 }
797
798 E = evalue_dup(E);
799 evalue_substitute(E, subs);
800 reduce_evalue(E);
801
802 for (j = 0; j < dim; ++j)
803 evalue_free(subs[j]);
804 free(subs);
805
806 if (new_dim-new_nparam > 0) {
807 sum = sum_over_polytope(P, E, new_dim-new_nparam, data, options);
808 evalue_free(E);
809 Polyhedron_Free(P);
810 } else {
811 sum = ALLOC(evalue);
812 value_init(sum->d);
813 sum->x.p = new_enode(partition, 2, new_dim);
814 EVALUE_SET_DOMAIN(sum->x.p->arr[0], P);
815 value_clear(sum->x.p->arr[1].d);
816 sum->x.p->arr[1] = *E;
817 free(E);
818 }
819
820 if (CP) {
821 evalue_backsubstitute(sum, CP, options->MaxRays);
822 Matrix_Free(CP);
823 }
824
825 if (T)
826 Matrix_Free(T);
827
828 return sum;
829 }
830
831 /* Splinter P into period slices along the first variable x = period y + c,
832 * 0 <= c < perdiod, * ensuring no fractional part contains the first variable,
833 * and sum over all slices.
834 */
835 static evalue *sum_over_polytope_slices(Polyhedron *P, evalue *E,
836 unsigned nvar,
837 Value period,
838 struct Bernoulli_data *data,
839 struct barvinok_options *options)
840 {
841 evalue *sum = evalue_zero();
842 evalue **subs;
843 unsigned dim = P->Dimension;
844 Matrix *T;
845 Value i;
846 Value one;
847 int j;
848
849 value_init(i);
850 value_init(one);
851 value_set_si(one, 1);
852
853 subs = ALLOCN(evalue *, dim);
854
855 T = Matrix_Alloc(dim+1, dim+1);
856 value_assign(T->p[0][0], period);
857 for (j = 1; j < dim+1; ++j)
858 value_set_si(T->p[j][j], 1);
859
860 for (j = 0; j < dim; ++j)
861 subs[j] = evalue_var(j);
862 evalue_mul(subs[0], period);
863
864 for (value_set_si(i, 0); value_lt(i, period); value_increment(i, i)) {
865 evalue *tmp;
866 Polyhedron *S = DomainPreimage(P, T, options->MaxRays);
867 evalue *e = evalue_dup(E);
868 evalue_substitute(e, subs);
869 reduce_evalue(e);
870
871 if (S->NbEq)
872 tmp = sum_over_polytope_with_equalities(S, e, nvar, data, options);
873 else
874 tmp = sum_over_polytope_base(S, e, nvar, data, options);
875
876 assert(tmp);
877 eadd(tmp, sum);
878 evalue_free(tmp);
879
880 Domain_Free(S);
881 evalue_free(e);
882
883 value_increment(T->p[0][dim], T->p[0][dim]);
884 evalue_add_constant(subs[0], one);
885 }
886
887 value_clear(i);
888 value_clear(one);
889 Matrix_Free(T);
890 for (j = 0; j < dim; ++j)
891 evalue_free(subs[j]);
892 free(subs);
893
894 reduce_evalue(sum);
895 return sum;
896 }
897
898 static evalue *sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
899 struct Bernoulli_data *data,
900 struct barvinok_options *options)
901 {
902 Polyhedron *P_orig = P;
903 evalue *E_orig = E;
904 Value period;
905 evalue *sum;
906 int exact = options->approximation_method == BV_APPROX_NONE;
907
908 if (P->NbEq)
909 return sum_over_polytope_with_equalities(P, E, nvar, data, options);
910
911 value_init(period);
912
913 move_best_to_front(&P, &E, nvar, exact ? &period : NULL);
914 if (exact && value_notone_p(period))
915 sum = sum_over_polytope_slices(P, E, nvar, period, data, options);
916 else
917 sum = sum_over_polytope_base(P, E, nvar, data, options);
918
919 if (P != P_orig)
920 Polyhedron_Free(P);
921 if (E != E_orig)
922 evalue_free(E);
923
924 value_clear(period);
925
926 return sum;
927 }
928
929 evalue *Bernoulli_sum_evalue(evalue *e, unsigned nvar,
930 struct barvinok_options *options)
931 {
932 struct Bernoulli_data data;
933 int i, j;
934 evalue *sum = evalue_zero();
935
936 if (EVALUE_IS_ZERO(*e))
937 return sum;
938
939 if (nvar == 0) {
940 eadd(e, sum);
941 return sum;
942 }
943
944 assert(value_zero_p(e->d));
945 assert(e->x.p->type == partition);
946
947 data.size = 16;
948 data.s = ALLOCN(struct evalue_section, data.size);
949 data.options = options;
950
951 for (i = 0; i < e->x.p->size/2; ++i) {
952 Polyhedron *D;
953 for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
954 Polyhedron *next = D->next;
955 evalue *tmp;
956 D->next = NULL;
957
958 tmp = sum_over_polytope(D, &e->x.p->arr[2*i+1], nvar, &data, options);
959 assert(tmp);
960 eadd(tmp, sum);
961 evalue_free(tmp);
962
963 D->next = next;
964 }
965 }
966
967 free(data.s);
968
969 reduce_evalue(sum);
970 return sum;
971 }
972
973 evalue *Bernoulli_sum(Polyhedron *P, Polyhedron *C,
974 struct barvinok_options *options)
975 {
976 Polyhedron *CA, *D;
977 evalue e;
978 evalue *sum;
979
980 if (emptyQ(P) || emptyQ(C))
981 return evalue_zero();
982
983 CA = align_context(C, P->Dimension, options->MaxRays);
984 D = DomainIntersection(P, CA, options->MaxRays);
985 Domain_Free(CA);
986
987 if (emptyQ(D)) {
988 Domain_Free(D);
989 return evalue_zero();
990 }
991
992 value_init(e.d);
993 e.x.p = new_enode(partition, 2, P->Dimension);
994 EVALUE_SET_DOMAIN(e.x.p->arr[0], D);
995 evalue_set_si(&e.x.p->arr[1], 1, 1);
996 sum = Bernoulli_sum_evalue(&e, P->Dimension - C->Dimension, options);
997 free_evalue_refs(&e);
998 return sum;
999 }
lattice point enumeration library