From c3550e872b04c72c0b2e889229cc184352c1da3c Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Tue, 16 Oct 2007 20:33:40 +0200 Subject: [PATCH] Replace incremental infinite set counter by "regular" infinite set counter The incremental infinite set counter was incorrect since it would only keep the negative terms in the final specialization step. It therefore only worked for 1-dimensional problems, where the final specialization step is the only specialization step, and so it wasn't actually incremental at all. --- counter.cc | 136 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ counter.h | 34 ++++++++++++++++ dpoly.cc | 1 + genfun.cc | 4 +- reducer.cc | 44 -------------------- reducer.h | 32 --------------- testlib.cc | 72 ++++++++++++++++++++++++++++++++ 7 files changed, 246 insertions(+), 77 deletions(-) diff --git a/counter.cc b/counter.cc index 3ad8463..87ede3d 100644 --- a/counter.cc +++ b/counter.cc @@ -154,3 +154,139 @@ void tcounter::add_lattice_points(int sign) else mpq_add(count, count, tcount); } + + +/* + * Set lambda to a random vector that has a positive inner product + * with all the rays of the context { x | A x + b >= 0 }. + * + * To do so, we take d rows A' from the constraint matrix A. + * For every ray, we have + * A' r >= 0 + * We compute a random positive row vector x' and set x = x' A'. + * We then have, for each ray r, + * x r = x' A' r >= 0 + * Although we can take any d rows from A, we choose linearly + * independent rows from A to avoid the elements of the transformed + * random vector to have simple linear relations, which would + * increase the risk of the vector being orthogonal to one of + * powers in the denominator of one of the terms in the generating + * function. + */ +void infinite_counter::init(Polyhedron *context) +{ + Matrix *M, *H, *Q, *U; + mat_ZZ A; + + randomvector(context, lambda, context->Dimension); + + M = Matrix_Alloc(context->NbConstraints, context->Dimension); + for (int i = 0; i < context->NbConstraints; ++i) + Vector_Copy(context->Constraint[i]+1, M->p[i], context->Dimension); + left_hermite(M, &H, &Q, &U); + Matrix_Free(Q); + Matrix_Free(U); + + for (int col = 0, row = 0; col < H->NbColumns; ++col, ++row) { + for (; row < H->NbRows; ++row) + if (value_notzero_p(H->p[row][col])) + break; + assert(row < H->NbRows); + Vector_Copy(M->p[row], M->p[col], M->NbColumns); + } + matrix2zz(M, A, context->Dimension, context->Dimension); + Matrix_Free(H); + Matrix_Free(M); + + for (int i = 0; i < lambda.length(); ++i) + lambda[i] = abs(lambda[i]); + lambda = lambda * A; +} + +static ZZ LCM(const ZZ& a, const ZZ& b) +{ + return a * b / GCD(a, b); +} + +/* Normalize the powers in the denominator to be positive + * and return -1 is the sign has to be changed. + */ +static int normalized_sign(vec_ZZ& num, vec_ZZ& den) +{ + int change = 0; + + for (int j = 0; j < den.length(); ++j) { + if (den[j] > 0) + change ^= 1; + else { + den[j] = abs(den[j]); + for (int k = 0; k < num.length(); ++k) + num[k] += den[j]; + } + } + return change ? -1 : 1; +} + +void infinite_counter::reduce(const vec_QQ& c, const mat_ZZ& num, + const mat_ZZ& den_f) +{ + mpq_t factor; + mpq_init(factor); + unsigned len = den_f.NumRows(); + + ZZ lcm = c[0].d; + for (int i = 1; i < c.length(); ++i) + lcm = LCM(lcm, c[i].d); + + vec_ZZ coeff; + coeff.SetLength(c.length()); + for (int i = 0; i < c.length(); ++i) + coeff[i] = c[i].n * lcm/c[i].d; + + if (len == 0) { + for (int i = 0; i < c.length(); ++i) { + zz2value(coeff[i], tz); + value_addto(mpq_numref(factor), mpq_numref(factor), tz); + } + zz2value(lcm, tz); + value_assign(mpq_denref(factor), tz); + mpq_add(count[0], count[0], factor); + mpq_clear(factor); + return; + } + + vec_ZZ num_s = num * lambda; + vec_ZZ den_s = den_f * lambda; + for (int i = 0; i < den_s.length(); ++i) + assert(den_s[i] != 0); + int sign = normalized_sign(num_s, den_s); + if (sign < 0) + coeff = -coeff; + + dpoly n(len); + zz2value(num_s[0], tz); + add_falling_powers(n, tz); + zz2value(coeff[0], tz); + n *= tz; + for (int i = 1; i < c.length(); ++i) { + dpoly t(len); + zz2value(num_s[i], tz); + add_falling_powers(t, tz); + zz2value(coeff[i], tz); + t *= tz; + n += t; + } + zz2value(den_s[0], tz); + dpoly d(len, tz, 1); + for (int i = 1; i < len; ++i) { + zz2value(den_s[i], tz); + dpoly fact(len, tz, 1); + d *= fact; + } + value_set_si(mpq_numref(factor), 1); + zz2value(lcm, tz); + value_assign(mpq_denref(factor), tz); + n.div(d, count, factor); + + mpq_clear(factor); +} diff --git a/counter.h b/counter.h index fe04530..8ff7698 100644 --- a/counter.h +++ b/counter.h @@ -76,3 +76,37 @@ struct tcounter : public counter_base { virtual void add_lattice_points(int sign); }; + +/* A counter for possibly infinite sets. + * Rather than just keeping track of the constant term + * of the Laurent expansions, we also keep track of the + * coefficients of negative powers. + * If any of these is non-zero, then the counted set is infinite. + */ +struct infinite_counter { + /* an array of coefficients; count[i] is the coeffient of + * the term with power -i. + */ + vec_ZZ lambda; + mpq_t *count; + unsigned maxlen; + Value tz; + + infinite_counter(unsigned dim, unsigned maxlen) : maxlen(maxlen) { + count = new mpq_t[maxlen+1]; + for (int i = 0; i <= maxlen; ++i) + mpq_init(count[i]); + value_init(tz); + } + + void init(Polyhedron *context); + + void reduce(const vec_QQ& c, const mat_ZZ& num, const mat_ZZ& den_f); + + ~infinite_counter() { + for (int i = 0; i <= maxlen; ++i) + mpq_clear(count[i]); + delete [] count; + value_clear(tz); + } +}; diff --git a/dpoly.cc b/dpoly.cc index 22150dd..23668fc 100644 --- a/dpoly.cc +++ b/dpoly.cc @@ -119,6 +119,7 @@ void dpoly::div(const dpoly& d, mpq_t *count, const mpq_t& factor) value_multiply(mpq_numref(tmp), coeff->p[len-1 - i], mpq_numref(factor)); value_multiply(mpq_denref(tmp), denom->p[len-1 - i], mpq_denref(factor)); mpq_add(count[i], count[i], tmp); + mpq_canonicalize(count[i]); } mpq_clear(tmp); diff --git a/genfun.cc b/genfun.cc index a52df96..309bd89 100644 --- a/genfun.cc +++ b/genfun.cc @@ -5,6 +5,7 @@ #include #include #include "conversion.h" +#include "counter.h" #include "genfun_constructor.h" #include "mat_util.h" #include "matrix_read.h" @@ -899,7 +900,8 @@ bool gen_fun::summate(Value *sum) const if ((*i)->d.power.NumRows() > maxlen) maxlen = (*i)->d.power.NumRows(); - infinite_icounter cnt((*term.begin())->d.power.NumCols(), maxlen); + infinite_counter cnt((*term.begin())->d.power.NumCols(), maxlen); + cnt.init(context); for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) cnt.reduce((*i)->n.coeff, (*i)->n.power, (*i)->d.power); diff --git a/reducer.cc b/reducer.cc index bae9db6..ff01bc2 100644 --- a/reducer.cc +++ b/reducer.cc @@ -463,47 +463,3 @@ void icounter::base(const QQ& c, const vec_ZZ& num, const mat_ZZ& den_f) } mpq_add(count, count, tcount); } - -void infinite_icounter::base(const QQ& c, const vec_ZZ& num, const mat_ZZ& den_f) -{ - Value tmp; - mpq_t factor; - mpq_init(factor); - value_init(tmp); - zz2value(c.n, tmp); - value_assign(mpq_numref(factor), tmp); - zz2value(c.d, tmp); - value_assign(mpq_denref(factor), tmp); - value_clear(tmp); - - unsigned len = den_f.NumRows(); // number of factors in den - - if (len > 0) { - int r; - vec_ZZ den_s; - den_s.SetLength(len); - assert(num.length() == 1); - ZZ num_s = num[0]; - - for (r = 0; r < len; ++r) - den_s[r] = den_f[r][0]; - int sign = (len % 2) ? -1 : 1; - - zz2value(num_s, tz); - dpoly n(len, tz); - zz2value(den_s[0], tz); - dpoly D(len, tz, 1); - for (int k = 1; k < len; ++k) { - zz2value(den_s[k], tz); - dpoly fact(len, tz, 1); - D *= fact; - } - - if (sign == -1) - value_oppose(mpq_numref(factor), mpq_numref(factor)); - n.div(D, count, factor); - } else - mpq_add(count[0], count[0], factor); - - mpq_clear(factor); -} diff --git a/reducer.h b/reducer.h index 9e01bce..fc7226a 100644 --- a/reducer.h +++ b/reducer.h @@ -127,36 +127,4 @@ struct icounter : public ireducer { } }; -/* An incremental counter for possibly infinite sets. - * Rather than just keeping track of the constant term - * of the Laurent expansions, we also keep track of the - * coefficients of negative powers. - * If any of these is non-zero, then the counted set is infinite. - */ -struct infinite_icounter : public ireducer { - /* an array of coefficients; count[i] is the coeffient of - * the term with power -i. - */ - mpq_t *count; - unsigned len; - Value tz; - - infinite_icounter(unsigned dim, unsigned maxlen) : ireducer(dim), len(maxlen+1) { - /* Not sure whether it works for dim != 1 */ - assert(dim == 1); - count = new mpq_t[len]; - for (int i = 0; i < len; ++i) - mpq_init(count[i]); - lower = 1; - value_init(tz); - } - ~infinite_icounter() { - for (int i = 0; i < len; ++i) - mpq_clear(count[i]); - delete [] count; - value_clear(tz); - } - virtual void base(const QQ& c, const vec_ZZ& num, const mat_ZZ& den_f); -}; - #endif diff --git a/testlib.cc b/testlib.cc index c3253f6..2418d42 100644 --- a/testlib.cc +++ b/testlib.cc @@ -174,6 +174,77 @@ static Matrix *matrix_read_from_str(const char *s) return Matrix_Read(str); } +static int test_infinite_counter(struct barvinok_options *options) +{ + Matrix *M = matrix_read_from_str("1 3\n 1 1 0\n"); + Polyhedron *ctx = Constraints2Polyhedron(M, options->MaxRays); + Matrix_Free(M); + + /* (1 -1/2 x^5 - 1/2 x^7)/(1-x) */ + infinite_counter *cnt = new infinite_counter(1, 1); + cnt->init(ctx); + vec_QQ n_coeff; + mat_ZZ n_power; + mat_ZZ d_power; + set_from_string(n_coeff, "[1/1 -1/2 -1/2]"); + set_from_string(n_power, "[[0][5][7]]"); + set_from_string(d_power, "[[1]]"); + cnt->reduce(n_coeff, n_power, d_power); + assert(value_cmp_si(mpq_numref(cnt->count[0]), 6) == 0); + assert(value_cmp_si(mpq_denref(cnt->count[0]), 1) == 0); + assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) == 0); + assert(value_cmp_si(mpq_denref(cnt->count[1]), 1) == 0); + delete cnt; + Polyhedron_Free(ctx); + + M = matrix_read_from_str("2 4\n 1 1 0 0\n 1 0 1 0\n"); + ctx = Constraints2Polyhedron(M, options->MaxRays); + Matrix_Free(M); + + /* (1 - xy)/((1-x)(1-xy)) */ + cnt = new infinite_counter(2, 3); + cnt->init(ctx); + set_from_string(n_coeff, "[1/1 -1/1]"); + set_from_string(n_power, "[[0 0][1 1]]"); + set_from_string(d_power, "[[1 0][1 1]]"); + cnt->reduce(n_coeff, n_power, d_power); + assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) != 0); + assert(value_cmp_si(mpq_numref(cnt->count[2]), 0) == 0); + assert(value_cmp_si(mpq_denref(cnt->count[2]), 1) == 0); + assert(value_cmp_si(mpq_numref(cnt->count[3]), 0) == 0); + assert(value_cmp_si(mpq_denref(cnt->count[3]), 1) == 0); + delete cnt; + + cnt = new infinite_counter(2, 2); + cnt->init(ctx); + set_from_string(n_coeff, "[-1/2 1/1 -1/3]"); + set_from_string(n_power, "[[2 6][3 6]]"); + d_power.SetDims(0, 2); + cnt->reduce(n_coeff, n_power, d_power); + assert(value_cmp_si(mpq_numref(cnt->count[0]), 1) == 0); + assert(value_cmp_si(mpq_denref(cnt->count[0]), 6) == 0); + assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) == 0); + assert(value_cmp_si(mpq_denref(cnt->count[1]), 1) == 0); + assert(value_cmp_si(mpq_numref(cnt->count[2]), 0) == 0); + assert(value_cmp_si(mpq_denref(cnt->count[2]), 1) == 0); + delete cnt; + + cnt = new infinite_counter(2, 2); + cnt->init(ctx); + set_from_string(n_coeff, "[1/1]"); + set_from_string(n_power, "[[0 11]]"); + set_from_string(d_power, "[[0 1]]"); + cnt->reduce(n_coeff, n_power, d_power); + assert(value_cmp_si(mpq_numref(cnt->count[1]), 0) != 0); + assert(value_cmp_si(mpq_numref(cnt->count[2]), 0) == 0); + assert(value_cmp_si(mpq_denref(cnt->count[2]), 1) == 0); + delete cnt; + + Polyhedron_Free(ctx); + + return 0; +} + static int test_series(struct barvinok_options *options) { Matrix *M = matrix_read_from_str( @@ -392,6 +463,7 @@ int main(int argc, char **argv) test_specialization(options); test_lattice_points(options); test_icounter(options); + test_infinite_counter(options); test_series(options); test_todd(options); test_bernoulli(options); -- 2.11.4.GIT