From 01aa2e7a19d8c2bafb650b9772129f59caf4257e Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Mon, 19 Mar 2007 23:07:18 +0100 Subject: [PATCH] barvinok_enumerate: more accurate polynomial approximation These polynomial approximations of quasi-polynomials are described by Benoit Meister's and should be more accurate than his "faster" approximation that had already been implemented. --- Makefile.am | 2 + barvinok.cc | 17 +++-- barvinok/options.h | 3 +- barvinok/util.h | 1 + options.c | 6 +- scale.c | 202 +++++++++++++++++++++++++++++++++++++++++++++++++++++ scale.h | 12 ++++ util.c | 27 +++++-- 8 files changed, 254 insertions(+), 16 deletions(-) create mode 100644 scale.c create mode 100644 scale.h diff --git a/Makefile.am b/Makefile.am index 7cc7026..a2fa810 100644 --- a/Makefile.am +++ b/Makefile.am @@ -66,6 +66,8 @@ libbarvinok_la_SOURCES = \ reducer.cc \ reducer.h \ remove_equalities.h \ + scale.c \ + scale.h \ scarf.cc \ mat_util.cc \ mat_util.h \ diff --git a/barvinok.cc b/barvinok.cc index 6ec7122..ffe984b 100644 --- a/barvinok.cc +++ b/barvinok.cc @@ -23,6 +23,7 @@ extern "C" { #include "reduce_domain.h" #include "genfun_constructor.h" #include "remove_equalities.h" +#include "scale.h" #ifndef HAVE_PARAM_POLYHEDRON_SCALE_INTEGER extern "C" void Param_Polyhedron_Scale_Integer(Param_Polyhedron *PP, Polyhedron **P, @@ -1857,9 +1858,10 @@ static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C, barvinok_options *options) { unsigned nparam = C->Dimension; - bool pre_approx = options->approximation_method == BV_APPROX_SCALE; + bool scale_fast = options->approximation_method == BV_APPROX_SCALE_FAST; + bool scale = options->approximation_method == BV_APPROX_SCALE; - if (P->Dimension - nparam == 1 && !pre_approx) + if (P->Dimension - nparam == 1 && !scale_fast) return ParamLine_Length(P, C, options); Param_Polyhedron *PP = NULL; @@ -1872,7 +1874,7 @@ static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C, Value det; Polyhedron *T; - if (options->approximation_method == BV_APPROX_SCALE) { + if (scale || scale_fast) { if (options->polynomial_approximation == BV_APPROX_SIGN_UPPER) P = Polyhedron_Inflate(P, nparam, options->MaxRays); if (options->polynomial_approximation == BV_APPROX_SIGN_LOWER) { @@ -1912,10 +1914,13 @@ out: nparam = CT->NbRows - 1; } - if (pre_approx) { + if (scale || scale_fast) { value_init(det); Polyhedron *T = P; - Param_Polyhedron_Scale_Integer(PP, &T, &det, options->MaxRays); + if (scale_fast) + Param_Polyhedron_Scale_Integer(PP, &T, &det, options->MaxRays); + else + Param_Polyhedron_Scale_Integer_Slow(PP, &T, &det, options->MaxRays); if (P != Porig) Polyhedron_Free(P); P = T; @@ -1993,7 +1998,7 @@ try_again: delete [] s; delete [] fVD; - if (pre_approx) { + if (scale || scale_fast) { evalue_div(eres, det); value_clear(det); } diff --git a/barvinok/options.h b/barvinok/options.h index 0cda29e..39997b3 100644 --- a/barvinok/options.h +++ b/barvinok/options.h @@ -46,7 +46,8 @@ struct barvinok_options { int polynomial_approximation; #define BV_APPROX_NONE 0 #define BV_APPROX_DROP 1 - #define BV_APPROX_SCALE 2 + #define BV_APPROX_SCALE_FAST 2 + #define BV_APPROX_SCALE 3 int approximation_method; /* basis reduction options */ diff --git a/barvinok/util.h b/barvinok/util.h index 3e0833d..8ca31ad 100644 --- a/barvinok/util.h +++ b/barvinok/util.h @@ -39,6 +39,7 @@ Polyhedron *Polyhedron_Read(unsigned MaxRays); Polyhedron* Polyhedron_Polar(Polyhedron *P, unsigned NbMaxRays); void Polyhedron_Polarize(Polyhedron *P); Polyhedron* supporting_cone(Polyhedron *P, int v); +unsigned char *supporting_constraints(Polyhedron *P, Param_Vertices *v, int *n); Polyhedron* supporting_cone_p(Polyhedron *P, Param_Vertices *v); Polyhedron* triangulate_cone(Polyhedron *P, unsigned NbMaxCons); void check_triangulization(Polyhedron *P, Polyhedron *T); diff --git a/options.c b/options.c index 9013ccd..bae6ac4 100644 --- a/options.c +++ b/options.c @@ -99,7 +99,7 @@ struct argp_option barvinok_argp_options[] = { { "table", BV_OPT_TABLE, 0, 0 }, { "specialization", BV_OPT_SPECIALIZATION, "[bf|df|random]", 0 }, { "polynomial-approximation", BV_OPT_POLAPPROX, "lower|upper", 1 }, - { "approximation-method", BV_OPT_APPROX, "scale|drop", 0, + { "approximation-method", BV_OPT_APPROX, "scale|scale-fast|drop", 0, "method to use in polynomial approximation [default: drop]" }, { "gbr", BV_OPT_GBR, "[cdd]", 0, "solver to use for basis reduction" }, @@ -140,7 +140,7 @@ error_t barvinok_parse_opt(int key, char *arg, struct argp_state *state) if (!arg) { options->polynomial_approximation = BV_APPROX_SIGN_APPROX; if (options->approximation_method == BV_APPROX_NONE) - options->approximation_method = BV_APPROX_SCALE; + options->approximation_method = BV_APPROX_SCALE_FAST; } else { if (!strcmp(arg, "lower")) options->polynomial_approximation = BV_APPROX_SIGN_LOWER; @@ -153,6 +153,8 @@ error_t barvinok_parse_opt(int key, char *arg, struct argp_state *state) case BV_OPT_APPROX: if (!strcmp(arg, "scale")) options->approximation_method = BV_APPROX_SCALE; + else if (!strcmp(arg, "scale-fast")) + options->approximation_method = BV_APPROX_SCALE_FAST; else if (!strcmp(arg, "drop")) options->approximation_method = BV_APPROX_DROP; break; diff --git a/scale.c b/scale.c new file mode 100644 index 0000000..d1d8b62 --- /dev/null +++ b/scale.c @@ -0,0 +1,202 @@ +#include +#include "scale.h" + +static Lattice *extract_lattice(Matrix *M, unsigned nvar) +{ + int row, col; + Matrix *H, *Q, *U, *Li; + Lattice *L; + int ok; + + left_hermite(M, &H, &Q, &U); + Matrix_Free(U); + + Li = Matrix_Alloc(nvar+1, nvar+1); + L = Matrix_Alloc(nvar+1, nvar+1); + value_set_si(Li->p[nvar][nvar], 1); + + for (row = 0, col = 0; col < nvar; ++col) { + while (value_zero_p(H->p[row][col])) + ++row; + Vector_Copy(Q->p[row], Li->p[col], nvar); + } + Matrix_Free(H); + Matrix_Free(Q); + + ok = Matrix_Inverse(Li, L); + assert(ok); + Matrix_Free(Li); + + return L; +} + +/* Returns the smallest (wrt inclusion) lattice that contains both X and Y */ +static Lattice *LatticeJoin(Lattice *X, Lattice *Y) +{ + int i; + int dim = X->NbRows-1; + Value lcm; + Value tmp; + Lattice *L; + Matrix *M, *H, *U, *Q; + + assert(X->NbColumns-1 == dim); + assert(Y->NbRows-1 == dim); + assert(Y->NbColumns-1 == dim); + + value_init(lcm); + value_init(tmp); + + M = Matrix_Alloc(dim, 2*dim); + value_lcm(X->p[dim][dim], Y->p[dim][dim], &lcm); + + value_division(tmp, lcm, X->p[dim][dim]); + for (i = 0; i < dim; ++i) + Vector_Scale(X->p[i], M->p[i], tmp, dim); + value_division(tmp, lcm, Y->p[dim][dim]); + for (i = 0; i < dim; ++i) + Vector_Scale(Y->p[i], M->p[i]+dim, tmp, dim); + + left_hermite(M, &H, &Q, &U); + Matrix_Free(M); + Matrix_Free(Q); + Matrix_Free(U); + + L = Matrix_Alloc(dim+1, dim+1); + value_assign(L->p[dim][dim], lcm); + for (i = 0; i < dim; ++i) + Vector_Copy(H->p[i], L->p[i], dim); + Matrix_Free(H); + + value_clear(tmp); + value_clear(lcm); + return L; +} + +static void Param_Vertex_Common_Denominator(Param_Vertices *V) +{ + unsigned dim; + Value lcm; + int i; + + assert(V->Vertex->NbRows > 0); + dim = V->Vertex->NbColumns-2; + + value_init(lcm); + + value_assign(lcm, V->Vertex->p[0][dim+1]); + for (i = 1; i < V->Vertex->NbRows; ++i) + value_lcm(V->Vertex->p[i][dim+1], lcm, &lcm); + + for (i = 0; i < V->Vertex->NbRows; ++i) { + if (value_eq(V->Vertex->p[i][dim+1], lcm)) + continue; + value_division(V->Vertex->p[i][dim+1], lcm, V->Vertex->p[i][dim+1]); + Vector_Scale(V->Vertex->p[i], V->Vertex->p[i], + V->Vertex->p[i][dim+1], dim+1); + value_assign(V->Vertex->p[i][dim+1], lcm); + } + + value_clear(lcm); +} + +static void Param_Vertex_Image(Param_Vertices *V, Matrix *T) +{ + unsigned nvar = V->Vertex->NbRows; + unsigned nparam = V->Vertex->NbColumns - 2; + Matrix *Vertex; + int i; + + Param_Vertex_Common_Denominator(V); + Vertex = Matrix_Alloc(V->Vertex->NbRows, V->Vertex->NbColumns); + Matrix_Product(T, V->Vertex, Vertex); + for (i = 0; i < nvar; ++i) { + value_assign(Vertex->p[i][nparam+1], V->Vertex->p[i][nparam+1]); + Vector_Normalize(Vertex->p[i], nparam+2); + } + Matrix_Free(V->Vertex); + V->Vertex = Vertex; +} + +/* Scales the parametric polyhedron with constraints *P and vertices PP + * such that the number of integer points can be represented by a polynomial. + * Both *P and P->Vertex are adapted according to the scaling. + * The scaling factor is returned in *det. + * The enumerator of the scaled parametric polyhedron should be divided + * by this number to obtain an approximation of the enumerator of the + * original parametric polyhedron. + * + * The algorithm is described in "Approximating Ehrhart Polynomials using + * affine transformations" by B. Meister. + */ +void Param_Polyhedron_Scale_Integer_Slow(Param_Polyhedron *PP, Polyhedron **P, + Value *det, unsigned MaxRays) +{ + Param_Vertices *V; + unsigned dim = (*P)->Dimension; + unsigned nparam; + unsigned nvar; + Lattice *L = NULL, *Li; + Matrix *T; + Matrix *expansion; + int i; + int ok; + + if (!PP->nbV) + return; + + nparam = PP->V->Vertex->NbColumns - 2; + nvar = dim - nparam; + + for (V = PP->V; V; V = V->next) { + Lattice *L2; + Matrix *M; + int i, j, n; + unsigned char *supporting; + + supporting = supporting_constraints(*P, V, &n); + M = Matrix_Alloc(n, (*P)->Dimension); + for (i = 0, j = 0; i < (*P)->NbConstraints; ++i) + if (supporting[i]) + Vector_Copy((*P)->Constraint[i]+1, M->p[j++], (*P)->Dimension); + free(supporting); + L2 = extract_lattice(M, nvar); + Matrix_Free(M); + + if (!L) + L = L2; + else { + Lattice *L3 = LatticeJoin(L, L2); + Matrix_Free(L); + Matrix_Free(L2); + L = L3; + } + } + + /* apply the variable expansion to the polyhedron (constraints) */ + expansion = Matrix_Alloc(nvar + nparam + 1, nvar + nparam + 1); + for (i = 0; i < nvar; ++i) + Vector_Copy(L->p[i], expansion->p[i], nvar); + for (i = nvar; i < nvar+nparam+1; ++i) + value_assign(expansion->p[i][i], L->p[nvar][nvar]); + + *P = Polyhedron_Preimage(*P, expansion, MaxRays); + Matrix_Free(expansion); + + /* apply the variable expansion to the parametric vertices */ + Li = Matrix_Alloc(nvar+1, nvar+1); + ok = Matrix_Inverse(L, Li); + assert(ok); + Matrix_Free(L); + assert(value_one_p(Li->p[nvar][nvar])); + T = Matrix_Alloc(nvar, nvar); + value_set_si(*det, 1); + for (i = 0; i < nvar; ++i) { + value_multiply(*det, *det, Li->p[i][i]); + Vector_Copy(Li->p[i], T->p[i], nvar); + } + Matrix_Free(Li); + for (V = PP->V; V; V = V->next) + Param_Vertex_Image(V, T); + Matrix_Free(T); +} diff --git a/scale.h b/scale.h new file mode 100644 index 0000000..24027a7 --- /dev/null +++ b/scale.h @@ -0,0 +1,12 @@ +#include + +#if defined(__cplusplus) +extern "C" { +#endif + +void Param_Polyhedron_Scale_Integer_Slow(Param_Polyhedron *PP, Polyhedron **P, + Value *det, unsigned MaxRays); + +#if defined(__cplusplus) +} +#endif diff --git a/util.c b/util.c index 3e0a912..aee7439 100644 --- a/util.c +++ b/util.c @@ -174,15 +174,14 @@ void value_lcm(const Value i, const Value j, Value* lcm) value_clear(aux); } -Polyhedron* supporting_cone_p(Polyhedron *P, Param_Vertices *v) +unsigned char *supporting_constraints(Polyhedron *P, Param_Vertices *v, int *n) { - Matrix *M; Value lcm, tmp, tmp2; - unsigned char *supporting = (unsigned char *)malloc(P->NbConstraints); unsigned dim = P->Dimension + 2; unsigned nparam = v->Vertex->NbColumns - 2; unsigned nvar = dim - nparam - 2; - int i, n, j; + unsigned char *supporting = (unsigned char *)malloc(P->NbConstraints); + int i, j; Vector *row; assert(supporting); @@ -192,7 +191,7 @@ Polyhedron* supporting_cone_p(Polyhedron *P, Param_Vertices *v) value_init(tmp); value_init(tmp2); value_set_si(lcm, 1); - for (i = 0, n = 0; i < P->NbConstraints; ++i) { + for (i = 0, *n = 0; i < P->NbConstraints; ++i) { Vector_Set(row->p, 0, nparam+1); for (j = 0 ; j < nvar; ++j) { value_set_si(tmp, 1); @@ -213,13 +212,27 @@ Polyhedron* supporting_cone_p(Polyhedron *P, Param_Vertices *v) if (value_notzero_p(row->p[j])) break; if ((supporting[i] = (j == nparam + 1))) - ++n; + ++*n; } - assert(n >= nvar); + assert(*n >= nvar); value_clear(tmp); value_clear(tmp2); value_clear(lcm); Vector_Free(row); + + return supporting; +} + +Polyhedron* supporting_cone_p(Polyhedron *P, Param_Vertices *v) +{ + Matrix *M; + unsigned dim = P->Dimension + 2; + unsigned nparam = v->Vertex->NbColumns - 2; + unsigned nvar = dim - nparam - 2; + int i, n, j; + unsigned char *supporting; + + supporting = supporting_constraints(P, v, &n); M = Matrix_Alloc(n, nvar+2); assert(M); for (i = 0, j = 0; i < P->NbConstraints; ++i) -- 2.11.4.GIT