?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
evalue.c: reduce_evalue: reduce floor with constant argument
[barvinok/uuh.git] / bernoulli.c
blobc40c5a29efbd000b0625b853bd5eb5cdac397bc1
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 static evalue *shifted_copy(const evalue *src)
155 {
156 evalue *e = ALLOC(evalue);
157 value_init(e->d);
158 evalue_copy(e, src);
159 evalue_shift_variables(e, -1);
160 return e;
161 }
162
163 /* Computes the argument for the Faulhaber polynomial.
164 * If we are computing a polynomial approximation (exact == 0),
165 * then this is simply arg/denom.
166 * Otherwise, if neg == 0, then we are dealing with an upper bound
167 * and we need to compute floor(arg/denom) = arg/denom - { arg/denom }
168 * If neg == 1, then we are dealing with a lower bound
169 * and we need to compute ceil(arg/denom) = arg/denom + { -arg/denom }
170 *
171 * Modifies arg (if exact == 1).
172 */
173 static evalue *power_sums_base(Vector *arg, Value denom, int neg, int exact)
174 {
175 evalue *fract;
176 evalue *base = affine2evalue(arg->p, denom, arg->Size-1);
177
178 if (!exact || value_one_p(denom))
179 return base;
180
181 if (neg)
182 Vector_Oppose(arg->p, arg->p, arg->Size);
183
184 fract = fractional_part(arg->p, denom, arg->Size-1, NULL);
185 if (!neg) {
186 value_set_si(arg->p[0], -1);
187 evalue_mul(fract, arg->p[0]);
188 }
189 eadd(fract, base);
190 evalue_free(fract);
191
192 return base;
193 }
194
195 static evalue *power_sums(struct poly_list *faulhaber, const evalue *poly,
196 Vector *arg, Value denom, int neg, int alt_neg,
197 int exact)
198 {
199 int i;
200 evalue *base = power_sums_base(arg, denom, neg, exact);
201 evalue *sum = evalue_zero();
202
203 for (i = 1; i < poly->x.p->size; ++i) {
204 evalue *term = evalue_polynomial(faulhaber->poly[i], base);
205 evalue *factor = shifted_copy(&poly->x.p->arr[i]);
206 emul(factor, term);
207 if (alt_neg && (i % 2))
208 evalue_negate(term);
209 eadd(term, sum);
210 evalue_free(factor);
211 evalue_free(term);
212 }
213 if (neg)
214 evalue_negate(sum);
215 evalue_free(base);
216
217 return sum;
218 }
219
220 /* Given a constraint (cst_affine) a x + b y + c >= 0, compate a constraint (c)
221 * +- (b y + c) +- a >=,> 0
222 * ^ ^ ^
223 * | | strict
224 * sign_affine sign_cst
225 */
226 static void bound_constraint(Value *c, unsigned dim,
227 Value *cst_affine,
228 int sign_affine, int sign_cst, int strict)
229 {
230 if (sign_affine == -1)
231 Vector_Oppose(cst_affine+1, c, dim+1);
232 else
233 Vector_Copy(cst_affine+1, c, dim+1);
234
235 if (sign_cst == -1)
236 value_subtract(c[dim], c[dim], cst_affine[0]);
237 else if (sign_cst == 1)
238 value_addto(c[dim], c[dim], cst_affine[0]);
239
240 if (strict)
241 value_decrement(c[dim], c[dim]);
242 }
243
244 struct Bernoulli_data {
245 struct barvinok_options *options;
246 struct evalue_section *s;
247 int size;
248 int ns;
249 const evalue *e;
250 };
251
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->approximation_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
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->approximation_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 */
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
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->approximation_method == BV_APPROX_NONE;
325 int cases = exact ? 3 : 5;
326
327 if (cases * n <= data->size)
328 return;
329
330 data->size = cases * (n + 16);
331 data->s = REALLOCN(data->s, struct evalue_section, data->size);
332 }
333
334 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data);
335
336 /*
337 * This function requires that either the lower or the upper bound
338 * represented by the constraints "lower" and "upper" is an integer
339 * affine expression.
340 * An affine substitution is performed to make this bound exactly
341 * zero, ensuring that in the recursive call to Bernoulli_cb,
342 * only one of the "cases" will apply.
343 */
344 static void transform_to_single_case(Matrix *M, Value *lower, Value *upper,
345 struct Bernoulli_data *data)
346 {
347 unsigned dim = M->NbColumns-2;
348 Vector *shadow;
349 Value a, b;
350 evalue **subs;
351 const evalue *e = data->e;
352 evalue *t;
353 int i;
354
355 value_init(a);
356 value_init(b);
357 subs = ALLOCN(evalue *, dim+1);
358 for (i = 0; i < dim; ++i)
359 subs[1+i] = evalue_var(1+i);
360 shadow = Vector_Alloc(2 * (2+dim+1));
361 if (value_one_p(lower[1])) {
362 /* Replace i by i + l' when b = 1 */
363 value_set_si(shadow->p[0], 1);
364 Vector_Oppose(lower+2, shadow->p+1, dim+1);
365 subs[0] = affine2evalue(shadow->p, shadow->p[0], dim+1);
366 /* new lower
367 * i >= 0
368 * new upper
369 * (-a i + u') + a (-l') >= 0
370 */
371 value_assign(shadow->p[2+dim+1+1], upper[1]);
372 value_oppose(a, upper[1]);
373 value_set_si(b, 1);
374 Vector_Combine(upper+2, lower+2, shadow->p+2+dim+1+2,
375 b, a, dim+1);
376 upper = shadow->p+2+dim+1;
377 lower = shadow->p;
378 value_set_si(lower[1], 1);
379 Vector_Set(lower+2, 0, dim+1);
380 } else {
381 /* Replace i by i + u' when a = 1 */
382 value_set_si(shadow->p[0], 1);
383 Vector_Copy(upper+2, shadow->p+1, dim+1);
384 subs[0] = affine2evalue(shadow->p, shadow->p[0], dim+1);
385 /* new lower
386 * (b i - l') + b u' >= 0
387 * new upper
388 * -i >= 0
389 */
390 value_assign(shadow->p[1], lower[1]);
391 value_set_si(a, 1);
392 value_assign(b, lower[1]);
393 Vector_Combine(upper+2, lower+2, shadow->p+2,
394 b, a, dim+1);
395 upper = shadow->p+2+dim+1;
396 lower = shadow->p;
397 value_set_si(upper[1], -1);
398 Vector_Set(upper+2, 0, dim+1);
399 }
400 value_clear(a);
401 value_clear(b);
402
403 t = evalue_dup(data->e);
404 evalue_substitute(t, subs);
405 reduce_evalue(t);
406 data->e = t;
407 for (i = 0; i < dim+1; ++i)
408 evalue_free(subs[i]);
409 free(subs);
410
411 Bernoulli_cb(M, lower, upper, data);
412
413 evalue_free(t);
414 data->e = e;
415 Vector_Free(shadow);
416 }
417
418 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
419 {
420 struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
421 Matrix *M2;
422 Polyhedron *T;
423 const evalue *factor = NULL;
424 evalue *linear = NULL;
425 int constant = 0;
426 Value tmp;
427 unsigned dim = M->NbColumns-2;
428 Vector *row;
429 int exact = data->options->approximation_method == BV_APPROX_NONE;
430 int cases = exact ? 3 : 5;
431
432 assert(lower);
433 assert(upper);
434 assert(data->ns + cases <= data->size);
435
436 M2 = Matrix_Copy(M);
437 T = Constraints2Polyhedron(M2, data->options->MaxRays);
438 Matrix_Free(M2);
439
440 POL_ENSURE_VERTICES(T);
441 if (emptyQ(T)) {
442 Polyhedron_Free(T);
443 return;
444 }
445
446 constant = value_notzero_p(data->e->d) ||
447 data->e->x.p->type != polynomial ||
448 data->e->x.p->pos != 1;
449 if (!constant && (value_one_p(lower[1]) || value_mone_p(upper[1]))) {
450 int single_case;
451 int lower_cst, upper_cst;
452
453 lower_cst = First_Non_Zero(lower+2, dim) == -1;
454 upper_cst = First_Non_Zero(upper+2, dim) == -1;
455 single_case =
456 (lower_cst && value_negz_p(lower[2+dim])) ||
457 (upper_cst && value_negz_p(upper[2+dim])) ||
458 (lower_cst && upper_cst &&
459 value_posz_p(lower[2+dim]) && value_posz_p(upper[2+dim]));
460 if (!single_case) {
461 transform_to_single_case(M, lower, upper, data);
462 Polyhedron_Free(T);
463 return;
464 }
465 }
466
467 assert(lower != upper);
468
469 row = Vector_Alloc(dim+1);
470 value_init(tmp);
471 if (value_notzero_p(data->e->d)) {
472 factor = data->e;
473 constant = 1;
474 } else {
475 if (data->e->x.p->type == polynomial && data->e->x.p->pos == 1)
476 factor = shifted_copy(&data->e->x.p->arr[0]);
477 else {
478 factor = shifted_copy(data->e);
479 constant = 1;
480 }
481 }
482 linear = linear_term(factor, lower, upper, row, tmp, exact);
483
484 if (constant) {
485 data->s[data->ns].E = linear;
486 data->s[data->ns].D = T;
487 ++data->ns;
488 data->options->stats->bernoulli_sums++;
489 } else {
490 evalue *poly_u = NULL, *poly_l = NULL;
491 Polyhedron *D;
492 struct poly_list *faulhaber;
493 assert(data->e->x.p->type == polynomial);
494 assert(data->e->x.p->pos == 1);
495 faulhaber = faulhaber_compute(data->e->x.p->size-1);
496 /* lower is the constraint
497 * b i - l' >= 0 i >= l'/b = l
498 * upper is the constraint
499 * -a i + u' >= 0 i <= u'/a = u
500 */
501 M2 = Matrix_Alloc(3, 2+T->Dimension);
502 value_set_si(M2->p[0][0], 1);
503 value_set_si(M2->p[1][0], 1);
504 value_set_si(M2->p[2][0], 1);
505 /* Case 1:
506 * 1 <= l l' - b >= 0
507 * 0 < l l' - 1 >= 0 if exact
508 */
509 if (exact)
510 bound_constraint(M2->p[0]+1, T->Dimension, lower+1, -1, 0, 1);
511 else
512 bound_constraint(M2->p[0]+1, T->Dimension, lower+1, -1, -1, 0);
513 D = AddConstraints(M2->p_Init, 1, T, data->options->MaxRays);
514 POL_ENSURE_VERTICES(D);
515 if (emptyQ2(D))
516 Polyhedron_Free(D);
517 else {
518 evalue *extra;
519 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
520 faulhaber, data);
521 Vector_Oppose(lower+2, row->p, dim+1);
522 extra = power_sums(faulhaber, data->e, row, lower[1], 1, 0, exact);
523 eadd(poly_u, extra);
524 eadd(linear, extra);
525
526 data->s[data->ns].E = extra;
527 data->s[data->ns].D = D;
528 ++data->ns;
529 data->options->stats->bernoulli_sums++;
530 }
531
532 /* Case 2:
533 * 1 <= -u -u' - a >= 0
534 * 0 < -u -u' - 1 >= 0 if exact
535 */
536 if (exact)
537 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, -1, 0, 1);
538 else
539 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, -1, 1, 0);
540 D = AddConstraints(M2->p_Init, 1, T, data->options->MaxRays);
541 POL_ENSURE_VERTICES(D);
542 if (emptyQ2(D))
543 Polyhedron_Free(D);
544 else {
545 evalue *extra;
546 poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
547 Vector_Oppose(upper+2, row->p, dim+1);
548 value_oppose(tmp, upper[1]);
549 extra = power_sums(faulhaber, data->e, row, tmp, 1, 1, exact);
550 eadd(poly_l, extra);
551 eadd(linear, extra);
552
553 data->s[data->ns].E = extra;
554 data->s[data->ns].D = D;
555 ++data->ns;
556 data->options->stats->bernoulli_sums++;
557 }
558
559 /* Case 3:
560 * u >= 0 u' >= 0
561 * -l >= 0 -l' >= 0
562 */
563 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, 0, 0);
564 bound_constraint(M2->p[1]+1, T->Dimension, lower+1, 1, 0, 0);
565 D = AddConstraints(M2->p_Init, 2, T, data->options->MaxRays);
566 POL_ENSURE_VERTICES(D);
567 if (emptyQ2(D))
568 Polyhedron_Free(D);
569 else {
570 poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
571 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
572 faulhaber, data);
573
574 data->s[data->ns].E = ALLOC(evalue);
575 value_init(data->s[data->ns].E->d);
576 evalue_copy(data->s[data->ns].E, poly_u);
577 eadd(poly_l, data->s[data->ns].E);
578 eadd(linear, data->s[data->ns].E);
579 data->s[data->ns].D = D;
580 ++data->ns;
581 data->options->stats->bernoulli_sums++;
582 }
583
584 if (!exact) {
585 /* Case 4:
586 * l < 1 -l' + b - 1 >= 0
587 * 0 < l l' - 1 >= 0
588 * u >= 1 u' - a >= 0
589 */
590 bound_constraint(M2->p[0]+1, T->Dimension, lower+1, 1, 1, 1);
591 bound_constraint(M2->p[1]+1, T->Dimension, lower+1, -1, 0, 1);
592 bound_constraint(M2->p[2]+1, T->Dimension, upper+1, 1, 1, 0);
593 D = AddConstraints(M2->p_Init, 3, T, data->options->MaxRays);
594 POL_ENSURE_VERTICES(D);
595 if (emptyQ2(D))
596 Polyhedron_Free(D);
597 else {
598 poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
599 faulhaber, data);
600 eadd(linear, poly_u);
601 data->s[data->ns].E = poly_u;
602 poly_u = NULL;
603 data->s[data->ns].D = D;
604 ++data->ns;
605 data->options->stats->bernoulli_sums++;
606 }
607
608 /* Case 5:
609 * -u < 1 u' + a - 1 >= 0
610 * 0 < -u -u' - 1 >= 0
611 * l <= -1 -l' - b >= 0
612 */
613 bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, -1, 1);
614 bound_constraint(M2->p[1]+1, T->Dimension, upper+1, -1, 0, 1);
615 bound_constraint(M2->p[2]+1, T->Dimension, lower+1, 1, -1, 0);
616 D = AddConstraints(M2->p_Init, 3, T, data->options->MaxRays);
617 POL_ENSURE_VERTICES(D);
618 if (emptyQ2(D))
619 Polyhedron_Free(D);
620 else {
621 poly_l = compute_poly_l(poly_l, lower, row, dim,
622 faulhaber, data);
623 eadd(linear, poly_l);
624 data->s[data->ns].E = poly_l;
625 poly_l = NULL;
626 data->s[data->ns].D = D;
627 ++data->ns;
628 data->options->stats->bernoulli_sums++;
629 }
630 }
631
632 Matrix_Free(M2);
633 Polyhedron_Free(T);
634 if (poly_l)
635 evalue_free(poly_l);
636 if (poly_u)
637 evalue_free(poly_u);
638 evalue_free(linear);
639 }
640 if (factor != data->e)
641 evalue_free((evalue *)factor);
642 value_clear(tmp);
643 Vector_Free(row);
644 }
645
646 /*
647 * Move the variable at position pos to the front by exchanging
648 * that variable with the first variable.
649 */
650 static void more_var_to_front(Polyhedron **P_p , evalue **E_p, int pos)
651 {
652 Polyhedron *P = *P_p;
653 evalue *E = *E_p;
654 unsigned dim = P->Dimension;
655
656 assert(pos != 0);
657
658 P = Polyhedron_Copy(P);
659 Polyhedron_ExchangeColumns(P, 1, 1+pos);
660 *P_p = P;
661
662 if (value_zero_p(E->d)) {
663 int j;
664 evalue **subs;
665
666 subs = ALLOCN(evalue *, dim);
667 for (j = 0; j < dim; ++j)
668 subs[j] = evalue_var(j);
669 E = subs[0];
670 subs[0] = subs[pos];
671 subs[pos] = E;
672 E = evalue_dup(*E_p);
673 evalue_substitute(E, subs);
674 for (j = 0; j < dim; ++j)
675 evalue_free(subs[j]);
676 free(subs);
677 *E_p = E;
678 }
679 }
680
681 static int type_offset(enode *p)
682 {
683 return p->type == fractional ? 1 :
684 p->type == flooring ? 1 : 0;
685 }
686
687 static void adjust_periods(evalue *e, unsigned nvar, Vector *periods)
688 {
689 for (; value_zero_p(e->d); e = &e->x.p->arr[0]) {
690 int pos;
691 assert(e->x.p->type == polynomial);
692 assert(e->x.p->size == 2);
693 assert(value_notzero_p(e->x.p->arr[1].d));
694
695 pos = e->x.p->pos - 1;
696 if (pos >= nvar)
697 break;
698
699 value_lcm(periods->p[pos], periods->p[pos], e->x.p->arr[1].d);
700 }
701 }
702
703 static void fractional_periods_r(evalue *e, unsigned nvar, Vector *periods)
704 {
705 int i;
706
707 if (value_notzero_p(e->d))
708 return;
709
710 assert(e->x.p->type == polynomial || e->x.p->type == fractional);
711
712 if (e->x.p->type == fractional)
713 adjust_periods(&e->x.p->arr[0], nvar, periods);
714
715 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
716 fractional_periods_r(&e->x.p->arr[i], nvar, periods);
717 }
718
719 /*
720 * For each of the nvar variables, compute the lcm of the
721 * denominators of the coefficients of that variable in
722 * any of the fractional parts.
723 */
724 static Vector *fractional_periods(evalue *e, unsigned nvar)
725 {
726 int i;
727 Vector *periods = Vector_Alloc(nvar);
728
729 for (i = 0; i < nvar; ++i)
730 value_set_si(periods->p[i], 1);
731
732 fractional_periods_r(e, nvar, periods);
733
734 return periods;
735 }
736
737 /* Move "best" variable to sum over into the first position,
738 * possibly changing *P_p and *E_p.
739 *
740 * If there are any fractional parts (period != NULL), then we
741 * first look for a variable with the smallest lcm of denominators
742 * inside a fractional part. This denominator is assigned to period
743 * and corresponds to the number of "splinters" we need to construct
744 * at this level.
745 *
746 * Of those with this denominator (all if period == NULL or if there
747 * are no fractional parts), we select the variable with the smallest
748 * maximal coefficient, as this coefficient will become a denominator
749 * in new fractional parts (for an exact computation), which may
750 * lead to splintering in the next step.
751 */
752 static void move_best_to_front(Polyhedron **P_p, evalue **E_p, unsigned nvar,
753 Value *period)
754 {
755 Polyhedron *P = *P_p;
756 evalue *E = *E_p;
757 int i, j, best_i = -1;
758 Vector *periods;
759 Value best, max;
760
761 assert(nvar >= 1);
762
763 if (period) {
764 periods = fractional_periods(E, nvar);
765 value_assign(*period, periods->p[0]);
766 for (i = 1; i < nvar; ++i)
767 if (value_lt(periods->p[i], *period))
768 value_assign(*period, periods->p[i]);
769 }
770
771 value_init(best);
772 value_init(max);
773
774 for (i = 0; i < nvar; ++i) {
775 if (period && value_ne(*period, periods->p[i]))
776 continue;
777
778 value_set_si(max, 0);
779
780 for (j = 0; j < P->NbConstraints; ++j) {
781 if (value_zero_p(P->Constraint[j][1+i]))
782 continue;
783 if (First_Non_Zero(P->Constraint[j]+1, i) == -1 &&
784 First_Non_Zero(P->Constraint[j]+1+i+1, nvar-i-1) == -1)
785 continue;
786 if (value_abs_gt(P->Constraint[j][1+i], max))
787 value_absolute(max, P->Constraint[j][1+i]);
788 }
789
790 if (best_i == -1 || value_lt(max, best)) {
791 value_assign(best, max);
792 best_i = i;
793 }
794 }
795
796 value_clear(best);
797 value_clear(max);
798
799 if (period)
800 Vector_Free(periods);
801
802 if (best_i != 0)
803 more_var_to_front(P_p, E_p, best_i);
804
805 return;
806 }
807
808 static evalue *sum_over_polytope_base(Polyhedron *P, evalue *E, unsigned nvar,
809 struct Bernoulli_data *data,
810 struct barvinok_options *options)
811 {
812 evalue *res;
813
814 assert(P->NbEq == 0);
815
816 data->ns = 0;
817 data->e = E;
818
819 for_each_lower_upper_bound(P, Bernoulli_init, Bernoulli_cb, data);
820
821 res = evalue_from_section_array(data->s, data->ns);
822
823 if (nvar > 1) {
824 evalue *tmp = Bernoulli_sum_evalue(res, nvar-1, options);
825 evalue_free(res);
826 res = tmp;
827 }
828
829 return res;
830 }
831
832 static evalue *sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
833 struct Bernoulli_data *data,
834 struct barvinok_options *options);
835
836 static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E,
837 unsigned nvar,
838 struct Bernoulli_data *data,
839 struct barvinok_options *options)
840 {
841 unsigned dim = P->Dimension;
842 unsigned new_dim, new_nparam;
843 Matrix *T = NULL, *CP = NULL;
844 evalue **subs;
845 evalue *sum;
846 int j;
847
848 if (emptyQ(P))
849 return evalue_zero();
850
851 assert(P->NbEq > 0);
852
853 remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays);
854
855 if (emptyQ(P)) {
856 Polyhedron_Free(P);
857 return evalue_zero();
858 }
859
860 new_nparam = CP ? CP->NbColumns-1 : dim - nvar;
861 new_dim = T ? T->NbColumns-1 : nvar + new_nparam;
862
863 /* We can avoid these substitutions if E is a constant */
864 subs = ALLOCN(evalue *, dim);
865 for (j = 0; j < nvar; ++j) {
866 if (T)
867 subs[j] = affine2evalue(T->p[j], T->p[nvar+new_nparam][new_dim],
868 new_dim);
869 else
870 subs[j] = evalue_var(j);
871 }
872 for (j = 0; j < dim-nvar; ++j) {
873 if (CP)
874 subs[nvar+j] = affine2evalue(CP->p[j], CP->p[dim-nvar][new_nparam],
875 new_nparam);
876 else
877 subs[nvar+j] = evalue_var(j);
878 evalue_shift_variables(subs[nvar+j], new_dim-new_nparam);
879 }
880
881 E = evalue_dup(E);
882 evalue_substitute(E, subs);
883 reduce_evalue(E);
884
885 for (j = 0; j < dim; ++j)
886 evalue_free(subs[j]);
887 free(subs);
888
889 if (new_dim-new_nparam > 0) {
890 sum = sum_over_polytope(P, E, new_dim-new_nparam, data, options);
891 evalue_free(E);
892 Polyhedron_Free(P);
893 } else {
894 sum = ALLOC(evalue);
895 value_init(sum->d);
896 sum->x.p = new_enode(partition, 2, new_dim);
897 EVALUE_SET_DOMAIN(sum->x.p->arr[0], P);
898 value_clear(sum->x.p->arr[1].d);
899 sum->x.p->arr[1] = *E;
900 free(E);
901 }
902
903 if (CP) {
904 evalue_backsubstitute(sum, CP, options->MaxRays);
905 Matrix_Free(CP);
906 }
907
908 if (T)
909 Matrix_Free(T);
910
911 return sum;
912 }
913
914 /* Splinter P into period slices along the first variable x = period y + c,
915 * 0 <= c < perdiod, * ensuring no fractional part contains the first variable,
916 * and sum over all slices.
917 */
918 static evalue *sum_over_polytope_slices(Polyhedron *P, evalue *E,
919 unsigned nvar,
920 Value period,
921 struct Bernoulli_data *data,
922 struct barvinok_options *options)
923 {
924 evalue *sum = evalue_zero();
925 evalue **subs;
926 unsigned dim = P->Dimension;
927 Matrix *T;
928 Value i;
929 Value one;
930 int j;
931
932 value_init(i);
933 value_init(one);
934 value_set_si(one, 1);
935
936 subs = ALLOCN(evalue *, dim);
937
938 T = Matrix_Alloc(dim+1, dim+1);
939 value_assign(T->p[0][0], period);
940 for (j = 1; j < dim+1; ++j)
941 value_set_si(T->p[j][j], 1);
942
943 for (j = 0; j < dim; ++j)
944 subs[j] = evalue_var(j);
945 evalue_mul(subs[0], period);
946
947 for (value_set_si(i, 0); value_lt(i, period); value_increment(i, i)) {
948 evalue *tmp;
949 Polyhedron *S = DomainPreimage(P, T, options->MaxRays);
950 evalue *e = evalue_dup(E);
951 evalue_substitute(e, subs);
952 reduce_evalue(e);
953
954 if (S->NbEq)
955 tmp = sum_over_polytope_with_equalities(S, e, nvar, data, options);
956 else
957 tmp = sum_over_polytope_base(S, e, nvar, data, options);
958
959 assert(tmp);
960 eadd(tmp, sum);
961 evalue_free(tmp);
962
963 Domain_Free(S);
964 evalue_free(e);
965
966 value_increment(T->p[0][dim], T->p[0][dim]);
967 evalue_add_constant(subs[0], one);
968 }
969
970 value_clear(i);
971 value_clear(one);
972 Matrix_Free(T);
973 for (j = 0; j < dim; ++j)
974 evalue_free(subs[j]);
975 free(subs);
976
977 reduce_evalue(sum);
978 return sum;
979 }
980
981 static evalue *sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
982 struct Bernoulli_data *data,
983 struct barvinok_options *options)
984 {
985 Polyhedron *P_orig = P;
986 evalue *E_orig = E;
987 Value period;
988 evalue *sum;
989 int exact = options->approximation_method == BV_APPROX_NONE;
990
991 if (P->NbEq)
992 return sum_over_polytope_with_equalities(P, E, nvar, data, options);
993
994 value_init(period);
995
996 move_best_to_front(&P, &E, nvar, exact ? &period : NULL);
997 if (exact && value_notone_p(period))
998 sum = sum_over_polytope_slices(P, E, nvar, period, data, options);
999 else
1000 sum = sum_over_polytope_base(P, E, nvar, data, options);
1001
1002 if (P != P_orig)
1003 Polyhedron_Free(P);
1004 if (E != E_orig)
1005 evalue_free(E);
1006
1007 value_clear(period);
1008
1009 return sum;
1010 }
1011
1012 evalue *Bernoulli_sum_evalue(evalue *e, unsigned nvar,
1013 struct barvinok_options *options)
1014 {
1015 struct Bernoulli_data data;
1016 int i, j;
1017 evalue *sum = evalue_zero();
1018
1019 if (EVALUE_IS_ZERO(*e))
1020 return sum;
1021
1022 if (nvar == 0) {
1023 eadd(e, sum);
1024 return sum;
1025 }
1026
1027 assert(value_zero_p(e->d));
1028 assert(e->x.p->type == partition);
1029
1030 data.size = 16;
1031 data.s = ALLOCN(struct evalue_section, data.size);
1032 data.options = options;
1033
1034 for (i = 0; i < e->x.p->size/2; ++i) {
1035 Polyhedron *D;
1036 for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
1037 Polyhedron *next = D->next;
1038 evalue *tmp;
1039 D->next = NULL;
1040
1041 tmp = sum_over_polytope(D, &e->x.p->arr[2*i+1], nvar, &data, options);
1042 assert(tmp);
1043 eadd(tmp, sum);
1044 evalue_free(tmp);
1045
1046 D->next = next;
1047 }
1048 }
1049
1050 free(data.s);
1051
1052 reduce_evalue(sum);
1053 return sum;
1054 }
patches by uuh to consider upstream