From 31e2f7163410497feec8ab4f06a65f65c3b7e36c Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Sat, 26 Apr 2008 21:08:25 +0200 Subject: [PATCH] summate.c: handle equalities for all summation algorithms Before, we would only to it for the nested sums. --- bernoulli.c | 84 ++------------------------------------------------------ summate.c | 91 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- summate.h | 4 +++ 3 files changed, 94 insertions(+), 85 deletions(-) diff --git a/bernoulli.c b/bernoulli.c index a1835d8..c47c40a 100644 --- a/bernoulli.c +++ b/bernoulli.c @@ -5,6 +5,7 @@ #include "bernoulli.h" #include "lattice_point.h" #include "section_array.h" +#include "summate.h" #define ALLOC(type) (type*)malloc(sizeof(type)) #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type)) @@ -818,84 +819,6 @@ static evalue *sum_over_polytope_base(Polyhedron *P, evalue *E, unsigned nvar, return res; } -static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E, - unsigned nvar, - struct evalue_section_array *sections, - struct barvinok_options *options) -{ - 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(); - - assert(P->NbEq > 0); - - remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays); - - if (emptyQ(P)) { - Polyhedron_Free(P); - return evalue_zero(); - } - - new_nparam = CP ? CP->NbColumns-1 : dim - nvar; - 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], 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); - - if (new_dim-new_nparam > 0) { - sum = bernoulli_summate(P, E, new_dim-new_nparam, sections, options); - evalue_free(E); - Polyhedron_Free(P); - } else { - sum = ALLOC(evalue); - value_init(sum->d); - sum->x.p = new_enode(partition, 2, new_dim); - EVALUE_SET_DOMAIN(sum->x.p->arr[0], P); - value_clear(sum->x.p->arr[1].d); - sum->x.p->arr[1] = *E; - free(E); - } - - if (CP) { - evalue_backsubstitute(sum, CP, options->MaxRays); - Matrix_Free(CP); - } - - if (T) - Matrix_Free(T); - - return sum; -} - /* Splinter P into period slices along the first variable x = period y + c, * 0 <= c < perdiod, * ensuring no fractional part contains the first variable, * and sum over all slices. @@ -937,7 +860,7 @@ static evalue *sum_over_polytope_slices(Polyhedron *P, evalue *E, reduce_evalue(e); if (S->NbEq) - tmp = sum_over_polytope_with_equalities(S, e, nvar, sections, options); + tmp = barvinok_sum_over_polytope(S, e, nvar, sections, options); else tmp = sum_over_polytope_base(S, e, nvar, sections, options); @@ -973,9 +896,6 @@ evalue *bernoulli_summate(Polyhedron *P, evalue *E, unsigned nvar, evalue *sum; int exact = options->approximation_method == BV_APPROX_NONE; - if (P->NbEq) - return sum_over_polytope_with_equalities(P, E, nvar, sections, options); - value_init(period); move_best_to_front(&P, &E, nvar, exact ? &period : NULL); diff --git a/summate.c b/summate.c index 881578a..13746b9 100644 --- a/summate.c +++ b/summate.c @@ -5,10 +5,95 @@ #include "summate.h" #include "section_array.h" -static evalue *sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar, +#define ALLOC(type) (type*)malloc(sizeof(type)) +#define ALLOCN(type,n) (type*)malloc((n) * sizeof(type)) + +static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E, + unsigned nvar, struct evalue_section_array *sections, struct barvinok_options *options) { + 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(); + + assert(P->NbEq > 0); + + remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays); + + if (emptyQ(P)) { + Polyhedron_Free(P); + return evalue_zero(); + } + + new_nparam = CP ? CP->NbColumns-1 : dim - nvar; + 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], 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); + + if (new_dim-new_nparam > 0) { + sum = barvinok_sum_over_polytope(P, E, new_dim-new_nparam, + sections, options); + evalue_free(E); + Polyhedron_Free(P); + } else { + sum = ALLOC(evalue); + value_init(sum->d); + sum->x.p = new_enode(partition, 2, new_dim); + EVALUE_SET_DOMAIN(sum->x.p->arr[0], P); + value_clear(sum->x.p->arr[1].d); + sum->x.p->arr[1] = *E; + free(E); + } + + if (CP) { + evalue_backsubstitute(sum, CP, options->MaxRays); + Matrix_Free(CP); + } + + if (T) + Matrix_Free(T); + + return sum; +} + +evalue *barvinok_sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar, + struct evalue_section_array *sections, + struct barvinok_options *options) +{ + if (P->NbEq) + return sum_over_polytope_with_equalities(P, E, nvar, sections, options); + if (options->summation == BV_SUM_EULER) return euler_summate(P, E, nvar, options); else if (options->summation == BV_SUM_LAURENT) @@ -42,8 +127,8 @@ evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options) evalue *tmp; D->next = NULL; - tmp = sum_over_polytope(D, &e->x.p->arr[2*i+1], nvar, - §ions, options); + tmp = barvinok_sum_over_polytope(D, &e->x.p->arr[2*i+1], nvar, + §ions, options); assert(tmp); eadd(tmp, sum); evalue_free(tmp); diff --git a/summate.h b/summate.h index f4d983f..dd08ee8 100644 --- a/summate.h +++ b/summate.h @@ -1,4 +1,5 @@ #include +#include "section_array.h" struct barvinok_options; @@ -9,6 +10,9 @@ extern "C" { evalue *box_summate(Polyhedron *P, evalue *E, unsigned nvar, unsigned MaxRays); evalue *barvinok_summate_unweighted(Polyhedron *P, Polyhedron *C, struct barvinok_options *options); +evalue *barvinok_sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar, + struct evalue_section_array *sections, + struct barvinok_options *options); #if defined(__cplusplus) } -- 2.11.4.GIT