configure.ac: rename BV_LDFLAGS to BV_LIB_LDFLAGS
[barvinok.git] / summate.c
blobf0c000434b2fb241487d34de0389ee9c59245c09
1 #include <isl/space.h>
2 #include <isl/set.h>
3 #include <isl/map.h>
4 #include <isl/union_set.h>
5 #include <isl/union_map.h>
6 #include <isl/polynomial.h>
7 #include <isl_set_polylib.h>
8 #include <barvinok/polylib.h>
9 #include <barvinok/options.h>
10 #include <barvinok/util.h>
11 #include "bernoulli.h"
12 #include "euler.h"
13 #include "laurent.h"
14 #include "laurent_old.h"
15 #include "summate.h"
16 #include "section_array.h"
17 #include "remove_equalities.h"
19 extern evalue *evalue_outer_floor(evalue *e);
20 extern int evalue_replace_floor(evalue *e, const evalue *floor, int var);
21 extern void evalue_drop_floor(evalue *e, const evalue *floor);
23 #define ALLOC(type) (type*)malloc(sizeof(type))
24 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
26 /* Apply the variable transformation specified by T and CP on
27 * the polynomial e. T expresses the old variables in terms
28 * of the new variables (and optionally also the new parameters),
29 * while CP expresses the old parameters in terms of the new
30 * parameters.
32 static void transform_polynomial(evalue *E, Matrix *T, Matrix *CP,
33 unsigned nvar, unsigned nparam,
34 unsigned new_nvar, unsigned new_nparam)
36 int j;
37 evalue **subs;
39 subs = ALLOCN(evalue *, nvar+nparam);
41 for (j = 0; j < nvar; ++j) {
42 if (T)
43 subs[j] = affine2evalue(T->p[j], T->p[T->NbRows-1][T->NbColumns-1],
44 T->NbColumns-1);
45 else
46 subs[j] = evalue_var(j);
48 for (j = 0; j < nparam; ++j) {
49 if (CP)
50 subs[nvar+j] = affine2evalue(CP->p[j], CP->p[nparam][new_nparam],
51 new_nparam);
52 else
53 subs[nvar+j] = evalue_var(j);
54 evalue_shift_variables(subs[nvar+j], 0, new_nvar);
57 evalue_substitute(E, subs);
58 reduce_evalue(E);
60 for (j = 0; j < nvar+nparam; ++j)
61 evalue_free(subs[j]);
62 free(subs);
65 /* Compute the sum of the quasi-polynomial E
66 * over a 0D (non-empty, but possibly parametric) polytope P.
68 * Consumes P and E.
70 * We simply return a partition evalue with P as domain and E as value.
72 static evalue *sum_over_polytope_0D(Polyhedron *P, evalue *E)
74 evalue *sum;
76 sum = ALLOC(evalue);
77 value_init(sum->d);
78 sum->x.p = new_enode(partition, 2, P->Dimension);
79 EVALUE_SET_DOMAIN(sum->x.p->arr[0], P);
80 value_clear(sum->x.p->arr[1].d);
81 sum->x.p->arr[1] = *E;
82 free(E);
84 return sum;
87 static evalue *sum_with_equalities(Polyhedron *P, evalue *E,
88 unsigned nvar, struct evalue_section_array *sections,
89 struct barvinok_options *options,
90 evalue *(*base)(Polyhedron *P, evalue *E, unsigned nvar,
91 struct evalue_section_array *sections,
92 struct barvinok_options *options))
94 unsigned dim = P->Dimension;
95 unsigned new_dim, new_nparam;
96 Matrix *T = NULL, *CP = NULL;
97 evalue *sum;
99 if (emptyQ(P))
100 return evalue_zero();
102 assert(P->NbEq > 0);
104 remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays);
106 if (emptyQ(P)) {
107 Polyhedron_Free(P);
108 return evalue_zero();
111 new_nparam = CP ? CP->NbColumns-1 : dim - nvar;
112 new_dim = T ? T->NbColumns-1 : nvar + new_nparam;
114 /* We can avoid these substitutions if E is a constant */
115 E = evalue_dup(E);
116 transform_polynomial(E, T, CP, nvar, dim-nvar,
117 new_dim-new_nparam, new_nparam);
119 if (new_dim-new_nparam > 0) {
120 sum = base(P, E, new_dim-new_nparam, sections, options);
121 evalue_free(E);
122 Polyhedron_Free(P);
123 } else {
124 sum = sum_over_polytope_0D(P, E);
127 if (CP) {
128 evalue_backsubstitute(sum, CP, options->MaxRays);
129 Matrix_Free(CP);
132 if (T)
133 Matrix_Free(T);
135 return sum;
138 static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E,
139 unsigned nvar, struct evalue_section_array *sections,
140 struct barvinok_options *options)
142 return sum_with_equalities(P, E, nvar, sections, options,
143 &barvinok_sum_over_polytope);
146 static evalue *sum_base(Polyhedron *P, evalue *E, unsigned nvar,
147 struct barvinok_options *options);
149 static evalue *sum_base_wrap(Polyhedron *P, evalue *E, unsigned nvar,
150 struct evalue_section_array *sections, struct barvinok_options *options)
152 return sum_base(P, E, nvar, options);
155 static evalue *sum_base_with_equalities(Polyhedron *P, evalue *E,
156 unsigned nvar, struct barvinok_options *options)
158 return sum_with_equalities(P, E, nvar, NULL, options, &sum_base_wrap);
161 /* The substitutions in sum_step_polynomial may have reintroduced equalities
162 * (in particular, one of the floor expressions may be equal to one of
163 * the variables), so we need to check for them again.
165 static evalue *sum_base(Polyhedron *P, evalue *E, unsigned nvar,
166 struct barvinok_options *options)
168 if (P->NbEq)
169 return sum_base_with_equalities(P, E, nvar, options);
170 if (options->summation == BV_SUM_EULER)
171 return euler_summate(P, E, nvar, options);
172 else if (options->summation == BV_SUM_LAURENT)
173 return laurent_summate(P, E, nvar, options);
174 else if (options->summation == BV_SUM_LAURENT_OLD)
175 return laurent_summate_old(P, E, nvar, options);
176 assert(0);
179 /* Count the number of non-zero terms in e when viewed as a polynomial
180 * in only the first nvar variables. "count" is the number counted
181 * so far.
183 static int evalue_count_terms(const evalue *e, unsigned nvar, int count)
185 int i;
187 if (EVALUE_IS_ZERO(*e))
188 return count;
190 if (value_zero_p(e->d))
191 assert(e->x.p->type == polynomial);
192 if (value_notzero_p(e->d) || e->x.p->pos >= nvar+1)
193 return count+1;
195 for (i = 0; i < e->x.p->size; ++i)
196 count = evalue_count_terms(&e->x.p->arr[i], nvar, count);
198 return count;
201 /* Create placeholder structure for unzipping.
202 * A "polynomial" is created with size terms in variable pos,
203 * with each term having itself as coefficient.
205 static evalue *create_placeholder(int size, int pos)
207 int i, j;
208 evalue *E = ALLOC(evalue);
209 value_init(E->d);
210 E->x.p = new_enode(polynomial, size, pos+1);
211 for (i = 0; i < size; ++i) {
212 E->x.p->arr[i].x.p = new_enode(polynomial, i+1, pos+1);
213 for (j = 0; j < i; ++j)
214 evalue_set_si(&E->x.p->arr[i].x.p->arr[j], 0, 1);
215 evalue_set_si(&E->x.p->arr[i].x.p->arr[i], 1, 1);
217 return E;
220 /* Interchange each non-zero term in e (when viewed as a polynomial
221 * in only the first nvar variables) with a placeholder in ph (created
222 * by create_placeholder), resulting in two polynomials in the
223 * placeholder variable such that for each non-zero term in e
224 * there is a power of the placeholder variable such that the factors
225 * in the first nvar variables form the coefficient of that power in
226 * the first polynomial (e) and the factors in the remaining variables
227 * form the coefficient of that power in the second polynomial (ph).
229 static int evalue_unzip_terms(evalue *e, evalue *ph, unsigned nvar, int count)
231 int i;
233 if (EVALUE_IS_ZERO(*e))
234 return count;
236 if (value_zero_p(e->d))
237 assert(e->x.p->type == polynomial);
238 if (value_notzero_p(e->d) || e->x.p->pos >= nvar+1) {
239 evalue t = *e;
240 *e = ph->x.p->arr[count];
241 ph->x.p->arr[count] = t;
242 return count+1;
245 for (i = 0; i < e->x.p->size; ++i)
246 count = evalue_unzip_terms(&e->x.p->arr[i], ph, nvar, count);
248 return count;
251 /* Remove n variables at pos (0-based) from the polyhedron P.
252 * Each of these variables is assumed to be completely free,
253 * i.e., there is a line in the polyhedron corresponding to
254 * each of these variables.
256 static Polyhedron *Polyhedron_Remove_Columns(Polyhedron *P, unsigned pos,
257 unsigned n)
259 int i, j;
260 unsigned NbConstraints = 0;
261 unsigned NbRays = 0;
262 Polyhedron *Q;
264 if (n == 0)
265 return P;
267 assert(pos <= P->Dimension);
269 if (POL_HAS(P, POL_INEQUALITIES))
270 NbConstraints = P->NbConstraints;
271 if (POL_HAS(P, POL_POINTS))
272 NbRays = P->NbRays - n;
274 Q = Polyhedron_Alloc(P->Dimension - n, NbConstraints, NbRays);
275 if (POL_HAS(P, POL_INEQUALITIES)) {
276 Q->NbEq = P->NbEq;
277 for (i = 0; i < P->NbConstraints; ++i) {
278 Vector_Copy(P->Constraint[i], Q->Constraint[i], 1+pos);
279 Vector_Copy(P->Constraint[i]+1+pos+n, Q->Constraint[i]+1+pos,
280 Q->Dimension-pos+1);
283 if (POL_HAS(P, POL_POINTS)) {
284 Q->NbBid = P->NbBid - n;
285 for (i = 0; i < n; ++i)
286 value_set_si(Q->Ray[i][1+pos+i], 1);
287 for (i = 0, j = 0; i < P->NbRays; ++i) {
288 int line = First_Non_Zero(P->Ray[i], 1+P->Dimension+1);
289 assert(line != -1);
290 if (line-1 >= pos && line-1 < pos+n) {
291 ++j;
292 assert(j <= n);
293 continue;
295 assert(i-j < Q->NbRays);
296 Vector_Copy(P->Ray[i], Q->Ray[i-j], 1+pos);
297 Vector_Copy(P->Ray[i]+1+pos+n, Q->Ray[i-j]+1+pos,
298 Q->Dimension-pos+1);
301 POL_SET(Q, POL_VALID);
302 if (POL_HAS(P, POL_INEQUALITIES))
303 POL_SET(Q, POL_INEQUALITIES);
304 if (POL_HAS(P, POL_POINTS))
305 POL_SET(Q, POL_POINTS);
306 if (POL_HAS(P, POL_VERTICES))
307 POL_SET(Q, POL_VERTICES);
308 return Q;
311 /* Remove n variables at pos (0-based) from the union of polyhedra P.
312 * Each of these variables is assumed to be completely free,
313 * i.e., there is a line in the polyhedron corresponding to
314 * each of these variables.
316 static Polyhedron *Domain_Remove_Columns(Polyhedron *P, unsigned pos,
317 unsigned n)
319 Polyhedron *R;
320 Polyhedron **next = &R;
322 for (; P; P = P->next) {
323 *next = Polyhedron_Remove_Columns(P, pos, n);
324 next = &(*next)->next;
326 return R;
329 /* Drop n parameters starting at first from partition evalue e */
330 static void drop_parameters(evalue *e, int first, int n)
332 int i;
334 if (EVALUE_IS_ZERO(*e))
335 return;
337 assert(value_zero_p(e->d) && e->x.p->type == partition);
338 for (i = 0; i < e->x.p->size/2; ++i) {
339 Polyhedron *P = EVALUE_DOMAIN(e->x.p->arr[2*i]);
340 Polyhedron *Q = Domain_Remove_Columns(P, first, n);
341 EVALUE_SET_DOMAIN(e->x.p->arr[2*i], Q);
342 Domain_Free(P);
343 evalue_shift_variables(&e->x.p->arr[2*i+1], first, -n);
345 e->x.p->pos -= n;
348 static void extract_term_into(const evalue *src, int var, int exp, evalue *dst)
350 int i;
352 if (value_notzero_p(src->d) ||
353 src->x.p->type != polynomial ||
354 src->x.p->pos > var+1) {
355 if (exp == 0)
356 evalue_copy(dst, src);
357 else
358 evalue_set_si(dst, 0, 1);
359 return;
362 if (src->x.p->pos == var+1) {
363 if (src->x.p->size > exp)
364 evalue_copy(dst, &src->x.p->arr[exp]);
365 else
366 evalue_set_si(dst, 0, 1);
367 return;
370 dst->x.p = new_enode(polynomial, src->x.p->size, src->x.p->pos);
371 for (i = 0; i < src->x.p->size; ++i)
372 extract_term_into(&src->x.p->arr[i], var, exp,
373 &dst->x.p->arr[i]);
376 /* Extract the coefficient of var^exp.
378 static evalue *extract_term(const evalue *e, int var, int exp)
380 int i;
381 evalue *res;
383 if (EVALUE_IS_ZERO(*e))
384 return evalue_zero();
386 assert(value_zero_p(e->d) && e->x.p->type == partition);
387 res = ALLOC(evalue);
388 value_init(res->d);
389 res->x.p = new_enode(partition, e->x.p->size, e->x.p->pos);
390 for (i = 0; i < e->x.p->size/2; ++i) {
391 EVALUE_SET_DOMAIN(res->x.p->arr[2*i],
392 Domain_Copy(EVALUE_DOMAIN(e->x.p->arr[2*i])));
393 extract_term_into(&e->x.p->arr[2*i+1], var, exp,
394 &res->x.p->arr[2*i+1]);
395 reduce_evalue(&res->x.p->arr[2*i+1]);
397 return res;
400 /* Insert n free variables at pos (0-based) in the polyhedron P.
402 static Polyhedron *Polyhedron_Insert_Columns(Polyhedron *P, unsigned pos,
403 unsigned n)
405 int i;
406 unsigned NbConstraints = 0;
407 unsigned NbRays = 0;
408 Polyhedron *Q;
410 if (!P)
411 return P;
412 if (n == 0)
413 return P;
415 assert(pos <= P->Dimension);
417 if (POL_HAS(P, POL_INEQUALITIES))
418 NbConstraints = P->NbConstraints;
419 if (POL_HAS(P, POL_POINTS))
420 NbRays = P->NbRays + n;
422 Q = Polyhedron_Alloc(P->Dimension+n, NbConstraints, NbRays);
423 if (POL_HAS(P, POL_INEQUALITIES)) {
424 Q->NbEq = P->NbEq;
425 for (i = 0; i < P->NbConstraints; ++i) {
426 Vector_Copy(P->Constraint[i], Q->Constraint[i], 1+pos);
427 Vector_Copy(P->Constraint[i]+1+pos, Q->Constraint[i]+1+pos+n,
428 P->Dimension-pos+1);
431 if (POL_HAS(P, POL_POINTS)) {
432 Q->NbBid = P->NbBid + n;
433 for (i = 0; i < n; ++i)
434 value_set_si(Q->Ray[i][1+pos+i], 1);
435 for (i = 0; i < P->NbRays; ++i) {
436 Vector_Copy(P->Ray[i], Q->Ray[n+i], 1+pos);
437 Vector_Copy(P->Ray[i]+1+pos, Q->Ray[n+i]+1+pos+n,
438 P->Dimension-pos+1);
441 POL_SET(Q, POL_VALID);
442 if (POL_HAS(P, POL_INEQUALITIES))
443 POL_SET(Q, POL_INEQUALITIES);
444 if (POL_HAS(P, POL_POINTS))
445 POL_SET(Q, POL_POINTS);
446 if (POL_HAS(P, POL_VERTICES))
447 POL_SET(Q, POL_VERTICES);
448 return Q;
451 /* Perform summation of e over a list of 1 or more factors F, with context C.
452 * nvar is the total number of variables in the remaining factors.
453 * extra is the number of placeholder parameters introduced in e,
454 * but not (yet) in F or C.
456 * If there is only one factor left, F is intersected with the
457 * context C, the placeholder variables are added, and then
458 * e is summed over the resulting parametric polytope.
460 * If there is more than one factor left, we create two polynomials
461 * in a new placeholder variable (which is placed after the regular
462 * parameters, but before any previously introduced placeholder
463 * variables) that has the factors of the variables in the first
464 * factor of F and the factor of the remaining variables of
465 * each term as its coefficients.
466 * These two polynomials are then summed over their domains
467 * and afterwards the results are combined and the placeholder
468 * variable is removed again.
470 static evalue *sum_factors(Polyhedron *F, Polyhedron *C, evalue *e,
471 unsigned nvar, unsigned extra,
472 struct barvinok_options *options)
474 unsigned nparam = C->Dimension;
475 unsigned F_var = F->Dimension - C->Dimension;
476 int i, n;
477 evalue *s1, *s2, *s;
478 evalue *ph;
480 if (!F->next) {
481 Polyhedron *CA = align_context(C, nvar+nparam, options->MaxRays);
482 Polyhedron *P = DomainIntersection(F, CA, options->MaxRays);
483 Polyhedron *Q = Polyhedron_Insert_Columns(P, nvar+nparam, extra);
484 Polyhedron_Free(CA);
485 Polyhedron_Free(F);
486 Polyhedron_Free(P);
487 evalue *sum = sum_base(Q, e, nvar, options);
488 Polyhedron_Free(Q);
489 return sum;
492 n = evalue_count_terms(e, F_var, 0);
493 ph = create_placeholder(n, nvar+nparam);
494 evalue_shift_variables(e, nvar+nparam, 1);
495 evalue_unzip_terms(e, ph, F_var, 0);
496 evalue_shift_variables(e, nvar, -(nvar-F_var));
497 evalue_reorder_terms(ph);
498 evalue_shift_variables(ph, 0, -F_var);
500 s2 = sum_factors(F->next, C, ph, nvar-F_var, extra+1, options);
501 evalue_free(ph);
502 F->next = NULL;
503 s1 = sum_factors(F, C, e, F_var, extra+1, options);
505 if (n == 1) {
506 /* remove placeholder "polynomial" */
507 reduce_evalue(s1);
508 emul(s1, s2);
509 evalue_free(s1);
510 drop_parameters(s2, nparam, 1);
511 return s2;
514 s = evalue_zero();
515 for (i = 0; i < n; ++i) {
516 evalue *t1, *t2;
517 t1 = extract_term(s1, nparam, i);
518 t2 = extract_term(s2, nparam, i);
519 emul(t1, t2);
520 eadd(t2, s);
521 evalue_free(t1);
522 evalue_free(t2);
524 evalue_free(s1);
525 evalue_free(s2);
527 drop_parameters(s, nparam, 1);
528 return s;
531 /* Perform summation over a product of factors F, obtained using
532 * variable transformation T from the original problem specification.
534 * We first perform the corresponding transformation on the polynomial E,
535 * compute the common context over all factors and then perform
536 * the actual summation over the factors.
538 static evalue *sum_product(Polyhedron *F, evalue *E, Matrix *T, unsigned nparam,
539 struct barvinok_options *options)
541 int i;
542 Matrix *T2;
543 unsigned nvar = T->NbRows;
544 Polyhedron *C;
545 evalue *sum;
547 assert(nvar == T->NbColumns);
548 T2 = Matrix_Alloc(nvar+1, nvar+1);
549 for (i = 0; i < nvar; ++i)
550 Vector_Copy(T->p[i], T2->p[i], nvar);
551 value_set_si(T2->p[nvar][nvar], 1);
553 transform_polynomial(E, T2, NULL, nvar, nparam, nvar, nparam);
555 C = Factor_Context(F, nparam, options->MaxRays);
556 if (F->Dimension == nparam) {
557 Polyhedron *T = F;
558 F = F->next;
559 T->next = NULL;
560 Polyhedron_Free(T);
562 sum = sum_factors(F, C, E, nvar, 0, options);
564 Polyhedron_Free(C);
565 Matrix_Free(T2);
566 Matrix_Free(T);
567 return sum;
570 /* Add two constraints corresponding to floor = floor(e/d),
572 * e - d t >= 0
573 * -e + d t + d-1 >= 0
575 * e is assumed to be an affine expression.
577 Polyhedron *add_floor_var(Polyhedron *P, unsigned nvar, const evalue *floor,
578 struct barvinok_options *options)
580 int i;
581 unsigned dim = P->Dimension+1;
582 Matrix *M = Matrix_Alloc(P->NbConstraints+2, 2+dim);
583 Polyhedron *CP;
584 Value *d = &M->p[0][1+nvar];
585 evalue_extract_affine(floor, M->p[0]+1, M->p[0]+1+dim, d);
586 value_oppose(*d, *d);
587 value_set_si(M->p[0][0], 1);
588 value_set_si(M->p[1][0], 1);
589 Vector_Oppose(M->p[0]+1, M->p[1]+1, M->NbColumns-1);
590 value_subtract(M->p[1][1+dim], M->p[1][1+dim], *d);
591 value_decrement(M->p[1][1+dim], M->p[1][1+dim]);
593 for (i = 0; i < P->NbConstraints; ++i) {
594 Vector_Copy(P->Constraint[i], M->p[i+2], 1+nvar);
595 Vector_Copy(P->Constraint[i]+1+nvar, M->p[i+2]+1+nvar+1, dim-nvar-1+1);
598 CP = Constraints2Polyhedron(M, options->MaxRays);
599 Matrix_Free(M);
600 return CP;
603 static evalue *evalue_add(evalue *a, evalue *b)
605 if (!a)
606 return b;
607 if (!b)
608 return a;
609 eadd(a, b);
610 evalue_free(a);
611 return b;
614 /* Compute sum of a step-polynomial over a polytope by grouping
615 * terms containing the same floor-expressions and introducing
616 * new variables for each such expression.
617 * In particular, while there is any floor-expression left,
618 * the step-polynomial is split into a polynomial containing
619 * the expression, which is then converted to a new variable,
620 * and a polynomial not containing the expression.
622 static evalue *sum_step_polynomial(Polyhedron *P, evalue *E, unsigned nvar,
623 struct barvinok_options *options)
625 evalue *floor;
626 evalue *cur = E;
627 evalue *sum = NULL;
628 evalue *t;
630 while ((floor = evalue_outer_floor(cur))) {
631 Polyhedron *CP;
632 evalue *converted;
633 evalue *converted_floor;
635 /* Ignore floors that do not depend on variables. */
636 if (value_notzero_p(floor->d) || floor->x.p->pos >= nvar+1)
637 break;
639 converted = evalue_dup(cur);
640 converted_floor = evalue_dup(floor);
641 evalue_shift_variables(converted, nvar, 1);
642 evalue_shift_variables(converted_floor, nvar, 1);
643 evalue_replace_floor(converted, converted_floor, nvar);
644 CP = add_floor_var(P, nvar, converted_floor, options);
645 evalue_free(converted_floor);
646 t = sum_step_polynomial(CP, converted, nvar+1, options);
647 evalue_free(converted);
648 Polyhedron_Free(CP);
649 sum = evalue_add(t, sum);
651 if (cur == E)
652 cur = evalue_dup(cur);
653 evalue_drop_floor(cur, floor);
654 evalue_free(floor);
656 if (floor) {
657 evalue_floor2frac(cur);
658 evalue_free(floor);
661 if (EVALUE_IS_ZERO(*cur))
662 t = NULL;
663 else {
664 Matrix *T;
665 unsigned nparam = P->Dimension - nvar;
666 Polyhedron *F = Polyhedron_Factor(P, nparam, &T, options->MaxRays);
667 if (!F)
668 t = sum_base(P, cur, nvar, options);
669 else {
670 if (cur == E)
671 cur = evalue_dup(cur);
672 t = sum_product(F, cur, T, nparam, options);
676 if (E != cur)
677 evalue_free(cur);
679 return evalue_add(t, sum);
682 evalue *barvinok_sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
683 struct evalue_section_array *sections,
684 struct barvinok_options *options)
686 if (P->NbEq)
687 return sum_over_polytope_with_equalities(P, E, nvar, sections, options);
689 if (nvar == 0)
690 return sum_over_polytope_0D(Polyhedron_Copy(P), evalue_dup(E));
692 if (options->summation == BV_SUM_BERNOULLI)
693 return bernoulli_summate(P, E, nvar, sections, options);
694 else if (options->summation == BV_SUM_BOX)
695 return box_summate(P, E, nvar, options->MaxRays);
697 evalue_frac2floor2(E, 0);
699 return sum_step_polynomial(P, E, nvar, options);
702 evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options)
704 int i;
705 struct evalue_section_array sections;
706 evalue *sum;
708 assert(nvar >= 0);
709 if (nvar == 0 || EVALUE_IS_ZERO(*e))
710 return evalue_dup(e);
712 assert(value_zero_p(e->d));
713 assert(e->x.p->type == partition);
715 evalue_section_array_init(&sections);
716 sum = evalue_zero();
718 for (i = 0; i < e->x.p->size/2; ++i) {
719 Polyhedron *D;
720 for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
721 Polyhedron *next = D->next;
722 evalue *tmp;
723 D->next = NULL;
725 tmp = barvinok_sum_over_polytope(D, &e->x.p->arr[2*i+1], nvar,
726 &sections, options);
727 assert(tmp);
728 eadd(tmp, sum);
729 evalue_free(tmp);
731 D->next = next;
735 free(sections.s);
737 reduce_evalue(sum);
738 return sum;
741 static __isl_give isl_pw_qpolynomial *add_unbounded_guarded_qp(
742 __isl_take isl_pw_qpolynomial *sum,
743 __isl_take isl_basic_set *bset, __isl_take isl_qpolynomial *qp)
745 int zero;
747 if (!sum || !bset || !qp)
748 goto error;
750 zero = isl_qpolynomial_is_zero(qp);
751 if (zero < 0)
752 goto error;
754 if (!zero) {
755 isl_space *space;
756 isl_set *set;
757 isl_pw_qpolynomial *pwqp;
759 space = isl_pw_qpolynomial_get_domain_space(sum);
760 set = isl_set_from_basic_set(isl_basic_set_copy(bset));
761 set = isl_map_domain(isl_map_from_range(set));
762 set = isl_set_reset_space(set, isl_space_copy(space));
763 pwqp = isl_pw_qpolynomial_alloc(set,
764 isl_qpolynomial_nan_on_domain(space));
765 sum = isl_pw_qpolynomial_add(sum, pwqp);
768 isl_basic_set_free(bset);
769 isl_qpolynomial_free(qp);
770 return sum;
771 error:
772 isl_basic_set_free(bset);
773 isl_qpolynomial_free(qp);
774 isl_pw_qpolynomial_free(sum);
775 return NULL;
778 struct barvinok_summate_data {
779 isl_space *space;
780 isl_qpolynomial *qp;
781 isl_pw_qpolynomial *sum;
782 int n_in;
783 int wrapping;
784 evalue *e;
785 struct evalue_section_array sections;
786 struct barvinok_options *options;
789 static isl_stat add_basic_guarded_qp(__isl_take isl_basic_set *bset, void *user)
791 struct barvinok_summate_data *data = user;
792 Polyhedron *P;
793 evalue *tmp;
794 isl_pw_qpolynomial *pwqp;
795 int bounded;
796 unsigned nvar = isl_basic_set_dim(bset, isl_dim_set);
797 isl_space *space;
799 if (!bset)
800 return isl_stat_error;
802 bounded = isl_basic_set_is_bounded(bset);
803 if (bounded < 0)
804 goto error;
806 if (!bounded) {
807 data->sum = add_unbounded_guarded_qp(data->sum, bset,
808 isl_qpolynomial_copy(data->qp));
809 return isl_stat_ok;
812 space = isl_space_params(isl_basic_set_get_space(bset));
814 P = isl_basic_set_to_polylib(bset);
815 tmp = barvinok_sum_over_polytope(P, data->e, nvar,
816 &data->sections, data->options);
817 Polyhedron_Free(P);
818 assert(tmp);
819 pwqp = isl_pw_qpolynomial_from_evalue(space, tmp);
820 evalue_free(tmp);
821 pwqp = isl_pw_qpolynomial_reset_domain_space(pwqp,
822 isl_space_domain(isl_space_copy(data->space)));
823 data->sum = isl_pw_qpolynomial_add(data->sum, pwqp);
825 isl_basic_set_free(bset);
827 return isl_stat_ok;
828 error:
829 isl_basic_set_free(bset);
830 return isl_stat_error;
833 static isl_stat add_guarded_qp(__isl_take isl_set *set,
834 __isl_take isl_qpolynomial *qp, void *user)
836 isl_stat r;
837 struct barvinok_summate_data *data = user;
839 if (!set || !qp)
840 goto error;
842 data->qp = qp;
844 if (data->wrapping) {
845 unsigned nparam = isl_set_dim(set, isl_dim_param);
846 isl_qpolynomial *qp2 = isl_qpolynomial_copy(qp);
847 set = isl_set_move_dims(set, isl_dim_param, nparam,
848 isl_dim_set, 0, data->n_in);
849 qp2 = isl_qpolynomial_move_dims(qp2, isl_dim_param, nparam,
850 isl_dim_in, 0, data->n_in);
851 data->e = isl_qpolynomial_to_evalue(qp2);
852 isl_qpolynomial_free(qp2);
853 } else
854 data->e = isl_qpolynomial_to_evalue(qp);
855 if (!data->e)
856 goto error;
858 evalue_section_array_init(&data->sections);
860 set = isl_set_make_disjoint(set);
861 set = isl_set_compute_divs(set);
863 r = isl_set_foreach_basic_set(set, &add_basic_guarded_qp, data);
865 free(data->sections.s);
867 evalue_free(data->e);
869 isl_set_free(set);
870 isl_qpolynomial_free(qp);
872 return r;
873 error:
874 isl_set_free(set);
875 isl_qpolynomial_free(qp);
876 return isl_stat_error;
879 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sum(
880 __isl_take isl_pw_qpolynomial *pwqp)
882 isl_ctx *ctx;
883 struct barvinok_summate_data data;
884 int options_allocated = 0;
885 int nvar;
887 data.space = NULL;
888 data.options = NULL;
889 data.sum = NULL;
891 if (!pwqp)
892 return NULL;
894 nvar = isl_pw_qpolynomial_dim(pwqp, isl_dim_set);
896 data.space = isl_pw_qpolynomial_get_domain_space(pwqp);
897 if (!data.space)
898 goto error;
899 if (isl_space_is_params(data.space))
900 isl_die(isl_pw_qpolynomial_get_ctx(pwqp), isl_error_invalid,
901 "input polynomial has no domain", goto error);
902 data.wrapping = isl_space_is_wrapping(data.space);
903 if (data.wrapping) {
904 data.space = isl_space_unwrap(data.space);
905 data.n_in = isl_space_dim(data.space, isl_dim_in);
906 nvar = isl_space_dim(data.space, isl_dim_out);
907 } else
908 data.n_in = 0;
910 data.space = isl_space_domain(data.space);
911 if (nvar == 0)
912 return isl_pw_qpolynomial_reset_domain_space(pwqp, data.space);
914 data.space = isl_space_from_domain(data.space);
915 data.space = isl_space_add_dims(data.space, isl_dim_out, 1);
916 data.sum = isl_pw_qpolynomial_zero(isl_space_copy(data.space));
918 ctx = isl_pw_qpolynomial_get_ctx(pwqp);
919 data.options = isl_ctx_peek_barvinok_options(ctx);
920 if (!data.options) {
921 data.options = barvinok_options_new_with_defaults();
922 options_allocated = 1;
925 if (isl_pw_qpolynomial_foreach_lifted_piece(pwqp,
926 add_guarded_qp, &data) < 0)
927 goto error;
929 if (options_allocated)
930 barvinok_options_free(data.options);
932 isl_space_free(data.space);
934 isl_pw_qpolynomial_free(pwqp);
936 return data.sum;
937 error:
938 if (options_allocated)
939 barvinok_options_free(data.options);
940 isl_pw_qpolynomial_free(pwqp);
941 isl_space_free(data.space);
942 isl_pw_qpolynomial_free(data.sum);
943 return NULL;
946 static isl_stat pw_qpolynomial_sum(__isl_take isl_pw_qpolynomial *pwqp,
947 void *user)
949 isl_union_pw_qpolynomial **res = (isl_union_pw_qpolynomial **)user;
950 isl_pw_qpolynomial *sum;
951 isl_union_pw_qpolynomial *upwqp;
953 sum = isl_pw_qpolynomial_sum(pwqp);
954 upwqp = isl_union_pw_qpolynomial_from_pw_qpolynomial(sum);
955 *res = isl_union_pw_qpolynomial_add(*res, upwqp);
957 return isl_stat_ok;
960 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_sum(
961 __isl_take isl_union_pw_qpolynomial *upwqp)
963 isl_space *space;
964 isl_union_pw_qpolynomial *res;
966 space = isl_union_pw_qpolynomial_get_space(upwqp);
967 res = isl_union_pw_qpolynomial_zero(space);
968 if (isl_union_pw_qpolynomial_foreach_pw_qpolynomial(upwqp,
969 &pw_qpolynomial_sum, &res) < 0)
970 goto error;
971 isl_union_pw_qpolynomial_free(upwqp);
973 return res;
974 error:
975 isl_union_pw_qpolynomial_free(upwqp);
976 isl_union_pw_qpolynomial_free(res);
977 return NULL;
980 static int join_compatible(__isl_keep isl_space *space1,
981 __isl_keep isl_space *space2)
983 int m;
984 m = isl_space_has_equal_params(space1, space2);
985 if (m < 0 || !m)
986 return m;
987 return isl_space_tuple_is_equal(space1, isl_dim_out,
988 space2, isl_dim_in);
991 /* Compute the intersection of the range of the map and the domain
992 * of the piecewise quasipolynomial and then sum the associated
993 * quasipolynomial over all elements in this intersection.
995 * We first introduce some unconstrained dimensions in the
996 * piecewise quasipolynomial, intersect the resulting domain
997 * with the wrapped map and then compute the sum.
999 __isl_give isl_pw_qpolynomial *isl_map_apply_pw_qpolynomial(
1000 __isl_take isl_map *map, __isl_take isl_pw_qpolynomial *pwqp)
1002 isl_ctx *ctx;
1003 isl_set *dom;
1004 isl_space *map_space;
1005 isl_space *pwqp_space;
1006 unsigned n_in;
1007 int ok;
1009 ctx = isl_map_get_ctx(map);
1010 if (!ctx)
1011 goto error;
1013 map_space = isl_map_get_space(map);
1014 pwqp_space = isl_pw_qpolynomial_get_space(pwqp);
1015 ok = join_compatible(map_space, pwqp_space);
1016 isl_space_free(map_space);
1017 isl_space_free(pwqp_space);
1018 if (!ok)
1019 isl_die(ctx, isl_error_invalid, "incompatible dimensions",
1020 goto error);
1022 n_in = isl_map_dim(map, isl_dim_in);
1023 pwqp = isl_pw_qpolynomial_insert_dims(pwqp, isl_dim_in, 0, n_in);
1025 dom = isl_map_wrap(map);
1026 pwqp = isl_pw_qpolynomial_reset_domain_space(pwqp,
1027 isl_set_get_space(dom));
1029 pwqp = isl_pw_qpolynomial_intersect_domain(pwqp, dom);
1030 pwqp = isl_pw_qpolynomial_sum(pwqp);
1032 return pwqp;
1033 error:
1034 isl_map_free(map);
1035 isl_pw_qpolynomial_free(pwqp);
1036 return NULL;
1039 __isl_give isl_pw_qpolynomial *isl_set_apply_pw_qpolynomial(
1040 __isl_take isl_set *set, __isl_take isl_pw_qpolynomial *pwqp)
1042 isl_map *map;
1044 map = isl_map_from_range(set);
1045 pwqp = isl_map_apply_pw_qpolynomial(map, pwqp);
1046 pwqp = isl_pw_qpolynomial_project_domain_on_params(pwqp);
1047 return pwqp;
1050 struct barvinok_apply_data {
1051 isl_union_pw_qpolynomial *upwqp;
1052 isl_union_pw_qpolynomial *res;
1053 isl_map *map;
1056 static isl_stat pw_qpolynomial_apply(__isl_take isl_pw_qpolynomial *pwqp,
1057 void *user)
1059 isl_space *map_space;
1060 isl_space *pwqp_space;
1061 struct barvinok_apply_data *data = user;
1062 int ok;
1064 map_space = isl_map_get_space(data->map);
1065 pwqp_space = isl_pw_qpolynomial_get_space(pwqp);
1066 ok = join_compatible(map_space, pwqp_space);
1067 isl_space_free(map_space);
1068 isl_space_free(pwqp_space);
1070 if (ok) {
1071 isl_union_pw_qpolynomial *upwqp;
1073 pwqp = isl_map_apply_pw_qpolynomial(isl_map_copy(data->map),
1074 pwqp);
1075 upwqp = isl_union_pw_qpolynomial_from_pw_qpolynomial(pwqp);
1076 data->res = isl_union_pw_qpolynomial_add(data->res, upwqp);
1077 } else
1078 isl_pw_qpolynomial_free(pwqp);
1080 return isl_stat_ok;
1083 static isl_stat map_apply(__isl_take isl_map *map, void *user)
1085 struct barvinok_apply_data *data = user;
1086 isl_stat r;
1088 data->map = map;
1089 r = isl_union_pw_qpolynomial_foreach_pw_qpolynomial(data->upwqp,
1090 &pw_qpolynomial_apply, data);
1092 isl_map_free(map);
1093 return r;
1096 __isl_give isl_union_pw_qpolynomial *isl_union_map_apply_union_pw_qpolynomial(
1097 __isl_take isl_union_map *umap,
1098 __isl_take isl_union_pw_qpolynomial *upwqp)
1100 isl_space *space;
1101 struct barvinok_apply_data data;
1103 upwqp = isl_union_pw_qpolynomial_align_params(upwqp,
1104 isl_union_map_get_space(umap));
1105 umap = isl_union_map_align_params(umap,
1106 isl_union_pw_qpolynomial_get_space(upwqp));
1108 data.upwqp = upwqp;
1109 space = isl_union_pw_qpolynomial_get_space(upwqp);
1110 data.res = isl_union_pw_qpolynomial_zero(space);
1111 if (isl_union_map_foreach_map(umap, &map_apply, &data) < 0)
1112 goto error;
1114 isl_union_map_free(umap);
1115 isl_union_pw_qpolynomial_free(upwqp);
1117 return data.res;
1118 error:
1119 isl_union_map_free(umap);
1120 isl_union_pw_qpolynomial_free(upwqp);
1121 isl_union_pw_qpolynomial_free(data.res);
1122 return NULL;
1125 struct barvinok_apply_set_data {
1126 isl_union_pw_qpolynomial *upwqp;
1127 isl_union_pw_qpolynomial *res;
1128 isl_set *set;
1131 static isl_stat pw_qpolynomial_apply_set(__isl_take isl_pw_qpolynomial *pwqp,
1132 void *user)
1134 isl_space *set_space;
1135 isl_space *pwqp_space;
1136 struct barvinok_apply_set_data *data = user;
1137 int ok;
1139 set_space = isl_set_get_space(data->set);
1140 pwqp_space = isl_pw_qpolynomial_get_space(pwqp);
1141 ok = join_compatible(set_space, pwqp_space);
1142 isl_space_free(set_space);
1143 isl_space_free(pwqp_space);
1145 if (ok) {
1146 isl_union_pw_qpolynomial *upwqp;
1148 pwqp = isl_set_apply_pw_qpolynomial(isl_set_copy(data->set),
1149 pwqp);
1150 upwqp = isl_union_pw_qpolynomial_from_pw_qpolynomial(pwqp);
1151 data->res = isl_union_pw_qpolynomial_add(data->res, upwqp);
1152 } else
1153 isl_pw_qpolynomial_free(pwqp);
1155 return isl_stat_ok;
1158 static isl_stat set_apply(__isl_take isl_set *set, void *user)
1160 struct barvinok_apply_set_data *data = user;
1161 isl_stat r;
1163 data->set = set;
1164 r = isl_union_pw_qpolynomial_foreach_pw_qpolynomial(data->upwqp,
1165 &pw_qpolynomial_apply_set, data);
1167 isl_set_free(set);
1168 return r;
1171 __isl_give isl_union_pw_qpolynomial *isl_union_set_apply_union_pw_qpolynomial(
1172 __isl_take isl_union_set *uset,
1173 __isl_take isl_union_pw_qpolynomial *upwqp)
1175 isl_space *space;
1176 struct barvinok_apply_set_data data;
1178 upwqp = isl_union_pw_qpolynomial_align_params(upwqp,
1179 isl_union_set_get_space(uset));
1180 uset = isl_union_set_align_params(uset,
1181 isl_union_pw_qpolynomial_get_space(upwqp));
1183 data.upwqp = upwqp;
1184 space = isl_union_pw_qpolynomial_get_space(upwqp);
1185 data.res = isl_union_pw_qpolynomial_zero(space);
1186 if (isl_union_set_foreach_set(uset, &set_apply, &data) < 0)
1187 goto error;
1189 isl_union_set_free(uset);
1190 isl_union_pw_qpolynomial_free(upwqp);
1192 return data.res;
1193 error:
1194 isl_union_set_free(uset);
1195 isl_union_pw_qpolynomial_free(upwqp);
1196 isl_union_pw_qpolynomial_free(data.res);
1197 return NULL;
1200 evalue *evalue_sum(evalue *E, int nvar, unsigned MaxRays)
1202 evalue *sum;
1203 struct barvinok_options *options = barvinok_options_new_with_defaults();
1204 options->MaxRays = MaxRays;
1205 sum = barvinok_summate(E, nvar, options);
1206 barvinok_options_free(options);
1207 return sum;
1210 evalue *esum(evalue *e, int nvar)
1212 evalue *sum;
1213 struct barvinok_options *options = barvinok_options_new_with_defaults();
1214 sum = barvinok_summate(e, nvar, options);
1215 barvinok_options_free(options);
1216 return sum;
1219 /* Turn unweighted counting problem into "weighted" counting problem
1220 * with weight equal to 1 and call barvinok_summate on this weighted problem.
1222 evalue *barvinok_summate_unweighted(Polyhedron *P, Polyhedron *C,
1223 struct barvinok_options *options)
1225 Polyhedron *CA, *D;
1226 evalue e;
1227 evalue *sum;
1229 if (emptyQ(P) || emptyQ(C))
1230 return evalue_zero();
1232 CA = align_context(C, P->Dimension, options->MaxRays);
1233 D = DomainIntersection(P, CA, options->MaxRays);
1234 Domain_Free(CA);
1236 if (emptyQ(D)) {
1237 Domain_Free(D);
1238 return evalue_zero();
1241 value_init(e.d);
1242 e.x.p = new_enode(partition, 2, P->Dimension);
1243 EVALUE_SET_DOMAIN(e.x.p->arr[0], D);
1244 evalue_set_si(&e.x.p->arr[1], 1, 1);
1245 sum = barvinok_summate(&e, P->Dimension - C->Dimension, options);
1246 free_evalue_refs(&e);
1247 return sum;