evalue_isl.c: isl_pw_qpolynomial_from_eval: use isl_val
[barvinok.git] / summate.c
blobd0608a5689a87814d05374cc01b48e0aec14ed8a
1 #include <isl/map.h>
2 #include <isl_set_polylib.h>
3 #include <barvinok/options.h>
4 #include <barvinok/util.h>
5 #include "bernoulli.h"
6 #include "euler.h"
7 #include "laurent.h"
8 #include "laurent_old.h"
9 #include "summate.h"
10 #include "section_array.h"
11 #include "remove_equalities.h"
13 extern evalue *evalue_outer_floor(evalue *e);
14 extern int evalue_replace_floor(evalue *e, const evalue *floor, int var);
15 extern void evalue_drop_floor(evalue *e, const evalue *floor);
17 #define ALLOC(type) (type*)malloc(sizeof(type))
18 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
20 /* Apply the variable transformation specified by T and CP on
21 * the polynomial e. T expresses the old variables in terms
22 * of the new variables (and optionally also the new parameters),
23 * while CP expresses the old parameters in terms of the new
24 * parameters.
26 static void transform_polynomial(evalue *E, Matrix *T, Matrix *CP,
27 unsigned nvar, unsigned nparam,
28 unsigned new_nvar, unsigned new_nparam)
30 int j;
31 evalue **subs;
33 subs = ALLOCN(evalue *, nvar+nparam);
35 for (j = 0; j < nvar; ++j) {
36 if (T)
37 subs[j] = affine2evalue(T->p[j], T->p[T->NbRows-1][T->NbColumns-1],
38 T->NbColumns-1);
39 else
40 subs[j] = evalue_var(j);
42 for (j = 0; j < nparam; ++j) {
43 if (CP)
44 subs[nvar+j] = affine2evalue(CP->p[j], CP->p[nparam][new_nparam],
45 new_nparam);
46 else
47 subs[nvar+j] = evalue_var(j);
48 evalue_shift_variables(subs[nvar+j], 0, new_nvar);
51 evalue_substitute(E, subs);
52 reduce_evalue(E);
54 for (j = 0; j < nvar+nparam; ++j)
55 evalue_free(subs[j]);
56 free(subs);
59 /* Compute the sum of the quasi-polynomial E
60 * over a 0D (non-empty, but possibly parametric) polytope P.
62 * Consumes P and E.
64 * We simply return a partition evalue with P as domain and E as value.
66 static evalue *sum_over_polytope_0D(Polyhedron *P, evalue *E)
68 evalue *sum;
70 sum = ALLOC(evalue);
71 value_init(sum->d);
72 sum->x.p = new_enode(partition, 2, P->Dimension);
73 EVALUE_SET_DOMAIN(sum->x.p->arr[0], P);
74 value_clear(sum->x.p->arr[1].d);
75 sum->x.p->arr[1] = *E;
76 free(E);
78 return sum;
81 static evalue *sum_with_equalities(Polyhedron *P, evalue *E,
82 unsigned nvar, struct evalue_section_array *sections,
83 struct barvinok_options *options,
84 evalue *(*base)(Polyhedron *P, evalue *E, unsigned nvar,
85 struct evalue_section_array *sections,
86 struct barvinok_options *options))
88 unsigned dim = P->Dimension;
89 unsigned new_dim, new_nparam;
90 Matrix *T = NULL, *CP = NULL;
91 evalue *sum;
93 if (emptyQ(P))
94 return evalue_zero();
96 assert(P->NbEq > 0);
98 remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays);
100 if (emptyQ(P)) {
101 Polyhedron_Free(P);
102 return evalue_zero();
105 new_nparam = CP ? CP->NbColumns-1 : dim - nvar;
106 new_dim = T ? T->NbColumns-1 : nvar + new_nparam;
108 /* We can avoid these substitutions if E is a constant */
109 E = evalue_dup(E);
110 transform_polynomial(E, T, CP, nvar, dim-nvar,
111 new_dim-new_nparam, new_nparam);
113 if (new_dim-new_nparam > 0) {
114 sum = base(P, E, new_dim-new_nparam, sections, options);
115 evalue_free(E);
116 Polyhedron_Free(P);
117 } else {
118 sum = sum_over_polytope_0D(P, E);
121 if (CP) {
122 evalue_backsubstitute(sum, CP, options->MaxRays);
123 Matrix_Free(CP);
126 if (T)
127 Matrix_Free(T);
129 return sum;
132 static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E,
133 unsigned nvar, struct evalue_section_array *sections,
134 struct barvinok_options *options)
136 return sum_with_equalities(P, E, nvar, sections, options,
137 &barvinok_sum_over_polytope);
140 static evalue *sum_base(Polyhedron *P, evalue *E, unsigned nvar,
141 struct barvinok_options *options);
143 static evalue *sum_base_wrap(Polyhedron *P, evalue *E, unsigned nvar,
144 struct evalue_section_array *sections, struct barvinok_options *options)
146 return sum_base(P, E, nvar, options);
149 static evalue *sum_base_with_equalities(Polyhedron *P, evalue *E,
150 unsigned nvar, struct barvinok_options *options)
152 return sum_with_equalities(P, E, nvar, NULL, options, &sum_base_wrap);
155 /* The substitutions in sum_step_polynomial may have reintroduced equalities
156 * (in particular, one of the floor expressions may be equal to one of
157 * the variables), so we need to check for them again.
159 static evalue *sum_base(Polyhedron *P, evalue *E, unsigned nvar,
160 struct barvinok_options *options)
162 if (P->NbEq)
163 return sum_base_with_equalities(P, E, nvar, options);
164 if (options->summation == BV_SUM_EULER)
165 return euler_summate(P, E, nvar, options);
166 else if (options->summation == BV_SUM_LAURENT)
167 return laurent_summate(P, E, nvar, options);
168 else if (options->summation == BV_SUM_LAURENT_OLD)
169 return laurent_summate_old(P, E, nvar, options);
170 assert(0);
173 /* Count the number of non-zero terms in e when viewed as a polynomial
174 * in only the first nvar variables. "count" is the number counted
175 * so far.
177 static int evalue_count_terms(const evalue *e, unsigned nvar, int count)
179 int i;
181 if (EVALUE_IS_ZERO(*e))
182 return count;
184 if (value_zero_p(e->d))
185 assert(e->x.p->type == polynomial);
186 if (value_notzero_p(e->d) || e->x.p->pos >= nvar+1)
187 return count+1;
189 for (i = 0; i < e->x.p->size; ++i)
190 count = evalue_count_terms(&e->x.p->arr[i], nvar, count);
192 return count;
195 /* Create placeholder structure for unzipping.
196 * A "polynomial" is created with size terms in variable pos,
197 * with each term having itself as coefficient.
199 static evalue *create_placeholder(int size, int pos)
201 int i, j;
202 evalue *E = ALLOC(evalue);
203 value_init(E->d);
204 E->x.p = new_enode(polynomial, size, pos+1);
205 for (i = 0; i < size; ++i) {
206 E->x.p->arr[i].x.p = new_enode(polynomial, i+1, pos+1);
207 for (j = 0; j < i; ++j)
208 evalue_set_si(&E->x.p->arr[i].x.p->arr[j], 0, 1);
209 evalue_set_si(&E->x.p->arr[i].x.p->arr[i], 1, 1);
211 return E;
214 /* Interchange each non-zero term in e (when viewed as a polynomial
215 * in only the first nvar variables) with a placeholder in ph (created
216 * by create_placeholder), resulting in two polynomials in the
217 * placeholder variable such that for each non-zero term in e
218 * there is a power of the placeholder variable such that the factors
219 * in the first nvar variables form the coefficient of that power in
220 * the first polynomial (e) and the factors in the remaining variables
221 * form the coefficient of that power in the second polynomial (ph).
223 static int evalue_unzip_terms(evalue *e, evalue *ph, unsigned nvar, int count)
225 int i;
227 if (EVALUE_IS_ZERO(*e))
228 return count;
230 if (value_zero_p(e->d))
231 assert(e->x.p->type == polynomial);
232 if (value_notzero_p(e->d) || e->x.p->pos >= nvar+1) {
233 evalue t = *e;
234 *e = ph->x.p->arr[count];
235 ph->x.p->arr[count] = t;
236 return count+1;
239 for (i = 0; i < e->x.p->size; ++i)
240 count = evalue_unzip_terms(&e->x.p->arr[i], ph, nvar, count);
242 return count;
245 /* Remove n variables at pos (0-based) from the polyhedron P.
246 * Each of these variables is assumed to be completely free,
247 * i.e., there is a line in the polyhedron corresponding to
248 * each of these variables.
250 static Polyhedron *Polyhedron_Remove_Columns(Polyhedron *P, unsigned pos,
251 unsigned n)
253 int i, j;
254 unsigned NbConstraints = 0;
255 unsigned NbRays = 0;
256 Polyhedron *Q;
258 if (n == 0)
259 return P;
261 assert(pos <= P->Dimension);
263 if (POL_HAS(P, POL_INEQUALITIES))
264 NbConstraints = P->NbConstraints;
265 if (POL_HAS(P, POL_POINTS))
266 NbRays = P->NbRays - n;
268 Q = Polyhedron_Alloc(P->Dimension - n, NbConstraints, NbRays);
269 if (POL_HAS(P, POL_INEQUALITIES)) {
270 Q->NbEq = P->NbEq;
271 for (i = 0; i < P->NbConstraints; ++i) {
272 Vector_Copy(P->Constraint[i], Q->Constraint[i], 1+pos);
273 Vector_Copy(P->Constraint[i]+1+pos+n, Q->Constraint[i]+1+pos,
274 Q->Dimension-pos+1);
277 if (POL_HAS(P, POL_POINTS)) {
278 Q->NbBid = P->NbBid - n;
279 for (i = 0; i < n; ++i)
280 value_set_si(Q->Ray[i][1+pos+i], 1);
281 for (i = 0, j = 0; i < P->NbRays; ++i) {
282 int line = First_Non_Zero(P->Ray[i], 1+P->Dimension+1);
283 assert(line != -1);
284 if (line-1 >= pos && line-1 < pos+n) {
285 ++j;
286 assert(j <= n);
287 continue;
289 assert(i-j < Q->NbRays);
290 Vector_Copy(P->Ray[i], Q->Ray[i-j], 1+pos);
291 Vector_Copy(P->Ray[i]+1+pos+n, Q->Ray[i-j]+1+pos,
292 Q->Dimension-pos+1);
295 POL_SET(Q, POL_VALID);
296 if (POL_HAS(P, POL_INEQUALITIES))
297 POL_SET(Q, POL_INEQUALITIES);
298 if (POL_HAS(P, POL_POINTS))
299 POL_SET(Q, POL_POINTS);
300 if (POL_HAS(P, POL_VERTICES))
301 POL_SET(Q, POL_VERTICES);
302 return Q;
305 /* Remove n variables at pos (0-based) from the union of polyhedra P.
306 * Each of these variables is assumed to be completely free,
307 * i.e., there is a line in the polyhedron corresponding to
308 * each of these variables.
310 static Polyhedron *Domain_Remove_Columns(Polyhedron *P, unsigned pos,
311 unsigned n)
313 Polyhedron *R;
314 Polyhedron **next = &R;
316 for (; P; P = P->next) {
317 *next = Polyhedron_Remove_Columns(P, pos, n);
318 next = &(*next)->next;
320 return R;
323 /* Drop n parameters starting at first from partition evalue e */
324 static void drop_parameters(evalue *e, int first, int n)
326 int i;
328 if (EVALUE_IS_ZERO(*e))
329 return;
331 assert(value_zero_p(e->d) && e->x.p->type == partition);
332 for (i = 0; i < e->x.p->size/2; ++i) {
333 Polyhedron *P = EVALUE_DOMAIN(e->x.p->arr[2*i]);
334 Polyhedron *Q = Domain_Remove_Columns(P, first, n);
335 EVALUE_SET_DOMAIN(e->x.p->arr[2*i], Q);
336 Domain_Free(P);
337 evalue_shift_variables(&e->x.p->arr[2*i+1], first, -n);
339 e->x.p->pos -= n;
342 static void extract_term_into(const evalue *src, int var, int exp, evalue *dst)
344 int i;
346 if (value_notzero_p(src->d) ||
347 src->x.p->type != polynomial ||
348 src->x.p->pos > var+1) {
349 if (exp == 0)
350 evalue_copy(dst, src);
351 else
352 evalue_set_si(dst, 0, 1);
353 return;
356 if (src->x.p->pos == var+1) {
357 if (src->x.p->size > exp)
358 evalue_copy(dst, &src->x.p->arr[exp]);
359 else
360 evalue_set_si(dst, 0, 1);
361 return;
364 dst->x.p = new_enode(polynomial, src->x.p->size, src->x.p->pos);
365 for (i = 0; i < src->x.p->size; ++i)
366 extract_term_into(&src->x.p->arr[i], var, exp,
367 &dst->x.p->arr[i]);
370 /* Extract the coefficient of var^exp.
372 static evalue *extract_term(const evalue *e, int var, int exp)
374 int i;
375 evalue *res;
377 if (EVALUE_IS_ZERO(*e))
378 return evalue_zero();
380 assert(value_zero_p(e->d) && e->x.p->type == partition);
381 res = ALLOC(evalue);
382 value_init(res->d);
383 res->x.p = new_enode(partition, e->x.p->size, e->x.p->pos);
384 for (i = 0; i < e->x.p->size/2; ++i) {
385 EVALUE_SET_DOMAIN(res->x.p->arr[2*i],
386 Domain_Copy(EVALUE_DOMAIN(e->x.p->arr[2*i])));
387 extract_term_into(&e->x.p->arr[2*i+1], var, exp,
388 &res->x.p->arr[2*i+1]);
389 reduce_evalue(&res->x.p->arr[2*i+1]);
391 return res;
394 /* Insert n free variables at pos (0-based) in the polyhedron P.
396 static Polyhedron *Polyhedron_Insert_Columns(Polyhedron *P, unsigned pos,
397 unsigned n)
399 int i;
400 unsigned NbConstraints = 0;
401 unsigned NbRays = 0;
402 Polyhedron *Q;
404 if (!P)
405 return P;
406 if (n == 0)
407 return P;
409 assert(pos <= P->Dimension);
411 if (POL_HAS(P, POL_INEQUALITIES))
412 NbConstraints = P->NbConstraints;
413 if (POL_HAS(P, POL_POINTS))
414 NbRays = P->NbRays + n;
416 Q = Polyhedron_Alloc(P->Dimension+n, NbConstraints, NbRays);
417 if (POL_HAS(P, POL_INEQUALITIES)) {
418 Q->NbEq = P->NbEq;
419 for (i = 0; i < P->NbConstraints; ++i) {
420 Vector_Copy(P->Constraint[i], Q->Constraint[i], 1+pos);
421 Vector_Copy(P->Constraint[i]+1+pos, Q->Constraint[i]+1+pos+n,
422 P->Dimension-pos+1);
425 if (POL_HAS(P, POL_POINTS)) {
426 Q->NbBid = P->NbBid + n;
427 for (i = 0; i < n; ++i)
428 value_set_si(Q->Ray[i][1+pos+i], 1);
429 for (i = 0; i < P->NbRays; ++i) {
430 Vector_Copy(P->Ray[i], Q->Ray[n+i], 1+pos);
431 Vector_Copy(P->Ray[i]+1+pos, Q->Ray[n+i]+1+pos+n,
432 P->Dimension-pos+1);
435 POL_SET(Q, POL_VALID);
436 if (POL_HAS(P, POL_INEQUALITIES))
437 POL_SET(Q, POL_INEQUALITIES);
438 if (POL_HAS(P, POL_POINTS))
439 POL_SET(Q, POL_POINTS);
440 if (POL_HAS(P, POL_VERTICES))
441 POL_SET(Q, POL_VERTICES);
442 return Q;
445 /* Perform summation of e over a list of 1 or more factors F, with context C.
446 * nvar is the total number of variables in the remaining factors.
447 * extra is the number of placeholder parameters introduced in e,
448 * but not (yet) in F or C.
450 * If there is only one factor left, F is intersected with the
451 * context C, the placeholder variables are added, and then
452 * e is summed over the resulting parametric polytope.
454 * If there is more than one factor left, we create two polynomials
455 * in a new placeholder variable (which is placed after the regular
456 * parameters, but before any previously introduced placeholder
457 * variables) that has the factors of the variables in the first
458 * factor of F and the factor of the remaining variables of
459 * each term as its coefficients.
460 * These two polynomials are then summed over their domains
461 * and afterwards the results are combined and the placeholder
462 * variable is removed again.
464 static evalue *sum_factors(Polyhedron *F, Polyhedron *C, evalue *e,
465 unsigned nvar, unsigned extra,
466 struct barvinok_options *options)
468 Polyhedron *P = F;
469 unsigned nparam = C->Dimension;
470 unsigned F_var = F->Dimension - C->Dimension;
471 int i, n;
472 evalue *s1, *s2, *s;
473 evalue *ph;
475 if (!F->next) {
476 Polyhedron *CA = align_context(C, nvar+nparam, options->MaxRays);
477 Polyhedron *P = DomainIntersection(F, CA, options->MaxRays);
478 Polyhedron *Q = Polyhedron_Insert_Columns(P, nvar+nparam, extra);
479 Polyhedron_Free(CA);
480 Polyhedron_Free(F);
481 Polyhedron_Free(P);
482 evalue *sum = sum_base(Q, e, nvar, options);
483 Polyhedron_Free(Q);
484 return sum;
487 n = evalue_count_terms(e, F_var, 0);
488 ph = create_placeholder(n, nvar+nparam);
489 evalue_shift_variables(e, nvar+nparam, 1);
490 evalue_unzip_terms(e, ph, F_var, 0);
491 evalue_shift_variables(e, nvar, -(nvar-F_var));
492 evalue_reorder_terms(ph);
493 evalue_shift_variables(ph, 0, -F_var);
495 s2 = sum_factors(F->next, C, ph, nvar-F_var, extra+1, options);
496 evalue_free(ph);
497 F->next = NULL;
498 s1 = sum_factors(F, C, e, F_var, extra+1, options);
500 if (n == 1) {
501 /* remove placeholder "polynomial" */
502 reduce_evalue(s1);
503 emul(s1, s2);
504 evalue_free(s1);
505 drop_parameters(s2, nparam, 1);
506 return s2;
509 s = evalue_zero();
510 for (i = 0; i < n; ++i) {
511 evalue *t1, *t2;
512 t1 = extract_term(s1, nparam, i);
513 t2 = extract_term(s2, nparam, i);
514 emul(t1, t2);
515 eadd(t2, s);
516 evalue_free(t1);
517 evalue_free(t2);
519 evalue_free(s1);
520 evalue_free(s2);
522 drop_parameters(s, nparam, 1);
523 return s;
526 /* Perform summation over a product of factors F, obtained using
527 * variable transformation T from the original problem specification.
529 * We first perform the corresponding transformation on the polynomial E,
530 * compute the common context over all factors and then perform
531 * the actual summation over the factors.
533 static evalue *sum_product(Polyhedron *F, evalue *E, Matrix *T, unsigned nparam,
534 struct barvinok_options *options)
536 int i;
537 Matrix *T2;
538 unsigned nvar = T->NbRows;
539 Polyhedron *C;
540 evalue *sum;
542 assert(nvar == T->NbColumns);
543 T2 = Matrix_Alloc(nvar+1, nvar+1);
544 for (i = 0; i < nvar; ++i)
545 Vector_Copy(T->p[i], T2->p[i], nvar);
546 value_set_si(T2->p[nvar][nvar], 1);
548 transform_polynomial(E, T2, NULL, nvar, nparam, nvar, nparam);
550 C = Factor_Context(F, nparam, options->MaxRays);
551 if (F->Dimension == nparam) {
552 Polyhedron *T = F;
553 F = F->next;
554 T->next = NULL;
555 Polyhedron_Free(T);
557 sum = sum_factors(F, C, E, nvar, 0, options);
559 Polyhedron_Free(C);
560 Matrix_Free(T2);
561 Matrix_Free(T);
562 return sum;
565 /* Add two constraints corresponding to floor = floor(e/d),
567 * e - d t >= 0
568 * -e + d t + d-1 >= 0
570 * e is assumed to be an affine expression.
572 Polyhedron *add_floor_var(Polyhedron *P, unsigned nvar, const evalue *floor,
573 struct barvinok_options *options)
575 int i;
576 unsigned dim = P->Dimension+1;
577 Matrix *M = Matrix_Alloc(P->NbConstraints+2, 2+dim);
578 Polyhedron *CP;
579 Value *d = &M->p[0][1+nvar];
580 evalue_extract_affine(floor, M->p[0]+1, M->p[0]+1+dim, d);
581 value_oppose(*d, *d);
582 value_set_si(M->p[0][0], 1);
583 value_set_si(M->p[1][0], 1);
584 Vector_Oppose(M->p[0]+1, M->p[1]+1, M->NbColumns-1);
585 value_subtract(M->p[1][1+dim], M->p[1][1+dim], *d);
586 value_decrement(M->p[1][1+dim], M->p[1][1+dim]);
588 for (i = 0; i < P->NbConstraints; ++i) {
589 Vector_Copy(P->Constraint[i], M->p[i+2], 1+nvar);
590 Vector_Copy(P->Constraint[i]+1+nvar, M->p[i+2]+1+nvar+1, dim-nvar-1+1);
593 CP = Constraints2Polyhedron(M, options->MaxRays);
594 Matrix_Free(M);
595 return CP;
598 static evalue *evalue_add(evalue *a, evalue *b)
600 if (!a)
601 return b;
602 if (!b)
603 return a;
604 eadd(a, b);
605 evalue_free(a);
606 return b;
609 /* Compute sum of a step-polynomial over a polytope by grouping
610 * terms containing the same floor-expressions and introducing
611 * new variables for each such expression.
612 * In particular, while there is any floor-expression left,
613 * the step-polynomial is split into a polynomial containing
614 * the expression, which is then converted to a new variable,
615 * and a polynomial not containing the expression.
617 static evalue *sum_step_polynomial(Polyhedron *P, evalue *E, unsigned nvar,
618 struct barvinok_options *options)
620 evalue *floor;
621 evalue *cur = E;
622 evalue *sum = NULL;
623 evalue *t;
625 while ((floor = evalue_outer_floor(cur))) {
626 Polyhedron *CP;
627 evalue *converted;
628 evalue *converted_floor;
630 /* Ignore floors that do not depend on variables. */
631 if (value_notzero_p(floor->d) || floor->x.p->pos >= nvar+1)
632 break;
634 converted = evalue_dup(cur);
635 converted_floor = evalue_dup(floor);
636 evalue_shift_variables(converted, nvar, 1);
637 evalue_shift_variables(converted_floor, nvar, 1);
638 evalue_replace_floor(converted, converted_floor, nvar);
639 CP = add_floor_var(P, nvar, converted_floor, options);
640 evalue_free(converted_floor);
641 t = sum_step_polynomial(CP, converted, nvar+1, options);
642 evalue_free(converted);
643 Polyhedron_Free(CP);
644 sum = evalue_add(t, sum);
646 if (cur == E)
647 cur = evalue_dup(cur);
648 evalue_drop_floor(cur, floor);
649 evalue_free(floor);
651 if (floor) {
652 evalue_floor2frac(cur);
653 evalue_free(floor);
656 if (EVALUE_IS_ZERO(*cur))
657 t = NULL;
658 else {
659 Matrix *T;
660 unsigned nparam = P->Dimension - nvar;
661 Polyhedron *F = Polyhedron_Factor(P, nparam, &T, options->MaxRays);
662 if (!F)
663 t = sum_base(P, cur, nvar, options);
664 else {
665 if (cur == E)
666 cur = evalue_dup(cur);
667 t = sum_product(F, cur, T, nparam, options);
671 if (E != cur)
672 evalue_free(cur);
674 return evalue_add(t, sum);
677 evalue *barvinok_sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
678 struct evalue_section_array *sections,
679 struct barvinok_options *options)
681 if (P->NbEq)
682 return sum_over_polytope_with_equalities(P, E, nvar, sections, options);
684 if (nvar == 0)
685 return sum_over_polytope_0D(Polyhedron_Copy(P), evalue_dup(E));
687 if (options->summation == BV_SUM_BERNOULLI)
688 return bernoulli_summate(P, E, nvar, sections, options);
689 else if (options->summation == BV_SUM_BOX)
690 return box_summate(P, E, nvar, options->MaxRays);
692 evalue_frac2floor2(E, 0);
694 return sum_step_polynomial(P, E, nvar, options);
697 evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options)
699 int i;
700 struct evalue_section_array sections;
701 evalue *sum;
703 assert(nvar >= 0);
704 if (nvar == 0 || EVALUE_IS_ZERO(*e))
705 return evalue_dup(e);
707 assert(value_zero_p(e->d));
708 assert(e->x.p->type == partition);
710 evalue_section_array_init(&sections);
711 sum = evalue_zero();
713 for (i = 0; i < e->x.p->size/2; ++i) {
714 Polyhedron *D;
715 for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
716 Polyhedron *next = D->next;
717 evalue *tmp;
718 D->next = NULL;
720 tmp = barvinok_sum_over_polytope(D, &e->x.p->arr[2*i+1], nvar,
721 &sections, options);
722 assert(tmp);
723 eadd(tmp, sum);
724 evalue_free(tmp);
726 D->next = next;
730 free(sections.s);
732 reduce_evalue(sum);
733 return sum;
736 static __isl_give isl_pw_qpolynomial *add_unbounded_guarded_qp(
737 __isl_take isl_pw_qpolynomial *sum,
738 __isl_take isl_basic_set *bset, __isl_take isl_qpolynomial *qp)
740 int zero;
742 if (!sum || !bset || !qp)
743 goto error;
745 zero = isl_qpolynomial_is_zero(qp);
746 if (zero < 0)
747 goto error;
749 if (!zero) {
750 isl_space *dim;
751 isl_set *set;
752 isl_pw_qpolynomial *pwqp;
754 dim = isl_pw_qpolynomial_get_domain_space(sum);
755 set = isl_set_from_basic_set(isl_basic_set_copy(bset));
756 set = isl_map_domain(isl_map_from_range(set));
757 set = isl_set_reset_space(set, isl_space_copy(dim));
758 pwqp = isl_pw_qpolynomial_alloc(set, isl_qpolynomial_nan_on_domain(dim));
759 sum = isl_pw_qpolynomial_add(sum, pwqp);
762 isl_basic_set_free(bset);
763 isl_qpolynomial_free(qp);
764 return sum;
765 error:
766 isl_basic_set_free(bset);
767 isl_qpolynomial_free(qp);
768 isl_pw_qpolynomial_free(sum);
769 return NULL;
772 struct barvinok_summate_data {
773 isl_space *dim;
774 __isl_take isl_qpolynomial *qp;
775 isl_pw_qpolynomial *sum;
776 int n_in;
777 int wrapping;
778 evalue *e;
779 struct evalue_section_array sections;
780 struct barvinok_options *options;
783 static int add_basic_guarded_qp(__isl_take isl_basic_set *bset, void *user)
785 struct barvinok_summate_data *data = user;
786 Polyhedron *P;
787 evalue *tmp;
788 isl_pw_qpolynomial *pwqp;
789 int bounded;
790 unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
791 unsigned nvar = isl_basic_set_dim(bset, isl_dim_set);
792 isl_space *dim;
794 if (!bset)
795 return -1;
797 bounded = isl_basic_set_is_bounded(bset);
798 if (bounded < 0)
799 goto error;
801 if (!bounded) {
802 data->sum = add_unbounded_guarded_qp(data->sum, bset,
803 isl_qpolynomial_copy(data->qp));
804 return 0;
807 dim = isl_basic_set_get_space(bset);
808 dim = isl_space_domain(isl_space_from_range(dim));
810 P = isl_basic_set_to_polylib(bset);
811 tmp = barvinok_sum_over_polytope(P, data->e, nvar,
812 &data->sections, data->options);
813 Polyhedron_Free(P);
814 assert(tmp);
815 pwqp = isl_pw_qpolynomial_from_evalue(dim, tmp);
816 evalue_free(tmp);
817 pwqp = isl_pw_qpolynomial_reset_domain_space(pwqp,
818 isl_space_domain(isl_space_copy(data->dim)));
819 data->sum = isl_pw_qpolynomial_add(data->sum, pwqp);
821 isl_basic_set_free(bset);
823 return 0;
824 error:
825 isl_basic_set_free(bset);
826 return -1;
829 static int add_guarded_qp(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
830 void *user)
832 int r;
833 struct barvinok_summate_data *data = user;
835 if (!set || !qp)
836 goto error;
838 data->qp = qp;
840 if (data->wrapping) {
841 unsigned nparam = isl_set_dim(set, isl_dim_param);
842 isl_qpolynomial *qp2 = isl_qpolynomial_copy(qp);
843 set = isl_set_move_dims(set, isl_dim_param, nparam,
844 isl_dim_set, 0, data->n_in);
845 qp2 = isl_qpolynomial_move_dims(qp2, isl_dim_param, nparam,
846 isl_dim_in, 0, data->n_in);
847 data->e = isl_qpolynomial_to_evalue(qp2);
848 isl_qpolynomial_free(qp2);
849 } else
850 data->e = isl_qpolynomial_to_evalue(qp);
851 if (!data->e)
852 goto error;
854 evalue_section_array_init(&data->sections);
856 set = isl_set_make_disjoint(set);
857 set = isl_set_compute_divs(set);
859 r = isl_set_foreach_basic_set(set, &add_basic_guarded_qp, data);
861 free(data->sections.s);
863 evalue_free(data->e);
865 isl_set_free(set);
866 isl_qpolynomial_free(qp);
868 return r;
869 error:
870 isl_set_free(set);
871 isl_qpolynomial_free(qp);
872 return -1;
875 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sum(
876 __isl_take isl_pw_qpolynomial *pwqp)
878 isl_ctx *ctx;
879 struct barvinok_summate_data data;
880 int options_allocated = 0;
881 int nvar;
883 data.dim = NULL;
884 data.options = NULL;
885 data.sum = NULL;
887 if (!pwqp)
888 return NULL;
890 nvar = isl_pw_qpolynomial_dim(pwqp, isl_dim_set);
892 data.dim = isl_pw_qpolynomial_get_domain_space(pwqp);
893 if (!data.dim)
894 goto error;
895 if (isl_space_is_params(data.dim))
896 isl_die(isl_pw_qpolynomial_get_ctx(pwqp), isl_error_invalid,
897 "input polynomial has no domain", goto error);
898 data.wrapping = isl_space_is_wrapping(data.dim);
899 if (data.wrapping) {
900 data.dim = isl_space_unwrap(data.dim);
901 data.n_in = isl_space_dim(data.dim, isl_dim_in);
902 nvar = isl_space_dim(data.dim, isl_dim_out);
903 } else
904 data.n_in = 0;
906 data.dim = isl_space_domain(data.dim);
907 if (nvar == 0)
908 return isl_pw_qpolynomial_reset_domain_space(pwqp, data.dim);
910 data.dim = isl_space_from_domain(data.dim);
911 data.dim = isl_space_add_dims(data.dim, isl_dim_out, 1);
912 data.sum = isl_pw_qpolynomial_zero(isl_space_copy(data.dim));
914 ctx = isl_pw_qpolynomial_get_ctx(pwqp);
915 data.options = isl_ctx_peek_barvinok_options(ctx);
916 if (!data.options) {
917 data.options = barvinok_options_new_with_defaults();
918 options_allocated = 1;
921 if (isl_pw_qpolynomial_foreach_lifted_piece(pwqp,
922 add_guarded_qp, &data) < 0)
923 goto error;
925 if (options_allocated)
926 barvinok_options_free(data.options);
928 isl_space_free(data.dim);
930 isl_pw_qpolynomial_free(pwqp);
932 return data.sum;
933 error:
934 if (options_allocated)
935 barvinok_options_free(data.options);
936 isl_pw_qpolynomial_free(pwqp);
937 isl_space_free(data.dim);
938 isl_pw_qpolynomial_free(data.sum);
939 return NULL;
942 static int pw_qpolynomial_sum(__isl_take isl_pw_qpolynomial *pwqp, void *user)
944 isl_union_pw_qpolynomial **res = (isl_union_pw_qpolynomial **)user;
945 isl_pw_qpolynomial *sum;
947 sum = isl_pw_qpolynomial_sum(pwqp);
948 *res = isl_union_pw_qpolynomial_add_pw_qpolynomial(*res, sum);
950 return 0;
953 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_sum(
954 __isl_take isl_union_pw_qpolynomial *upwqp)
956 isl_space *dim;
957 isl_union_pw_qpolynomial *res;
959 dim = isl_union_pw_qpolynomial_get_space(upwqp);
960 res = isl_union_pw_qpolynomial_zero(dim);
961 if (isl_union_pw_qpolynomial_foreach_pw_qpolynomial(upwqp,
962 &pw_qpolynomial_sum, &res) < 0)
963 goto error;
964 isl_union_pw_qpolynomial_free(upwqp);
966 return res;
967 error:
968 isl_union_pw_qpolynomial_free(upwqp);
969 isl_union_pw_qpolynomial_free(res);
970 return NULL;
973 static int join_compatible(__isl_keep isl_space *dim1, __isl_keep isl_space *dim2)
975 int m;
976 m = isl_space_match(dim1, isl_dim_param, dim2, isl_dim_param);
977 if (m < 0 || !m)
978 return m;
979 return isl_space_tuple_match(dim1, isl_dim_out, dim2, isl_dim_in);
982 /* Compute the intersection of the range of the map and the domain
983 * of the piecewise quasipolynomial and then sum the associated
984 * quasipolynomial over all elements in this intersection.
986 * We first introduce some unconstrained dimensions in the
987 * piecewise quasipolynomial, intersect the resulting domain
988 * with the wrapped map and then compute the sum.
990 __isl_give isl_pw_qpolynomial *isl_map_apply_pw_qpolynomial(
991 __isl_take isl_map *map, __isl_take isl_pw_qpolynomial *pwqp)
993 isl_ctx *ctx;
994 isl_set *dom;
995 isl_space *map_dim;
996 isl_space *pwqp_dim;
997 unsigned n_in;
998 int ok;
1000 ctx = isl_map_get_ctx(map);
1001 if (!ctx)
1002 goto error;
1004 map_dim = isl_map_get_space(map);
1005 pwqp_dim = isl_pw_qpolynomial_get_space(pwqp);
1006 ok = join_compatible(map_dim, pwqp_dim);
1007 isl_space_free(map_dim);
1008 isl_space_free(pwqp_dim);
1009 if (!ok)
1010 isl_die(ctx, isl_error_invalid, "incompatible dimensions",
1011 goto error);
1013 n_in = isl_map_dim(map, isl_dim_in);
1014 pwqp = isl_pw_qpolynomial_insert_dims(pwqp, isl_dim_in, 0, n_in);
1016 dom = isl_map_wrap(map);
1017 pwqp = isl_pw_qpolynomial_reset_domain_space(pwqp,
1018 isl_set_get_space(dom));
1020 pwqp = isl_pw_qpolynomial_intersect_domain(pwqp, dom);
1021 pwqp = isl_pw_qpolynomial_sum(pwqp);
1023 return pwqp;
1024 error:
1025 isl_map_free(map);
1026 isl_pw_qpolynomial_free(pwqp);
1027 return NULL;
1030 __isl_give isl_pw_qpolynomial *isl_set_apply_pw_qpolynomial(
1031 __isl_take isl_set *set, __isl_take isl_pw_qpolynomial *pwqp)
1033 isl_map *map;
1035 map = isl_map_from_range(set);
1036 pwqp = isl_map_apply_pw_qpolynomial(map, pwqp);
1037 pwqp = isl_pw_qpolynomial_project_domain_on_params(pwqp);
1038 return pwqp;
1041 struct barvinok_apply_data {
1042 isl_union_pw_qpolynomial *upwqp;
1043 isl_union_pw_qpolynomial *res;
1044 isl_map *map;
1047 static int pw_qpolynomial_apply(__isl_take isl_pw_qpolynomial *pwqp, void *user)
1049 isl_space *map_dim;
1050 isl_space *pwqp_dim;
1051 struct barvinok_apply_data *data = user;
1052 int ok;
1054 map_dim = isl_map_get_space(data->map);
1055 pwqp_dim = isl_pw_qpolynomial_get_space(pwqp);
1056 ok = join_compatible(map_dim, pwqp_dim);
1057 isl_space_free(map_dim);
1058 isl_space_free(pwqp_dim);
1060 if (ok) {
1061 pwqp = isl_map_apply_pw_qpolynomial(isl_map_copy(data->map),
1062 pwqp);
1063 data->res = isl_union_pw_qpolynomial_add_pw_qpolynomial(
1064 data->res, pwqp);
1065 } else
1066 isl_pw_qpolynomial_free(pwqp);
1068 return 0;
1071 static int map_apply(__isl_take isl_map *map, void *user)
1073 struct barvinok_apply_data *data = user;
1074 int r;
1076 data->map = map;
1077 r = isl_union_pw_qpolynomial_foreach_pw_qpolynomial(data->upwqp,
1078 &pw_qpolynomial_apply, data);
1080 isl_map_free(map);
1081 return r;
1084 __isl_give isl_union_pw_qpolynomial *isl_union_map_apply_union_pw_qpolynomial(
1085 __isl_take isl_union_map *umap,
1086 __isl_take isl_union_pw_qpolynomial *upwqp)
1088 isl_space *dim;
1089 struct barvinok_apply_data data;
1091 upwqp = isl_union_pw_qpolynomial_align_params(upwqp,
1092 isl_union_map_get_space(umap));
1093 umap = isl_union_map_align_params(umap,
1094 isl_union_pw_qpolynomial_get_space(upwqp));
1096 data.upwqp = upwqp;
1097 dim = isl_union_pw_qpolynomial_get_space(upwqp);
1098 data.res = isl_union_pw_qpolynomial_zero(dim);
1099 if (isl_union_map_foreach_map(umap, &map_apply, &data) < 0)
1100 goto error;
1102 isl_union_map_free(umap);
1103 isl_union_pw_qpolynomial_free(upwqp);
1105 return data.res;
1106 error:
1107 isl_union_map_free(umap);
1108 isl_union_pw_qpolynomial_free(upwqp);
1109 isl_union_pw_qpolynomial_free(data.res);
1110 return NULL;
1113 struct barvinok_apply_set_data {
1114 isl_union_pw_qpolynomial *upwqp;
1115 isl_union_pw_qpolynomial *res;
1116 isl_set *set;
1119 static int pw_qpolynomial_apply_set(__isl_take isl_pw_qpolynomial *pwqp,
1120 void *user)
1122 isl_space *set_dim;
1123 isl_space *pwqp_dim;
1124 struct barvinok_apply_set_data *data = user;
1125 int ok;
1127 set_dim = isl_set_get_space(data->set);
1128 pwqp_dim = isl_pw_qpolynomial_get_space(pwqp);
1129 ok = join_compatible(set_dim, pwqp_dim);
1130 isl_space_free(set_dim);
1131 isl_space_free(pwqp_dim);
1133 if (ok) {
1134 pwqp = isl_set_apply_pw_qpolynomial(isl_set_copy(data->set),
1135 pwqp);
1136 data->res = isl_union_pw_qpolynomial_add_pw_qpolynomial(
1137 data->res, pwqp);
1138 } else
1139 isl_pw_qpolynomial_free(pwqp);
1141 return 0;
1144 static int set_apply(__isl_take isl_set *set, void *user)
1146 struct barvinok_apply_set_data *data = user;
1147 int r;
1149 data->set = set;
1150 r = isl_union_pw_qpolynomial_foreach_pw_qpolynomial(data->upwqp,
1151 &pw_qpolynomial_apply_set, data);
1153 isl_set_free(set);
1154 return r;
1157 __isl_give isl_union_pw_qpolynomial *isl_union_set_apply_union_pw_qpolynomial(
1158 __isl_take isl_union_set *uset,
1159 __isl_take isl_union_pw_qpolynomial *upwqp)
1161 isl_space *dim;
1162 struct barvinok_apply_set_data data;
1164 upwqp = isl_union_pw_qpolynomial_align_params(upwqp,
1165 isl_union_set_get_space(uset));
1166 uset = isl_union_set_align_params(uset,
1167 isl_union_pw_qpolynomial_get_space(upwqp));
1169 data.upwqp = upwqp;
1170 dim = isl_union_pw_qpolynomial_get_space(upwqp);
1171 data.res = isl_union_pw_qpolynomial_zero(dim);
1172 if (isl_union_set_foreach_set(uset, &set_apply, &data) < 0)
1173 goto error;
1175 isl_union_set_free(uset);
1176 isl_union_pw_qpolynomial_free(upwqp);
1178 return data.res;
1179 error:
1180 isl_union_set_free(uset);
1181 isl_union_pw_qpolynomial_free(upwqp);
1182 isl_union_pw_qpolynomial_free(data.res);
1183 return NULL;
1186 evalue *evalue_sum(evalue *E, int nvar, unsigned MaxRays)
1188 evalue *sum;
1189 struct barvinok_options *options = barvinok_options_new_with_defaults();
1190 options->MaxRays = MaxRays;
1191 sum = barvinok_summate(E, nvar, options);
1192 barvinok_options_free(options);
1193 return sum;
1196 evalue *esum(evalue *e, int nvar)
1198 evalue *sum;
1199 struct barvinok_options *options = barvinok_options_new_with_defaults();
1200 sum = barvinok_summate(e, nvar, options);
1201 barvinok_options_free(options);
1202 return sum;
1205 /* Turn unweighted counting problem into "weighted" counting problem
1206 * with weight equal to 1 and call barvinok_summate on this weighted problem.
1208 evalue *barvinok_summate_unweighted(Polyhedron *P, Polyhedron *C,
1209 struct barvinok_options *options)
1211 Polyhedron *CA, *D;
1212 evalue e;
1213 evalue *sum;
1215 if (emptyQ(P) || emptyQ(C))
1216 return evalue_zero();
1218 CA = align_context(C, P->Dimension, options->MaxRays);
1219 D = DomainIntersection(P, CA, options->MaxRays);
1220 Domain_Free(CA);
1222 if (emptyQ(D)) {
1223 Domain_Free(D);
1224 return evalue_zero();
1227 value_init(e.d);
1228 e.x.p = new_enode(partition, 2, P->Dimension);
1229 EVALUE_SET_DOMAIN(e.x.p->arr[0], D);
1230 evalue_set_si(&e.x.p->arr[1], 1, 1);
1231 sum = barvinok_summate(&e, P->Dimension - C->Dimension, options);
1232 free_evalue_refs(&e);
1233 return sum;