From 09a3229ac8c003c442fe2d42d7185c0f00f76f2e Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Sun, 4 May 2008 00:32:15 +0200 Subject: [PATCH] summate.c: barvinok_summate: handle products of polytopes for euler and laurent --- Makefile.am | 6 +- summate.c | 489 +++++++++++++++++++++++++++++++++++++++++++++--- tests/euler/EML-product | 11 ++ tests/evalue/product | 9 + 4 files changed, 481 insertions(+), 34 deletions(-) create mode 100644 tests/euler/EML-product create mode 100644 tests/evalue/product diff --git a/Makefile.am b/Makefile.am index 6804f7a..81357e0 100644 --- a/Makefile.am +++ b/Makefile.am @@ -381,10 +381,10 @@ check-evalue: @bv_barvinok_bound@ barvinok_summate$(EXEEXT) if test -f $$i; then \ echo $$i; \ if test -n "@bv_barvinok_bound@"; then \ - ./barvinok_bound$(EXEEXT) -T < $$i || exit; \ - echo $$i | ./test_bound$(EXEEXT) -q || exit; \ + ./barvinok_bound$(EXEEXT) -T -r30 < $$i || exit; \ + echo $$i | ./test_bound$(EXEEXT) -q -r30 || exit; \ fi; \ - ./barvinok_summate$(EXEEXT) -T < $$i || exit; \ + ./barvinok_summate$(EXEEXT) -T -r30 < $$i || exit; \ fi \ done check-euler: barvinok_summate$(EXEEXT) diff --git a/summate.c b/summate.c index db47b57..552f722 100644 --- a/summate.c +++ b/summate.c @@ -1,4 +1,5 @@ #include +#include #include "bernoulli.h" #include "euler.h" #include "laurent.h" @@ -12,6 +13,45 @@ extern void evalue_drop_floor(evalue *e, const evalue *floor); #define ALLOC(type) (type*)malloc(sizeof(type)) #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type)) +/* Apply the variable transformation specified by T and CP on + * the polynomial e. T expresses the old variables in terms + * of the new variables (and optionally also the new parameters), + * while CP expresses the old parameters in terms of the new + * parameters. + */ +static void transform_polynomial(evalue *E, Matrix *T, Matrix *CP, + unsigned nvar, unsigned nparam, + unsigned new_nvar, unsigned new_nparam) +{ + int j; + evalue **subs; + + subs = ALLOCN(evalue *, nvar+nparam); + + for (j = 0; j < nvar; ++j) { + if (T) + subs[j] = affine2evalue(T->p[j], T->p[T->NbRows-1][T->NbColumns-1], + T->NbColumns-1); + else + subs[j] = evalue_var(j); + } + for (j = 0; j < nparam; ++j) { + if (CP) + subs[nvar+j] = affine2evalue(CP->p[j], CP->p[nparam][new_nparam], + new_nparam); + else + subs[nvar+j] = evalue_var(j); + evalue_shift_variables(subs[nvar+j], 0, new_nvar); + } + + evalue_substitute(E, subs); + reduce_evalue(E); + + for (j = 0; j < nvar+nparam; ++j) + evalue_free(subs[j]); + free(subs); +} + static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E, unsigned nvar, struct evalue_section_array *sections, @@ -20,9 +60,7 @@ static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E, unsigned dim = P->Dimension; unsigned new_dim, new_nparam; Matrix *T = NULL, *CP = NULL; - evalue **subs; evalue *sum; - int j; if (emptyQ(P)) return evalue_zero(); @@ -40,30 +78,9 @@ static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E, new_dim = T ? T->NbColumns-1 : nvar + new_nparam; /* We can avoid these substitutions if E is a constant */ - subs = ALLOCN(evalue *, dim); - for (j = 0; j < nvar; ++j) { - if (T) - subs[j] = affine2evalue(T->p[j], T->p[nvar+new_nparam][new_dim], - new_dim); - else - subs[j] = evalue_var(j); - } - for (j = 0; j < dim-nvar; ++j) { - if (CP) - subs[nvar+j] = affine2evalue(CP->p[j], CP->p[dim-nvar][new_nparam], - new_nparam); - else - subs[nvar+j] = evalue_var(j); - evalue_shift_variables(subs[nvar+j], 0, new_dim-new_nparam); - } - E = evalue_dup(E); - evalue_substitute(E, subs); - reduce_evalue(E); - - for (j = 0; j < dim; ++j) - evalue_free(subs[j]); - free(subs); + transform_polynomial(E, T, CP, nvar, dim-nvar, + new_dim-new_nparam, new_nparam); if (new_dim-new_nparam > 0) { sum = barvinok_sum_over_polytope(P, E, new_dim-new_nparam, @@ -91,6 +108,410 @@ static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E, return sum; } +static evalue *sum_base(Polyhedron *P, evalue *E, unsigned nvar, + struct barvinok_options *options) +{ + if (options->summation == BV_SUM_EULER) + return euler_summate(P, E, nvar, options); + else if (options->summation == BV_SUM_LAURENT) + return laurent_summate(P, E, nvar, options); + assert(0); +} + +/* Count the number of non-zero terms in e when viewed as a polynomial + * in only the first nvar variables. "count" is the number counted + * so far. + */ +static int evalue_count_terms(const evalue *e, unsigned nvar, int count) +{ + int i; + + if (EVALUE_IS_ZERO(*e)) + return count; + + if (value_zero_p(e->d)) + assert(e->x.p->type == polynomial); + if (value_notzero_p(e->d) || e->x.p->pos >= nvar+1) + return count+1; + + for (i = 0; i < e->x.p->size; ++i) + count = evalue_count_terms(&e->x.p->arr[i], nvar, count); + + return count; +} + +/* Create placeholder structure for unzipping. + * A "polynomial" is created with size terms in variable pos, + * with each term having itself as coefficient. + */ +static evalue *create_placeholder(int size, int pos) +{ + int i, j; + evalue *E = ALLOC(evalue); + value_init(E->d); + E->x.p = new_enode(polynomial, size, pos+1); + for (i = 0; i < size; ++i) { + E->x.p->arr[i].x.p = new_enode(polynomial, i+1, pos+1); + for (j = 0; j < i; ++j) + evalue_set_si(&E->x.p->arr[i].x.p->arr[j], 0, 1); + evalue_set_si(&E->x.p->arr[i].x.p->arr[i], 1, 1); + } + return E; +} + +/* Interchange each non-zero term in e (when viewed as a polynomial + * in only the first nvar variables) with a placeholder in ph (created + * by create_placeholder), resulting in two polynomials in the + * placeholder variable such that for each non-zero term in e + * there is a power of the placeholder variable such that the factors + * in the first nvar variables form the coefficient of that power in + * the first polynomial (e) and the factors in the remaining variables + * form the coefficient of that power in the second polynomial (ph). + */ +static int evalue_unzip_terms(evalue *e, evalue *ph, unsigned nvar, int count) +{ + int i; + + if (EVALUE_IS_ZERO(*e)) + return count; + + if (value_zero_p(e->d)) + assert(e->x.p->type == polynomial); + if (value_notzero_p(e->d) || e->x.p->pos >= nvar+1) { + evalue t = *e; + *e = ph->x.p->arr[count]; + ph->x.p->arr[count] = t; + return count+1; + } + + for (i = 0; i < e->x.p->size; ++i) + count = evalue_unzip_terms(&e->x.p->arr[i], ph, nvar, count); + + return count; +} + +/* Remove n variables at pos (0-based) from the polyhedron P. + * Each of these variables is assumed to be completely free, + * i.e., there is a line in the polyhedron corresponding to + * each of these variables. + */ +static Polyhedron *Polyhedron_Remove_Columns(Polyhedron *P, unsigned pos, + unsigned n) +{ + int i, j; + unsigned NbConstraints = 0; + unsigned NbRays = 0; + Polyhedron *Q; + + if (!P) + return; + if (n == 0) + return; + + assert(pos <= P->Dimension); + + if (POL_HAS(P, POL_INEQUALITIES)) + NbConstraints = P->NbConstraints; + if (POL_HAS(P, POL_POINTS)) + NbRays = P->NbRays - n; + + Q = Polyhedron_Alloc(P->Dimension - n, NbConstraints, NbRays); + if (POL_HAS(P, POL_INEQUALITIES)) { + Q->NbEq = P->NbEq; + for (i = 0; i < P->NbConstraints; ++i) { + Vector_Copy(P->Constraint[i], Q->Constraint[i], 1+pos); + Vector_Copy(P->Constraint[i]+1+pos+n, Q->Constraint[i]+1+pos, + Q->Dimension-pos+1); + } + } + if (POL_HAS(P, POL_POINTS)) { + Q->NbBid = P->NbBid - n; + for (i = 0; i < n; ++i) + value_set_si(Q->Ray[i][1+pos+i], 1); + for (i = 0, j = 0; i < P->NbRays; ++i) { + int line = First_Non_Zero(P->Ray[i], 1+P->Dimension+1); + assert(line != -1); + if (line-1 >= pos && line-1 < pos+n) { + ++j; + assert(j <= n); + continue; + } + assert(i-j < Q->NbRays); + Vector_Copy(P->Ray[i], Q->Ray[i-j], 1+pos); + Vector_Copy(P->Ray[i]+1+pos+n, Q->Ray[i-j]+1+pos, + Q->Dimension-pos+1); + } + } + POL_SET(Q, POL_VALID); + if (POL_HAS(P, POL_INEQUALITIES)) + POL_SET(Q, POL_INEQUALITIES); + if (POL_HAS(P, POL_POINTS)) + POL_SET(Q, POL_POINTS); + if (POL_HAS(P, POL_VERTICES)) + POL_SET(Q, POL_VERTICES); + return Q; +} + +/* Remove n variables at pos (0-based) from the union of polyhedra P. + * Each of these variables is assumed to be completely free, + * i.e., there is a line in the polyhedron corresponding to + * each of these variables. + */ +static Polyhedron *Domain_Remove_Columns(Polyhedron *P, unsigned pos, + unsigned n) +{ + Polyhedron *R; + Polyhedron **next = &R; + + for (; P; P = P->next) { + *next = Polyhedron_Remove_Columns(P, pos, n); + next = &(*next)->next; + } + return R; +} + +/* Drop n parameters starting at first from partition evalue e */ +static void drop_parameters(evalue *e, int first, int n) +{ + int i; + + if (EVALUE_IS_ZERO(*e)) + return; + + assert(value_zero_p(e->d) && e->x.p->type == partition); + for (i = 0; i < e->x.p->size/2; ++i) { + Polyhedron *P = EVALUE_DOMAIN(e->x.p->arr[2*i]); + Polyhedron *Q = Domain_Remove_Columns(P, first, n); + EVALUE_SET_DOMAIN(e->x.p->arr[2*i], Q); + Domain_Free(P); + evalue_shift_variables(&e->x.p->arr[2*i+1], first, -n); + } + e->x.p->pos -= n; +} + +static void extract_term_into(const evalue *src, int var, int exp, evalue *dst) +{ + int i; + + if (value_notzero_p(src->d) || + src->x.p->type != polynomial || + src->x.p->pos > var+1) { + if (exp == 0) + evalue_copy(dst, src); + else + evalue_set_si(dst, 0, 1); + return; + } + + if (src->x.p->pos == var+1) { + if (src->x.p->size > exp) + evalue_copy(dst, &src->x.p->arr[exp]); + else + evalue_set_si(dst, 0, 1); + return; + } + + dst->x.p = new_enode(polynomial, src->x.p->size, src->x.p->pos); + for (i = 0; i < src->x.p->size; ++i) + extract_term_into(&src->x.p->arr[i], var, exp, + &dst->x.p->arr[i]); +} + +/* Extract the coefficient of var^exp. + */ +static evalue *extract_term(const evalue *e, int var, int exp) +{ + int i; + evalue *res; + + if (EVALUE_IS_ZERO(*e)) + return evalue_zero(); + + assert(value_zero_p(e->d) && e->x.p->type == partition); + res = ALLOC(evalue); + value_init(res->d); + res->x.p = new_enode(partition, e->x.p->size, e->x.p->pos); + for (i = 0; i < e->x.p->size/2; ++i) { + EVALUE_SET_DOMAIN(res->x.p->arr[2*i], + Domain_Copy(EVALUE_DOMAIN(e->x.p->arr[2*i]))); + extract_term_into(&e->x.p->arr[2*i+1], var, exp, + &res->x.p->arr[2*i+1]); + reduce_evalue(&res->x.p->arr[2*i+1]); + } + return res; +} + +/* Insert n free variables at pos (0-based) in the polyhedron P. + */ +static Polyhedron *Polyhedron_Insert_Columns(Polyhedron *P, unsigned pos, + unsigned n) +{ + int i; + unsigned NbConstraints = 0; + unsigned NbRays = 0; + Polyhedron *Q; + + if (!P) + return; + if (n == 0) + return; + + assert(pos <= P->Dimension); + + if (POL_HAS(P, POL_INEQUALITIES)) + NbConstraints = P->NbConstraints; + if (POL_HAS(P, POL_POINTS)) + NbRays = P->NbRays + n; + + Q = Polyhedron_Alloc(P->Dimension+n, NbConstraints, NbRays); + if (POL_HAS(P, POL_INEQUALITIES)) { + Q->NbEq = P->NbEq; + for (i = 0; i < P->NbConstraints; ++i) { + Vector_Copy(P->Constraint[i], Q->Constraint[i], 1+pos); + Vector_Copy(P->Constraint[i]+1+pos, Q->Constraint[i]+1+pos+n, + P->Dimension-pos+1); + } + } + if (POL_HAS(P, POL_POINTS)) { + Q->NbBid = P->NbBid + n; + for (i = 0; i < n; ++i) + value_set_si(Q->Ray[i][1+pos+i], 1); + for (i = 0; i < P->NbRays; ++i) { + Vector_Copy(P->Constraint[i], Q->Constraint[n+i], 1+pos); + Vector_Copy(P->Constraint[i]+1+pos, Q->Constraint[n+i]+1+pos+n, + P->Dimension-pos+1); + } + } + POL_SET(Q, POL_VALID); + if (POL_HAS(P, POL_INEQUALITIES)) + POL_SET(Q, POL_INEQUALITIES); + if (POL_HAS(P, POL_POINTS)) + POL_SET(Q, POL_POINTS); + if (POL_HAS(P, POL_VERTICES)) + POL_SET(Q, POL_VERTICES); + return Q; +} + +/* Perform summation of e over a list of 1 or more factors F, with context C. + * nvar is the total number of variables in the remaining factors. + * extra is the number of placeholder parameters introduced in e, + * but not (yet) in F or C. + * + * If there is only one factor left, F is intersected with the + * context C, the placeholder variables are added, and then + * e is summed over the resulting parametric polytope. + * + * If there is more than one factor left, we create to polynomials + * in a new placeholder variable (which is placed after the regular + * parameters, but before any previously introduced placeholder + * variables) that has the factors of the variables in the first + * factor of F and the factor of the remaining variables of + * each term as its coefficients. + * These two polynomials are then summed over their domains + * and afterwards the results are combined and the placeholder + * variable is removed again. + */ +static evalue *sum_factors(Polyhedron *F, Polyhedron *C, evalue *e, + unsigned nvar, unsigned extra, + struct barvinok_options *options) +{ + Polyhedron *P = F; + unsigned nparam = C->Dimension; + unsigned F_var = F->Dimension - C->Dimension; + int i, n; + evalue *s1, *s2, *s; + evalue *ph; + + if (!F->next) { + Polyhedron *CA = align_context(C, nvar+nparam, options->MaxRays); + Polyhedron *P = DomainIntersection(F, CA, options->MaxRays); + Polyhedron *Q = Polyhedron_Insert_Columns(P, nvar+nparam, extra); + Polyhedron_Free(CA); + Polyhedron_Free(F); + Polyhedron_Free(P); + evalue *sum = sum_base(Q, e, nvar, options); + Polyhedron_Free(Q); + return sum; + } + + n = evalue_count_terms(e, F_var, 0); + ph = create_placeholder(n, nvar+nparam); + evalue_shift_variables(e, nvar+nparam, 1); + evalue_unzip_terms(e, ph, F_var, 0); + evalue_shift_variables(e, nvar, -(nvar-F_var)); + evalue_reorder_terms(ph); + evalue_shift_variables(ph, 0, -F_var); + + s2 = sum_factors(F->next, C, ph, nvar-F_var, extra+1, options); + evalue_free(ph); + F->next = NULL; + s1 = sum_factors(F, C, e, F_var, extra+1, options); + + if (n == 1) { + /* remove placeholder "polynomial" */ + reduce_evalue(s1); + emul(s1, s2); + evalue_free(s1); + drop_parameters(s2, nparam, 1); + return s2; + } + + s = evalue_zero(); + for (i = 0; i < n; ++i) { + evalue *t1, *t2; + t1 = extract_term(s1, nparam, i); + t2 = extract_term(s2, nparam, i); + emul(t1, t2); + eadd(t2, s); + evalue_free(t1); + evalue_free(t2); + } + evalue_free(s1); + evalue_free(s2); + + drop_parameters(s, nparam, 1); + return s; +} + +/* Perform summation over a product of factors F, obtained using + * variable transformation T from the original problem specification. + * + * We first perform the corresponding transformation on the polynomial E, + * compute the common context over all factors and then perform + * the actual summation over the factors. + */ +static evalue *sum_product(Polyhedron *F, evalue *E, Matrix *T, unsigned nparam, + struct barvinok_options *options) +{ + int i; + Matrix *T2; + unsigned nvar = T->NbRows; + Polyhedron *C; + evalue *sum; + + assert(nvar == T->NbColumns); + T2 = Matrix_Alloc(nvar+1, nvar+1); + for (i = 0; i < nvar; ++i) + Vector_Copy(T->p[i], T2->p[i], nvar); + value_set_si(T2->p[nvar][nvar], 1); + + transform_polynomial(E, T2, NULL, nvar, nparam, nvar, nparam); + + C = Factor_Context(F, nparam, options->MaxRays); + if (F->Dimension == nparam) { + Polyhedron *T = F; + F = F->next; + T->next = NULL; + Polyhedron_Free(T); + } + sum = sum_factors(F, C, E, nvar, 0, options); + + Polyhedron_Free(C); + Matrix_Free(T2); + Matrix_Free(T); + return sum; +} + /* Add two constraints corresponding to floor = floor(e/d), * * e - d t >= 0 @@ -184,12 +605,18 @@ static evalue *sum_step_polynomial(Polyhedron *P, evalue *E, unsigned nvar, if (EVALUE_IS_ZERO(*cur)) t = NULL; - else if (options->summation == BV_SUM_EULER) - t = euler_summate(P, cur, nvar, options); - else if (options->summation == BV_SUM_LAURENT) - t = laurent_summate(P, cur, nvar, options); - else - assert(0); + else { + Matrix *T; + unsigned nparam = P->Dimension - nvar; + Polyhedron *F = Polyhedron_Factor(P, nparam, &T, options->MaxRays); + if (!F) + t = sum_base(P, cur, nvar, options); + else { + if (cur == E) + cur = evalue_dup(cur); + t = sum_product(F, cur, T, nparam, options); + } + } if (E != cur) evalue_free(cur); diff --git a/tests/euler/EML-product b/tests/euler/EML-product new file mode 100644 index 0000000..7c46fd3 --- /dev/null +++ b/tests/euler/EML-product @@ -0,0 +1,11 @@ +#variables x,y,x2,y2 + +3 x + 210y -43 >=0 +6795 x -2660y -1733 >= 0 +-1365 x +112y +5899 >= 0 +3 x2 + 210y2 -43 >=0 +6795 x2 -2660y2 -1733 >= 0 +-1365 x2 +112y2 +5899 >= 0 + +x^2 * y^2 * x2^2 * y2^2 * N + diff --git a/tests/evalue/product b/tests/evalue/product new file mode 100644 index 0000000..1f5d0c5 --- /dev/null +++ b/tests/evalue/product @@ -0,0 +1,9 @@ +#variables x,y,z + x >= 0 + -x + N >= 0 + y >= 0 + -y + N >= 0 + z >= 0 + -z + N >= 0 + +x^3*y^2 + N*y*z -- 2.11.4.GIT