update cloog to version 0.17.0
[barvinok.git] / summate.c
blobf61671c042c03b465d420ab0ebd94d2d2c996d37
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 static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E,
60 unsigned nvar,
61 struct evalue_section_array *sections,
62 struct barvinok_options *options)
64 unsigned dim = P->Dimension;
65 unsigned new_dim, new_nparam;
66 Matrix *T = NULL, *CP = NULL;
67 evalue *sum;
69 if (emptyQ(P))
70 return evalue_zero();
72 assert(P->NbEq > 0);
74 remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays);
76 if (emptyQ(P)) {
77 Polyhedron_Free(P);
78 return evalue_zero();
81 new_nparam = CP ? CP->NbColumns-1 : dim - nvar;
82 new_dim = T ? T->NbColumns-1 : nvar + new_nparam;
84 /* We can avoid these substitutions if E is a constant */
85 E = evalue_dup(E);
86 transform_polynomial(E, T, CP, nvar, dim-nvar,
87 new_dim-new_nparam, new_nparam);
89 if (new_dim-new_nparam > 0) {
90 sum = barvinok_sum_over_polytope(P, E, new_dim-new_nparam,
91 sections, options);
92 evalue_free(E);
93 Polyhedron_Free(P);
94 } else {
95 sum = ALLOC(evalue);
96 value_init(sum->d);
97 sum->x.p = new_enode(partition, 2, new_dim);
98 EVALUE_SET_DOMAIN(sum->x.p->arr[0], P);
99 value_clear(sum->x.p->arr[1].d);
100 sum->x.p->arr[1] = *E;
101 free(E);
104 if (CP) {
105 evalue_backsubstitute(sum, CP, options->MaxRays);
106 Matrix_Free(CP);
109 if (T)
110 Matrix_Free(T);
112 return sum;
115 static evalue *sum_base(Polyhedron *P, evalue *E, unsigned nvar,
116 struct barvinok_options *options)
118 if (options->summation == BV_SUM_EULER)
119 return euler_summate(P, E, nvar, options);
120 else if (options->summation == BV_SUM_LAURENT)
121 return laurent_summate(P, E, nvar, options);
122 else if (options->summation == BV_SUM_LAURENT_OLD)
123 return laurent_summate_old(P, E, nvar, options);
124 assert(0);
127 /* Count the number of non-zero terms in e when viewed as a polynomial
128 * in only the first nvar variables. "count" is the number counted
129 * so far.
131 static int evalue_count_terms(const evalue *e, unsigned nvar, int count)
133 int i;
135 if (EVALUE_IS_ZERO(*e))
136 return count;
138 if (value_zero_p(e->d))
139 assert(e->x.p->type == polynomial);
140 if (value_notzero_p(e->d) || e->x.p->pos >= nvar+1)
141 return count+1;
143 for (i = 0; i < e->x.p->size; ++i)
144 count = evalue_count_terms(&e->x.p->arr[i], nvar, count);
146 return count;
149 /* Create placeholder structure for unzipping.
150 * A "polynomial" is created with size terms in variable pos,
151 * with each term having itself as coefficient.
153 static evalue *create_placeholder(int size, int pos)
155 int i, j;
156 evalue *E = ALLOC(evalue);
157 value_init(E->d);
158 E->x.p = new_enode(polynomial, size, pos+1);
159 for (i = 0; i < size; ++i) {
160 E->x.p->arr[i].x.p = new_enode(polynomial, i+1, pos+1);
161 for (j = 0; j < i; ++j)
162 evalue_set_si(&E->x.p->arr[i].x.p->arr[j], 0, 1);
163 evalue_set_si(&E->x.p->arr[i].x.p->arr[i], 1, 1);
165 return E;
168 /* Interchange each non-zero term in e (when viewed as a polynomial
169 * in only the first nvar variables) with a placeholder in ph (created
170 * by create_placeholder), resulting in two polynomials in the
171 * placeholder variable such that for each non-zero term in e
172 * there is a power of the placeholder variable such that the factors
173 * in the first nvar variables form the coefficient of that power in
174 * the first polynomial (e) and the factors in the remaining variables
175 * form the coefficient of that power in the second polynomial (ph).
177 static int evalue_unzip_terms(evalue *e, evalue *ph, 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 evalue t = *e;
188 *e = ph->x.p->arr[count];
189 ph->x.p->arr[count] = t;
190 return count+1;
193 for (i = 0; i < e->x.p->size; ++i)
194 count = evalue_unzip_terms(&e->x.p->arr[i], ph, nvar, count);
196 return count;
199 /* Remove n variables at pos (0-based) from the polyhedron P.
200 * Each of these variables is assumed to be completely free,
201 * i.e., there is a line in the polyhedron corresponding to
202 * each of these variables.
204 static Polyhedron *Polyhedron_Remove_Columns(Polyhedron *P, unsigned pos,
205 unsigned n)
207 int i, j;
208 unsigned NbConstraints = 0;
209 unsigned NbRays = 0;
210 Polyhedron *Q;
212 if (n == 0)
213 return P;
215 assert(pos <= P->Dimension);
217 if (POL_HAS(P, POL_INEQUALITIES))
218 NbConstraints = P->NbConstraints;
219 if (POL_HAS(P, POL_POINTS))
220 NbRays = P->NbRays - n;
222 Q = Polyhedron_Alloc(P->Dimension - n, NbConstraints, NbRays);
223 if (POL_HAS(P, POL_INEQUALITIES)) {
224 Q->NbEq = P->NbEq;
225 for (i = 0; i < P->NbConstraints; ++i) {
226 Vector_Copy(P->Constraint[i], Q->Constraint[i], 1+pos);
227 Vector_Copy(P->Constraint[i]+1+pos+n, Q->Constraint[i]+1+pos,
228 Q->Dimension-pos+1);
231 if (POL_HAS(P, POL_POINTS)) {
232 Q->NbBid = P->NbBid - n;
233 for (i = 0; i < n; ++i)
234 value_set_si(Q->Ray[i][1+pos+i], 1);
235 for (i = 0, j = 0; i < P->NbRays; ++i) {
236 int line = First_Non_Zero(P->Ray[i], 1+P->Dimension+1);
237 assert(line != -1);
238 if (line-1 >= pos && line-1 < pos+n) {
239 ++j;
240 assert(j <= n);
241 continue;
243 assert(i-j < Q->NbRays);
244 Vector_Copy(P->Ray[i], Q->Ray[i-j], 1+pos);
245 Vector_Copy(P->Ray[i]+1+pos+n, Q->Ray[i-j]+1+pos,
246 Q->Dimension-pos+1);
249 POL_SET(Q, POL_VALID);
250 if (POL_HAS(P, POL_INEQUALITIES))
251 POL_SET(Q, POL_INEQUALITIES);
252 if (POL_HAS(P, POL_POINTS))
253 POL_SET(Q, POL_POINTS);
254 if (POL_HAS(P, POL_VERTICES))
255 POL_SET(Q, POL_VERTICES);
256 return Q;
259 /* Remove n variables at pos (0-based) from the union of polyhedra P.
260 * Each of these variables is assumed to be completely free,
261 * i.e., there is a line in the polyhedron corresponding to
262 * each of these variables.
264 static Polyhedron *Domain_Remove_Columns(Polyhedron *P, unsigned pos,
265 unsigned n)
267 Polyhedron *R;
268 Polyhedron **next = &R;
270 for (; P; P = P->next) {
271 *next = Polyhedron_Remove_Columns(P, pos, n);
272 next = &(*next)->next;
274 return R;
277 /* Drop n parameters starting at first from partition evalue e */
278 static void drop_parameters(evalue *e, int first, int n)
280 int i;
282 if (EVALUE_IS_ZERO(*e))
283 return;
285 assert(value_zero_p(e->d) && e->x.p->type == partition);
286 for (i = 0; i < e->x.p->size/2; ++i) {
287 Polyhedron *P = EVALUE_DOMAIN(e->x.p->arr[2*i]);
288 Polyhedron *Q = Domain_Remove_Columns(P, first, n);
289 EVALUE_SET_DOMAIN(e->x.p->arr[2*i], Q);
290 Domain_Free(P);
291 evalue_shift_variables(&e->x.p->arr[2*i+1], first, -n);
293 e->x.p->pos -= n;
296 static void extract_term_into(const evalue *src, int var, int exp, evalue *dst)
298 int i;
300 if (value_notzero_p(src->d) ||
301 src->x.p->type != polynomial ||
302 src->x.p->pos > var+1) {
303 if (exp == 0)
304 evalue_copy(dst, src);
305 else
306 evalue_set_si(dst, 0, 1);
307 return;
310 if (src->x.p->pos == var+1) {
311 if (src->x.p->size > exp)
312 evalue_copy(dst, &src->x.p->arr[exp]);
313 else
314 evalue_set_si(dst, 0, 1);
315 return;
318 dst->x.p = new_enode(polynomial, src->x.p->size, src->x.p->pos);
319 for (i = 0; i < src->x.p->size; ++i)
320 extract_term_into(&src->x.p->arr[i], var, exp,
321 &dst->x.p->arr[i]);
324 /* Extract the coefficient of var^exp.
326 static evalue *extract_term(const evalue *e, int var, int exp)
328 int i;
329 evalue *res;
331 if (EVALUE_IS_ZERO(*e))
332 return evalue_zero();
334 assert(value_zero_p(e->d) && e->x.p->type == partition);
335 res = ALLOC(evalue);
336 value_init(res->d);
337 res->x.p = new_enode(partition, e->x.p->size, e->x.p->pos);
338 for (i = 0; i < e->x.p->size/2; ++i) {
339 EVALUE_SET_DOMAIN(res->x.p->arr[2*i],
340 Domain_Copy(EVALUE_DOMAIN(e->x.p->arr[2*i])));
341 extract_term_into(&e->x.p->arr[2*i+1], var, exp,
342 &res->x.p->arr[2*i+1]);
343 reduce_evalue(&res->x.p->arr[2*i+1]);
345 return res;
348 /* Insert n free variables at pos (0-based) in the polyhedron P.
350 static Polyhedron *Polyhedron_Insert_Columns(Polyhedron *P, unsigned pos,
351 unsigned n)
353 int i;
354 unsigned NbConstraints = 0;
355 unsigned NbRays = 0;
356 Polyhedron *Q;
358 if (!P)
359 return P;
360 if (n == 0)
361 return P;
363 assert(pos <= P->Dimension);
365 if (POL_HAS(P, POL_INEQUALITIES))
366 NbConstraints = P->NbConstraints;
367 if (POL_HAS(P, POL_POINTS))
368 NbRays = P->NbRays + n;
370 Q = Polyhedron_Alloc(P->Dimension+n, NbConstraints, NbRays);
371 if (POL_HAS(P, POL_INEQUALITIES)) {
372 Q->NbEq = P->NbEq;
373 for (i = 0; i < P->NbConstraints; ++i) {
374 Vector_Copy(P->Constraint[i], Q->Constraint[i], 1+pos);
375 Vector_Copy(P->Constraint[i]+1+pos, Q->Constraint[i]+1+pos+n,
376 P->Dimension-pos+1);
379 if (POL_HAS(P, POL_POINTS)) {
380 Q->NbBid = P->NbBid + n;
381 for (i = 0; i < n; ++i)
382 value_set_si(Q->Ray[i][1+pos+i], 1);
383 for (i = 0; i < P->NbRays; ++i) {
384 Vector_Copy(P->Ray[i], Q->Ray[n+i], 1+pos);
385 Vector_Copy(P->Ray[i]+1+pos, Q->Ray[n+i]+1+pos+n,
386 P->Dimension-pos+1);
389 POL_SET(Q, POL_VALID);
390 if (POL_HAS(P, POL_INEQUALITIES))
391 POL_SET(Q, POL_INEQUALITIES);
392 if (POL_HAS(P, POL_POINTS))
393 POL_SET(Q, POL_POINTS);
394 if (POL_HAS(P, POL_VERTICES))
395 POL_SET(Q, POL_VERTICES);
396 return Q;
399 /* Perform summation of e over a list of 1 or more factors F, with context C.
400 * nvar is the total number of variables in the remaining factors.
401 * extra is the number of placeholder parameters introduced in e,
402 * but not (yet) in F or C.
404 * If there is only one factor left, F is intersected with the
405 * context C, the placeholder variables are added, and then
406 * e is summed over the resulting parametric polytope.
408 * If there is more than one factor left, we create two polynomials
409 * in a new placeholder variable (which is placed after the regular
410 * parameters, but before any previously introduced placeholder
411 * variables) that has the factors of the variables in the first
412 * factor of F and the factor of the remaining variables of
413 * each term as its coefficients.
414 * These two polynomials are then summed over their domains
415 * and afterwards the results are combined and the placeholder
416 * variable is removed again.
418 static evalue *sum_factors(Polyhedron *F, Polyhedron *C, evalue *e,
419 unsigned nvar, unsigned extra,
420 struct barvinok_options *options)
422 Polyhedron *P = F;
423 unsigned nparam = C->Dimension;
424 unsigned F_var = F->Dimension - C->Dimension;
425 int i, n;
426 evalue *s1, *s2, *s;
427 evalue *ph;
429 if (!F->next) {
430 Polyhedron *CA = align_context(C, nvar+nparam, options->MaxRays);
431 Polyhedron *P = DomainIntersection(F, CA, options->MaxRays);
432 Polyhedron *Q = Polyhedron_Insert_Columns(P, nvar+nparam, extra);
433 Polyhedron_Free(CA);
434 Polyhedron_Free(F);
435 Polyhedron_Free(P);
436 evalue *sum = sum_base(Q, e, nvar, options);
437 Polyhedron_Free(Q);
438 return sum;
441 n = evalue_count_terms(e, F_var, 0);
442 ph = create_placeholder(n, nvar+nparam);
443 evalue_shift_variables(e, nvar+nparam, 1);
444 evalue_unzip_terms(e, ph, F_var, 0);
445 evalue_shift_variables(e, nvar, -(nvar-F_var));
446 evalue_reorder_terms(ph);
447 evalue_shift_variables(ph, 0, -F_var);
449 s2 = sum_factors(F->next, C, ph, nvar-F_var, extra+1, options);
450 evalue_free(ph);
451 F->next = NULL;
452 s1 = sum_factors(F, C, e, F_var, extra+1, options);
454 if (n == 1) {
455 /* remove placeholder "polynomial" */
456 reduce_evalue(s1);
457 emul(s1, s2);
458 evalue_free(s1);
459 drop_parameters(s2, nparam, 1);
460 return s2;
463 s = evalue_zero();
464 for (i = 0; i < n; ++i) {
465 evalue *t1, *t2;
466 t1 = extract_term(s1, nparam, i);
467 t2 = extract_term(s2, nparam, i);
468 emul(t1, t2);
469 eadd(t2, s);
470 evalue_free(t1);
471 evalue_free(t2);
473 evalue_free(s1);
474 evalue_free(s2);
476 drop_parameters(s, nparam, 1);
477 return s;
480 /* Perform summation over a product of factors F, obtained using
481 * variable transformation T from the original problem specification.
483 * We first perform the corresponding transformation on the polynomial E,
484 * compute the common context over all factors and then perform
485 * the actual summation over the factors.
487 static evalue *sum_product(Polyhedron *F, evalue *E, Matrix *T, unsigned nparam,
488 struct barvinok_options *options)
490 int i;
491 Matrix *T2;
492 unsigned nvar = T->NbRows;
493 Polyhedron *C;
494 evalue *sum;
496 assert(nvar == T->NbColumns);
497 T2 = Matrix_Alloc(nvar+1, nvar+1);
498 for (i = 0; i < nvar; ++i)
499 Vector_Copy(T->p[i], T2->p[i], nvar);
500 value_set_si(T2->p[nvar][nvar], 1);
502 transform_polynomial(E, T2, NULL, nvar, nparam, nvar, nparam);
504 C = Factor_Context(F, nparam, options->MaxRays);
505 if (F->Dimension == nparam) {
506 Polyhedron *T = F;
507 F = F->next;
508 T->next = NULL;
509 Polyhedron_Free(T);
511 sum = sum_factors(F, C, E, nvar, 0, options);
513 Polyhedron_Free(C);
514 Matrix_Free(T2);
515 Matrix_Free(T);
516 return sum;
519 /* Add two constraints corresponding to floor = floor(e/d),
521 * e - d t >= 0
522 * -e + d t + d-1 >= 0
524 * e is assumed to be an affine expression.
526 Polyhedron *add_floor_var(Polyhedron *P, unsigned nvar, const evalue *floor,
527 struct barvinok_options *options)
529 int i;
530 unsigned dim = P->Dimension+1;
531 Matrix *M = Matrix_Alloc(P->NbConstraints+2, 2+dim);
532 Polyhedron *CP;
533 Value *d = &M->p[0][1+nvar];
534 evalue_extract_affine(floor, M->p[0]+1, M->p[0]+1+dim, d);
535 value_oppose(*d, *d);
536 value_set_si(M->p[0][0], 1);
537 value_set_si(M->p[1][0], 1);
538 Vector_Oppose(M->p[0]+1, M->p[1]+1, M->NbColumns-1);
539 value_subtract(M->p[1][1+dim], M->p[1][1+dim], *d);
540 value_decrement(M->p[1][1+dim], M->p[1][1+dim]);
542 for (i = 0; i < P->NbConstraints; ++i) {
543 Vector_Copy(P->Constraint[i], M->p[i+2], 1+nvar);
544 Vector_Copy(P->Constraint[i]+1+nvar, M->p[i+2]+1+nvar+1, dim-nvar-1+1);
547 CP = Constraints2Polyhedron(M, options->MaxRays);
548 Matrix_Free(M);
549 return CP;
552 static evalue *evalue_add(evalue *a, evalue *b)
554 if (!a)
555 return b;
556 if (!b)
557 return a;
558 eadd(a, b);
559 evalue_free(a);
560 return b;
563 /* Compute sum of a step-polynomial over a polytope by grouping
564 * terms containing the same floor-expressions and introducing
565 * new variables for each such expression.
566 * In particular, while there is any floor-expression left,
567 * the step-polynomial is split into a polynomial containing
568 * the expression, which is then converted to a new variable,
569 * and a polynomial not containing the expression.
571 static evalue *sum_step_polynomial(Polyhedron *P, evalue *E, unsigned nvar,
572 struct barvinok_options *options)
574 evalue *floor;
575 evalue *cur = E;
576 evalue *sum = NULL;
577 evalue *t;
579 while ((floor = evalue_outer_floor(cur))) {
580 Polyhedron *CP;
581 evalue *converted;
582 evalue *converted_floor;
584 /* Ignore floors that do not depend on variables. */
585 if (value_notzero_p(floor->d) || floor->x.p->pos >= nvar+1)
586 break;
588 converted = evalue_dup(cur);
589 converted_floor = evalue_dup(floor);
590 evalue_shift_variables(converted, nvar, 1);
591 evalue_shift_variables(converted_floor, nvar, 1);
592 evalue_replace_floor(converted, converted_floor, nvar);
593 CP = add_floor_var(P, nvar, converted_floor, options);
594 evalue_free(converted_floor);
595 t = sum_step_polynomial(CP, converted, nvar+1, options);
596 evalue_free(converted);
597 Polyhedron_Free(CP);
598 sum = evalue_add(t, sum);
600 if (cur == E)
601 cur = evalue_dup(cur);
602 evalue_drop_floor(cur, floor);
603 evalue_free(floor);
605 if (floor) {
606 evalue_floor2frac(cur);
607 evalue_free(floor);
610 if (EVALUE_IS_ZERO(*cur))
611 t = NULL;
612 else {
613 Matrix *T;
614 unsigned nparam = P->Dimension - nvar;
615 Polyhedron *F = Polyhedron_Factor(P, nparam, &T, options->MaxRays);
616 if (!F)
617 t = sum_base(P, cur, nvar, options);
618 else {
619 if (cur == E)
620 cur = evalue_dup(cur);
621 t = sum_product(F, cur, T, nparam, options);
625 if (E != cur)
626 evalue_free(cur);
628 return evalue_add(t, sum);
631 evalue *barvinok_sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
632 struct evalue_section_array *sections,
633 struct barvinok_options *options)
635 if (P->NbEq)
636 return sum_over_polytope_with_equalities(P, E, nvar, sections, options);
638 if (options->summation == BV_SUM_BERNOULLI)
639 return bernoulli_summate(P, E, nvar, sections, options);
640 else if (options->summation == BV_SUM_BOX)
641 return box_summate(P, E, nvar, options->MaxRays);
643 evalue_frac2floor2(E, 0);
645 return sum_step_polynomial(P, E, nvar, options);
648 evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options)
650 int i;
651 struct evalue_section_array sections;
652 evalue *sum;
654 assert(nvar >= 0);
655 if (nvar == 0 || EVALUE_IS_ZERO(*e))
656 return evalue_dup(e);
658 assert(value_zero_p(e->d));
659 assert(e->x.p->type == partition);
661 evalue_section_array_init(&sections);
662 sum = evalue_zero();
664 for (i = 0; i < e->x.p->size/2; ++i) {
665 Polyhedron *D;
666 for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
667 Polyhedron *next = D->next;
668 evalue *tmp;
669 D->next = NULL;
671 tmp = barvinok_sum_over_polytope(D, &e->x.p->arr[2*i+1], nvar,
672 &sections, options);
673 assert(tmp);
674 eadd(tmp, sum);
675 evalue_free(tmp);
677 D->next = next;
681 free(sections.s);
683 reduce_evalue(sum);
684 return sum;
687 static __isl_give isl_pw_qpolynomial *add_unbounded_guarded_qp(
688 __isl_take isl_pw_qpolynomial *sum,
689 __isl_take isl_basic_set *bset, __isl_take isl_qpolynomial *qp)
691 int zero;
693 if (!sum || !bset || !qp)
694 goto error;
696 zero = isl_qpolynomial_is_zero(qp);
697 if (zero < 0)
698 goto error;
700 if (!zero) {
701 isl_space *dim;
702 isl_set *set;
703 isl_pw_qpolynomial *pwqp;
705 dim = isl_pw_qpolynomial_get_space(sum);
706 set = isl_set_from_basic_set(isl_basic_set_copy(bset));
707 set = isl_map_domain(isl_map_from_range(set));
708 set = isl_set_reset_space(set, isl_space_copy(dim));
709 pwqp = isl_pw_qpolynomial_alloc(set, isl_qpolynomial_nan_on_domain(dim));
710 sum = isl_pw_qpolynomial_add(sum, pwqp);
713 isl_basic_set_free(bset);
714 isl_qpolynomial_free(qp);
715 return sum;
716 error:
717 isl_basic_set_free(bset);
718 isl_qpolynomial_free(qp);
719 isl_pw_qpolynomial_free(sum);
720 return NULL;
723 struct barvinok_summate_data {
724 isl_space *dim;
725 __isl_take isl_qpolynomial *qp;
726 isl_pw_qpolynomial *sum;
727 int n_in;
728 int wrapping;
729 evalue *e;
730 struct evalue_section_array sections;
731 struct barvinok_options *options;
734 static int add_basic_guarded_qp(__isl_take isl_basic_set *bset, void *user)
736 struct barvinok_summate_data *data = user;
737 Polyhedron *P;
738 evalue *tmp;
739 isl_pw_qpolynomial *pwqp;
740 int bounded;
741 unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
742 unsigned nvar = isl_basic_set_dim(bset, isl_dim_set);
743 isl_space *dim;
745 if (!bset)
746 return -1;
748 bounded = isl_basic_set_is_bounded(bset);
749 if (bounded < 0)
750 goto error;
752 if (!bounded) {
753 data->sum = add_unbounded_guarded_qp(data->sum, bset,
754 isl_qpolynomial_copy(data->qp));
755 return 0;
758 dim = isl_basic_set_get_space(bset);
759 dim = isl_space_domain(isl_space_from_range(dim));
761 P = isl_basic_set_to_polylib(bset);
762 tmp = barvinok_sum_over_polytope(P, data->e, nvar,
763 &data->sections, data->options);
764 Polyhedron_Free(P);
765 assert(tmp);
766 pwqp = isl_pw_qpolynomial_from_evalue(dim, tmp);
767 evalue_free(tmp);
768 pwqp = isl_pw_qpolynomial_reset_domain_space(pwqp,
769 isl_space_domain(isl_space_copy(data->dim)));
770 data->sum = isl_pw_qpolynomial_add(data->sum, pwqp);
772 isl_basic_set_free(bset);
774 return 0;
775 error:
776 isl_basic_set_free(bset);
777 return -1;
780 static int add_guarded_qp(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
781 void *user)
783 int r;
784 struct barvinok_summate_data *data = user;
786 if (!set || !qp)
787 goto error;
789 data->qp = qp;
791 if (data->wrapping) {
792 unsigned nparam = isl_set_dim(set, isl_dim_param);
793 isl_qpolynomial *qp2 = isl_qpolynomial_copy(qp);
794 set = isl_set_move_dims(set, isl_dim_param, nparam,
795 isl_dim_set, 0, data->n_in);
796 qp2 = isl_qpolynomial_move_dims(qp2, isl_dim_param, nparam,
797 isl_dim_in, 0, data->n_in);
798 data->e = isl_qpolynomial_to_evalue(qp2);
799 isl_qpolynomial_free(qp2);
800 } else
801 data->e = isl_qpolynomial_to_evalue(qp);
802 if (!data->e)
803 goto error;
805 evalue_section_array_init(&data->sections);
807 set = isl_set_make_disjoint(set);
808 set = isl_set_compute_divs(set);
810 r = isl_set_foreach_basic_set(set, &add_basic_guarded_qp, data);
812 free(data->sections.s);
814 evalue_free(data->e);
816 isl_set_free(set);
817 isl_qpolynomial_free(qp);
819 return r;
820 error:
821 isl_set_free(set);
822 isl_qpolynomial_free(qp);
823 return -1;
826 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sum(
827 __isl_take isl_pw_qpolynomial *pwqp)
829 isl_ctx *ctx;
830 struct barvinok_summate_data data;
831 int options_allocated = 0;
832 int nvar;
834 data.dim = NULL;
835 data.options = NULL;
837 if (!pwqp)
838 return NULL;
840 nvar = isl_pw_qpolynomial_dim(pwqp, isl_dim_set);
842 data.dim = isl_pw_qpolynomial_get_domain_space(pwqp);
843 data.wrapping = isl_space_is_wrapping(data.dim);
844 if (data.wrapping) {
845 data.dim = isl_space_unwrap(data.dim);
846 data.n_in = isl_space_dim(data.dim, isl_dim_in);
847 nvar = isl_space_dim(data.dim, isl_dim_out);
848 } else
849 data.n_in = 0;
851 data.dim = isl_space_domain(data.dim);
852 if (nvar == 0)
853 return isl_pw_qpolynomial_reset_domain_space(pwqp, data.dim);
855 data.dim = isl_space_from_domain(data.dim);
856 data.dim = isl_space_add_dims(data.dim, isl_dim_out, 1);
857 data.sum = isl_pw_qpolynomial_zero(isl_space_copy(data.dim));
859 ctx = isl_pw_qpolynomial_get_ctx(pwqp);
860 data.options = isl_ctx_peek_barvinok_options(ctx);
861 if (!data.options) {
862 data.options = barvinok_options_new_with_defaults();
863 options_allocated = 1;
866 if (isl_pw_qpolynomial_foreach_lifted_piece(pwqp,
867 add_guarded_qp, &data) < 0)
868 goto error;
870 if (options_allocated)
871 barvinok_options_free(data.options);
873 isl_space_free(data.dim);
875 isl_pw_qpolynomial_free(pwqp);
877 return data.sum;
878 error:
879 if (options_allocated)
880 barvinok_options_free(data.options);
881 isl_pw_qpolynomial_free(pwqp);
882 isl_space_free(data.dim);
883 isl_pw_qpolynomial_free(data.sum);
884 return NULL;
887 static int pw_qpolynomial_sum(__isl_take isl_pw_qpolynomial *pwqp, void *user)
889 isl_union_pw_qpolynomial **res = (isl_union_pw_qpolynomial **)user;
890 isl_pw_qpolynomial *sum;
892 sum = isl_pw_qpolynomial_sum(pwqp);
893 *res = isl_union_pw_qpolynomial_add_pw_qpolynomial(*res, sum);
895 return 0;
898 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_sum(
899 __isl_take isl_union_pw_qpolynomial *upwqp)
901 isl_space *dim;
902 isl_union_pw_qpolynomial *res;
904 dim = isl_union_pw_qpolynomial_get_space(upwqp);
905 res = isl_union_pw_qpolynomial_zero(dim);
906 if (isl_union_pw_qpolynomial_foreach_pw_qpolynomial(upwqp,
907 &pw_qpolynomial_sum, &res) < 0)
908 goto error;
909 isl_union_pw_qpolynomial_free(upwqp);
911 return res;
912 error:
913 isl_union_pw_qpolynomial_free(upwqp);
914 isl_union_pw_qpolynomial_free(res);
915 return NULL;
918 static int join_compatible(__isl_keep isl_space *dim1, __isl_keep isl_space *dim2)
920 int m;
921 m = isl_space_match(dim1, isl_dim_param, dim2, isl_dim_param);
922 if (m < 0 || !m)
923 return m;
924 return isl_space_tuple_match(dim1, isl_dim_out, dim2, isl_dim_in);
927 /* Compute the intersection of the range of the map and the domain
928 * of the piecewise quasipolynomial and then sum the associated
929 * quasipolynomial over all elements in this intersection.
931 * We first introduce some unconstrained dimensions in the
932 * piecewise quasipolynomial, intersect the resulting domain
933 * with the wrapped map and then compute the sum.
935 __isl_give isl_pw_qpolynomial *isl_map_apply_pw_qpolynomial(
936 __isl_take isl_map *map, __isl_take isl_pw_qpolynomial *pwqp)
938 isl_ctx *ctx;
939 isl_set *dom;
940 isl_space *map_dim;
941 isl_space *pwqp_dim;
942 unsigned n_in;
943 int ok;
945 ctx = isl_map_get_ctx(map);
946 if (!ctx)
947 goto error;
949 map_dim = isl_map_get_space(map);
950 pwqp_dim = isl_pw_qpolynomial_get_space(pwqp);
951 ok = join_compatible(map_dim, pwqp_dim);
952 isl_space_free(map_dim);
953 isl_space_free(pwqp_dim);
954 if (!ok)
955 isl_die(ctx, isl_error_invalid, "incompatible dimensions",
956 goto error);
958 n_in = isl_map_dim(map, isl_dim_in);
959 pwqp = isl_pw_qpolynomial_insert_dims(pwqp, isl_dim_in, 0, n_in);
961 dom = isl_map_wrap(map);
962 pwqp = isl_pw_qpolynomial_reset_domain_space(pwqp,
963 isl_set_get_space(dom));
965 pwqp = isl_pw_qpolynomial_intersect_domain(pwqp, dom);
966 pwqp = isl_pw_qpolynomial_sum(pwqp);
968 return pwqp;
969 error:
970 isl_map_free(map);
971 isl_pw_qpolynomial_free(pwqp);
972 return NULL;
975 __isl_give isl_pw_qpolynomial *isl_set_apply_pw_qpolynomial(
976 __isl_take isl_set *set, __isl_take isl_pw_qpolynomial *pwqp)
978 isl_map *map;
980 map = isl_map_from_range(set);
981 pwqp = isl_map_apply_pw_qpolynomial(map, pwqp);
982 pwqp = isl_pw_qpolynomial_project_domain_on_params(pwqp);
983 return pwqp;
986 struct barvinok_apply_data {
987 isl_union_pw_qpolynomial *upwqp;
988 isl_union_pw_qpolynomial *res;
989 isl_map *map;
992 static int pw_qpolynomial_apply(__isl_take isl_pw_qpolynomial *pwqp, void *user)
994 isl_space *map_dim;
995 isl_space *pwqp_dim;
996 struct barvinok_apply_data *data = user;
997 int ok;
999 map_dim = isl_map_get_space(data->map);
1000 pwqp_dim = isl_pw_qpolynomial_get_space(pwqp);
1001 ok = join_compatible(map_dim, pwqp_dim);
1002 isl_space_free(map_dim);
1003 isl_space_free(pwqp_dim);
1005 if (ok) {
1006 pwqp = isl_map_apply_pw_qpolynomial(isl_map_copy(data->map),
1007 pwqp);
1008 data->res = isl_union_pw_qpolynomial_add_pw_qpolynomial(
1009 data->res, pwqp);
1010 } else
1011 isl_pw_qpolynomial_free(pwqp);
1013 return 0;
1016 static int map_apply(__isl_take isl_map *map, void *user)
1018 struct barvinok_apply_data *data = user;
1019 int r;
1021 data->map = map;
1022 r = isl_union_pw_qpolynomial_foreach_pw_qpolynomial(data->upwqp,
1023 &pw_qpolynomial_apply, data);
1025 isl_map_free(map);
1026 return r;
1029 __isl_give isl_union_pw_qpolynomial *isl_union_map_apply_union_pw_qpolynomial(
1030 __isl_take isl_union_map *umap,
1031 __isl_take isl_union_pw_qpolynomial *upwqp)
1033 isl_space *dim;
1034 struct barvinok_apply_data data;
1036 upwqp = isl_union_pw_qpolynomial_align_params(upwqp,
1037 isl_union_map_get_space(umap));
1038 umap = isl_union_map_align_params(umap,
1039 isl_union_pw_qpolynomial_get_space(upwqp));
1041 data.upwqp = upwqp;
1042 dim = isl_union_pw_qpolynomial_get_space(upwqp);
1043 data.res = isl_union_pw_qpolynomial_zero(dim);
1044 if (isl_union_map_foreach_map(umap, &map_apply, &data) < 0)
1045 goto error;
1047 isl_union_map_free(umap);
1048 isl_union_pw_qpolynomial_free(upwqp);
1050 return data.res;
1051 error:
1052 isl_union_map_free(umap);
1053 isl_union_pw_qpolynomial_free(upwqp);
1054 isl_union_pw_qpolynomial_free(data.res);
1055 return NULL;
1058 struct barvinok_apply_set_data {
1059 isl_union_pw_qpolynomial *upwqp;
1060 isl_union_pw_qpolynomial *res;
1061 isl_set *set;
1064 static int pw_qpolynomial_apply_set(__isl_take isl_pw_qpolynomial *pwqp,
1065 void *user)
1067 isl_space *set_dim;
1068 isl_space *pwqp_dim;
1069 struct barvinok_apply_set_data *data = user;
1070 int ok;
1072 set_dim = isl_set_get_space(data->set);
1073 pwqp_dim = isl_pw_qpolynomial_get_space(pwqp);
1074 ok = join_compatible(set_dim, pwqp_dim);
1075 isl_space_free(set_dim);
1076 isl_space_free(pwqp_dim);
1078 if (ok) {
1079 pwqp = isl_set_apply_pw_qpolynomial(isl_set_copy(data->set),
1080 pwqp);
1081 data->res = isl_union_pw_qpolynomial_add_pw_qpolynomial(
1082 data->res, pwqp);
1083 } else
1084 isl_pw_qpolynomial_free(pwqp);
1086 return 0;
1089 static int set_apply(__isl_take isl_set *set, void *user)
1091 struct barvinok_apply_set_data *data = user;
1092 int r;
1094 data->set = set;
1095 r = isl_union_pw_qpolynomial_foreach_pw_qpolynomial(data->upwqp,
1096 &pw_qpolynomial_apply_set, data);
1098 isl_set_free(set);
1099 return r;
1102 __isl_give isl_union_pw_qpolynomial *isl_union_set_apply_union_pw_qpolynomial(
1103 __isl_take isl_union_set *uset,
1104 __isl_take isl_union_pw_qpolynomial *upwqp)
1106 isl_space *dim;
1107 struct barvinok_apply_set_data data;
1109 upwqp = isl_union_pw_qpolynomial_align_params(upwqp,
1110 isl_union_set_get_space(uset));
1111 uset = isl_union_set_align_params(uset,
1112 isl_union_pw_qpolynomial_get_space(upwqp));
1114 data.upwqp = upwqp;
1115 dim = isl_union_pw_qpolynomial_get_space(upwqp);
1116 data.res = isl_union_pw_qpolynomial_zero(dim);
1117 if (isl_union_set_foreach_set(uset, &set_apply, &data) < 0)
1118 goto error;
1120 isl_union_set_free(uset);
1121 isl_union_pw_qpolynomial_free(upwqp);
1123 return data.res;
1124 error:
1125 isl_union_set_free(uset);
1126 isl_union_pw_qpolynomial_free(upwqp);
1127 isl_union_pw_qpolynomial_free(data.res);
1128 return NULL;
1131 evalue *evalue_sum(evalue *E, int nvar, unsigned MaxRays)
1133 evalue *sum;
1134 struct barvinok_options *options = barvinok_options_new_with_defaults();
1135 options->MaxRays = MaxRays;
1136 sum = barvinok_summate(E, nvar, options);
1137 barvinok_options_free(options);
1138 return sum;
1141 evalue *esum(evalue *e, int nvar)
1143 evalue *sum;
1144 struct barvinok_options *options = barvinok_options_new_with_defaults();
1145 sum = barvinok_summate(e, nvar, options);
1146 barvinok_options_free(options);
1147 return sum;
1150 /* Turn unweighted counting problem into "weighted" counting problem
1151 * with weight equal to 1 and call barvinok_summate on this weighted problem.
1153 evalue *barvinok_summate_unweighted(Polyhedron *P, Polyhedron *C,
1154 struct barvinok_options *options)
1156 Polyhedron *CA, *D;
1157 evalue e;
1158 evalue *sum;
1160 if (emptyQ(P) || emptyQ(C))
1161 return evalue_zero();
1163 CA = align_context(C, P->Dimension, options->MaxRays);
1164 D = DomainIntersection(P, CA, options->MaxRays);
1165 Domain_Free(CA);
1167 if (emptyQ(D)) {
1168 Domain_Free(D);
1169 return evalue_zero();
1172 value_init(e.d);
1173 e.x.p = new_enode(partition, 2, P->Dimension);
1174 EVALUE_SET_DOMAIN(e.x.p->arr[0], D);
1175 evalue_set_si(&e.x.p->arr[1], 1, 1);
1176 sum = barvinok_summate(&e, P->Dimension - C->Dimension, options);
1177 free_evalue_refs(&e);
1178 return sum;