From 93cd92b37142ad0fd8baf7e1644ab5775ce9a7f5 Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Fri, 25 Apr 2008 09:55:06 +0200 Subject: [PATCH] Add Laurent expansion based summation --- Makefile.am | 2 + barvinok.cc | 3 + barvinok/options.h | 1 + configure.in | 2 + doc/barvinok.bib | 7 + doc/implementation.tex | 310 ++++++++++++++++- doc/mydefs.sty | 1 + laurent.cc | 823 ++++++++++++++++++++++++++++++++++++++++++++ laurent.h | 14 + m4/ac_cxx_gnucxx_hashmap.m4 | 39 +++ options.c | 10 +- testlib.cc | 35 ++ 12 files changed, 1238 insertions(+), 9 deletions(-) create mode 100644 laurent.cc create mode 100644 laurent.h create mode 100644 m4/ac_cxx_gnucxx_hashmap.m4 diff --git a/Makefile.am b/Makefile.am index 6fa24d7..28cb02b 100644 --- a/Makefile.am +++ b/Makefile.am @@ -117,6 +117,8 @@ libbarvinok_core_la_SOURCES = \ lattice_point.h \ lattice_width.c \ lattice_width.h \ + laurent.cc \ + laurent.h \ normalization.c \ normalization.h \ options.c \ diff --git a/barvinok.cc b/barvinok.cc index 7272ea8..ef308b4 100644 --- a/barvinok.cc +++ b/barvinok.cc @@ -20,6 +20,7 @@ #include "decomposer.h" #include "euler.h" #include "lattice_point.h" +#include "laurent.h" #include "reduce_domain.h" #include "remove_equalities.h" #include "scale.h" @@ -1569,6 +1570,8 @@ evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options) { if (options->summation == BV_SUM_EULER) return euler_summate(e, nvar, options); + else if (options->summation == BV_SUM_LAURENT) + return laurent_summate(e, nvar, options); else if (options->summation == BV_SUM_BERNOULLI) return Bernoulli_sum_evalue(e, nvar, options); else diff --git a/barvinok/options.h b/barvinok/options.h index 7401f00..699f71c 100644 --- a/barvinok/options.h +++ b/barvinok/options.h @@ -93,6 +93,7 @@ struct barvinok_options { #define BV_SUM_BARVINOK 0 #define BV_SUM_EULER 1 #define BV_SUM_BERNOULLI 2 + #define BV_SUM_LAURENT 3 int summation; #define BV_CHAMBERS_POLYLIB 0 diff --git a/configure.in b/configure.in index d772c0b..f3cc279 100644 --- a/configure.in +++ b/configure.in @@ -40,6 +40,8 @@ if test "x$bv_cv_shared_barvinok" != "xyes" ; then BV_LDFLAGS="$BV_LDFLAGS -static" fi +AC_CXX_GNUCXX_HASHMAP + AC_ARG_WITH(default-prefix, AS_HELP_STRING([--with-default-prefix=DIR], [Default installation prefix of optional packages])) diff --git a/doc/barvinok.bib b/doc/barvinok.bib index 0c07406..d5dc0cc 100644 --- a/doc/barvinok.bib +++ b/doc/barvinok.bib @@ -1005,3 +1005,10 @@ MRREVIEWER = {P. McMullen}, publisher = "Springer-Verlag", address = "Berlin", } + +@misc{Baldoni2008, + title = "Sum over lattice points of a polygon with iterated Laurent series. User's guide", + author = "Velleda Baldoni and Nicole Berline and Mich\'ele Vergne", + month = mar, + year = 2008, +} diff --git a/doc/implementation.tex b/doc/implementation.tex index 91682c2..094dadd 100644 --- a/doc/implementation.tex +++ b/doc/implementation.tex @@ -1477,6 +1477,7 @@ c_l = \frac{d_l}{b_0^{l+1}} $$ \subsection{Specialization through exponential substitution} +\label{s:exponential} This section draws heavily from \shortciteN{Koeppe2006experiments}. @@ -1635,16 +1636,16 @@ $$ $$ with $\alpha_i = \sum_{k=1}^{r} a_k^i$. -\begin{example} +\begin{example} \label{ex:todd} Consider the \rgf/ -\begin{multline*} +$$ \f T x = \frac{x_1^2}{(1-x_1^{-1})(1-x_1^{-1}x_2)} + \frac{x_2^2}{(1-x_2^{-1})(1-x_1 x_2^{-1})} -+ {} \\ ++ \frac1{(1-x_1)(1-x_2)} -\end{multline*} +$$ from \shortciteN[Example~39]{Verdoolaege2005PhD}. Since this is a 2-dimensional problem, we first compute the first 3 Todd polynomials (evaluated at $-1$), @@ -2328,10 +2329,12 @@ we have~\cite[(13)]{Berline2006local} \end{equation} with $y = \sps \xi r$ and $b(n,t)$ the \ai{Bernoulli polynomial}s defined by the generating series -$$ +\begin{equation} +\label{eq:bernoulli} +\Todd(t,y) = \frac{e^{ty} y}{e^y - 1} = \sum_{n=0}^\infty \frac{b(n,t)}{n!} y^n . -$$ +\end{equation} The constant terms of these Bernoulli polynomials are the \ai{Bernoulli number}s. @@ -3423,6 +3426,301 @@ the integer points of $P(n)$} \end{example} +\subsection{Summation through exponential substitution and Laurent expansions} + +This section was inspired by \shortciteN{Baldoni2008}. + +Let $f(\vec x)$ be the generating function of a polytope $P$, +i.e., +$$ +f(\vec x) = \sum_{\vec t \in P \cap \ZZ^d} \vec x^{\vec t} +. +$$ +Substituting $\vec x = \vec e^{\vec y}$, we obtain +$$ +f(\vec e^{\vec y}) = \sum_{\vec t \in P \cap \ZZ^d} e^{\sp t y} += +\sum_{\vec t \in P \cap \ZZ^d} + \sum_{\vec n \ge \vec 0} \frac{\vec t^{\vec n} \vec y^{\vec n}}{\vec n!} += +\sum_{\vec n \ge \vec 0} + \left(\sum_{\vec t \in P \cap \ZZ^d} \vec t^{\vec n}\right) + \frac{\vec y^{\vec n}}{\vec n!} +, +$$ +with $\vec n! = n_1! n_2! \cdots n_d!$. +We observe that the sum of the monomial $\vec t^{\vec n}$ +over the integer points in $P$ is equal to $\vec n!$ times the coefficient +of the $\vec y^{\vec n}$ term in the Taylor expansion of $f(\vec e^{\vec y})$. + +As in the case of unweighted counting (see \autoref{s:exponential}), +we have to add the coefficients +of these monomials in the Laurent expansions of the terms in \eqref{eq:rgf}. +However, unlike the case of unweighted counting, we cannot transform +this problem to a univariate problem and computing the coefficient +of a monomial in the Laurent expansions does not reduces to computing +the coefficient of a single higher-degree monomial in a Taylor expansion. + +Consider now one of the terms $g(\vec x) = f_{ik}(\vec x)$ in~\eqref{eq:rgf}, +$$ +g(\vec e^{\vec y}) = + \frac{e^{\sum_{j=1}^d s_j(\vec p) \sp{b_j}y}} + {\prod_{j=1}^{d}\left(1-e^{\sp{b_j}y}\right)} +, +$$ +with $\vec w_{ij}(\vec p) = \sum_{j=1}^d s_j(\vec p) \vec b_j$ written +in terms of the $\vec b_j$, which are assumed to form a basis +and where we have made explicit the only place where the +parameter $\vec p$ appear. +We rewrite this equation as +\begin{equation} +\label{eq:laurent} +g(\vec e^{\vec y}) = + \left(\prod_{j=1}^d \frac{-1}{\sp{b_j}y}\right) + \left(\prod_{j=1}^d \frac{-\sp{b_j}y \, e^{s_j(\vec p)\sp{b_j}y}} + {1-e^{\sp{b_j}y}} \right) +. +\end{equation} +The second factor is analytic and is a product of +generating functions +$\Todd(s_j(\vec p), \sp{b_j}y)$ +of \ai{Bernoulli polynomial}s~\eqref{eq:bernoulli}. +Plugging in these expressions, we find +\begin{align} +\Todd(s_j(\vec p), \sp{b_j}y) +&= \frac{-\sp{b_j}y e^{s_j(\vec p)\sp{b_j}y}} + {1-e^{\sp{b_j}y}} +\nonumber +\\ +&= \sum_{n=0}^\infty \frac{b(n,s_j(\vec p))}{n!} \sp{b_j}y^n +\nonumber +\\ +&= \sum_{\vec k \ge \vec 0} + \frac{b(\sum k_i,s_j(\vec p))}{(\sum k_i)!} + {\sum k_i \choose \vec k} \vec b_j^{\vec k} \vec y^{\vec k} +\nonumber +\\ +&= \sum_{\vec k \ge \vec 0} + \frac{b(\sum k_i,s_j(\vec p))}{\prod_i k_i!} + \vec b_j^{\vec k} \vec y^{\vec k} +\label{eq:laurent:todd} +, +\end{align} +with +$$ +{\sum k_i \choose \vec k} = +{\sum k_i \choose k_1, k_2, \ldots k_d} = +\frac{(\sum k_i)!}{\vec k!} = +\prod_{i=1}^d {\sum_{j=1}^i k_j \choose k_i} +$$ +the \ai{multinomial coefficient}s. +For the first factor, we compute the Laurent expansion +of each of its factors, +\begin{align} +\frac{-1}{\sp{b_j}y} +&= \frac{-1}{\sum_{k=f}^d b_{jk} y_k} +\nonumber +\\ +&= \frac{-1}{b_{jf} y_f\left(1 + \frac{\sum_{k=f+1}^d b_{jk} y_k} + {b_{jf} y_f}\right)} +\nonumber +\\ +&= \frac{-1}{b_{jf} y_f} + \sum_{n=0}^\infty (-1)^n \left(\frac{\sum_{k=f+1}^d b_{jk} y_k} + {b_{jf} y_f}\right)^n +\nonumber +\\ +&= \sum_{\vec n \ge \vec 0} + {\sum n_k \choose \vec n} + (-1)^{\sum n_k+1} + \frac{{\vec b'_j}^{\vec n}}{b_{jf}^{\sum n_k+1}} + \frac{{\vec y'}^{\vec n}}{y_f^{\sum n_k+1}} +\label{eq:laurent:reciprocal} +, +\end{align} +where $b_{jf}$ is the first non-zero coefficient of $\vec b_j$ +and the vector $\vec b_j'$ contains +the subsequent $d-f$ coefficients of $\vec b_j$. + +Given a polynomial +$$ +q(\vec y, \vec p) = \sum_{\vec m} \beta_{\vec m}(\vec p) \, \vec y^{\vec m} +$$ +that we wish to sum over the integer points of a polytope $P$, we perform +the following operations for each unimodular cone in the decomposition +of each vertex cone. +\begin{itemize} +\item For each $\vec m$ with $\beta_{\vec m}(\vec p) \ne 0$ +\begin{itemize} +\item Compute all sums +$\vec N = \sum_{j=1}^d (\vec 0, -\sum_k n_{jk}-1, \vec n_j)$ +of exponents from \eqref{eq:laurent:reciprocal} such that +$\vec N \le \vec m$ and compute the corresponding coefficient $\gamma_{\vec N}$ +in the product of Laurent series by enumerating all combinations +of $\vec n_j$ leading to the same $\vec N$. +Note that there are only a finite number of $\vec N$ satisfying this constraint +since $\sum N_k = -d$. +By reordering the variables such that the highest exponents occurs +for the first variable, the number of $\vec N$ can be reduced. +\item For each of these $\vec N$ +\begin{itemize} +\item Compute the coefficient $\delta_{\vec m - \vec N}(\vec p)$ of +$\vec y^{\vec m - \vec N}$ in the product of +Taylor expansions~\eqref{eq:laurent:todd}. +\end{itemize} +\end{itemize} +\item The contribution of this cone is the sum of +$$ +\vec m! \, \alpha \, \beta_{\vec m}(\vec p) \, \gamma_{\vec N} \, + \delta_{\vec m - \vec N}(\vec p) +$$ +over all considered $\vec m$ and $\vec N$. +\end{itemize} +Within each vertex cone computation, the coefficients +$\gamma_{\vec N}$ and $\delta_{\vec m - \vec N}(\vec p)$ +only need to be computed once. + +\begin{example} +Consider once more the \rgf/ +$$ +\f T x = +\frac{x_1^2}{(1-x_1^{-1})(1-x_1^{-1}x_2)} ++ +\frac{x_2^2}{(1-x_2^{-1})(1-x_1 x_2^{-1})} ++ +\frac1{(1-x_1)(1-x_2)} +$$ +from \shortciteN[Example~39]{Verdoolaege2005PhD} +and Example~\ref{ex:todd}. +Assume we want to compute +$$ +\sum_{\vec y \in T\cap \ZZ^2} y_1^2 + y_2^2 +. +$$ +We will need the following \ai{Bernoulli polynomials} +\begin{align*} +b(0,s) & = 1 +\\ +b(1,s) &= \frac 1 2 \left(-1 + 2 s \right) +\\ +b(2,s) &= \frac 1 6 \left(1 -6 s + 6 s^2 \right) +\\ +b(3,s) &= \frac 1 2 \left(s -3 s^2 + 2 s^3\right) +\\ +b(4,s) & = \frac 1{30} \left(-1 + 30 s^2 -60 s^3 + 30 s^4\right) +\end{align*} +For the first term, substitution yields +\begin{align*} +h(\vec y) +=& +\frac 1{y_1} +\frac 1{y_1-y_2} +\frac {y_1 e^{(-2)(-y_1)}}{1-e^{-y_1}} +\frac {(y_1-y2) e^{0(-y_1+y_2)}}{1-e^{-y_1+y_2}} +\\ +=& +\frac 1{y_1} +\left( + \frac 1{y_1} \left( + 1 + \frac{y_2}{y_1} + \frac{y_2^2}{y_1^2} + \cdots + \right) +\right) +\\ +& +\left( +1 + \frac{b(1,-2)}1 (-y_1) + \frac{b(2,-2)}2 (-y_1)^2 + + \frac{b(3,-2)}{3!} (-y_1)^3 + \frac{b(4,-2)}{4!} (-y_1)^4 + \cdots +\right) +\\ +& +\left( +1 + \frac{-1}2 (-y_1+y_2) + \frac{1}{12} (-y_1+y_2)^2 + + 0 (-y_1+y_2)^3 + + \frac{1}{720} (-y_1+y_2)^4 + \cdots +\right) +\end{align*} +We obtain the following results: +{\renewcommand\arraystretch{2} +$$ +\begin{array}{rrr@{}lrr@{}lc} +\vec m & \vec N & \gamma_{\vec N} & \vec y^{\vec N} & +\vec m - \vec N & \delta_{\vec m - \vec N} & \vec y^{\vec m - \vec N} & +\vec m! \alpha \beta_{\vec m} \gamma_{\vec N} \delta_{\vec m - \vec N} +\\ +(2,0) & (-2,0) & 1 & y_1^{-2} & (4,0) + & \displaystyle\frac{721}{240} & y_1^4 + & \displaystyle\frac{721}{120} +\\ +(0,2) & (-2,0) & 1 & y_1^{-2} & (2,2) + & \displaystyle\frac{179}{720} & y_1^2 y_2^2 + & \displaystyle\frac{179}{360} +\\ + & (-3,1) & 1 & y_1^{-3} y_2 & (3,1) + & \displaystyle-\frac{211}{120} & y_1^3 y_1 + & \displaystyle-\frac{211}{60} +\\ + & (-4,2) & 1 & y_1^{-4} y_2^2 & (4,0) + & \displaystyle\frac{721}{240} & y_1^4 + & \displaystyle\frac{721}{120} +\end{array} +$$ +}% +For the second term, we similarly obtain +{\renewcommand\arraystretch{2} +$$ +\begin{array}{rrr@{}lrr@{}lc} +\vec m & \vec N & \gamma_{\vec N} & \vec y^{\vec N} & +\vec m - \vec N & \delta_{\vec m - \vec N} & \vec y^{\vec m - \vec N} & +\vec m! \alpha \beta_{\vec m} \gamma_{\vec N} \delta_{\vec m - \vec N} +\\ +(2,0) & (-1,-1) & -1 & y_1^{-1} y_2^{-1} & (3,1) + & \displaystyle\frac{1}{180} & y_1^3 y_1 + & \displaystyle-\frac{1}{90} +\\ + & (-2,0) & -1 & y_1^{-2} & (4,0) + & \displaystyle-\frac{1}{720} & y_1^4 + & \displaystyle\frac{1}{360} +\\ +(0,2) & (-1,-1) & -1 & y_1^{-1} y_2^{-1} & (1,3) + & \displaystyle-\frac{211}{120} & y_1 y_2^3 + & \displaystyle\frac{211}{60} +\\ + & (-2,0) & -1 & y_1^{-3} y_2 & (2,2) + & \displaystyle\frac{179}{720} & y_1^3 y_2 + & \displaystyle-\frac{179}{360} +\\ + & (-3,1) & -1 & y_1^{-3} y_2 & (3,1) + & \displaystyle\frac{1}{180} & y_1^3 y_1 + & \displaystyle-\frac{1}{90} +\\ + & (-4,2) & -1 & y_1^{-4} y_2^2 & (4,0) + & \displaystyle-\frac{1}{720} & y_1^4 + & \displaystyle\frac{1}{360} +\end{array} +$$ +}% +Finally, for the third term, we obtain +$$ +\begin{array}{rrr@{}lrr@{}lc} +\vec m & \vec N & \gamma_{\vec N} & \vec y^{\vec N} & +\vec m - \vec N & \delta_{\vec m - \vec N} & \vec y^{\vec m - \vec N} & +\vec m! \alpha \beta_{\vec m} \gamma_{\vec N} \delta_{\vec m - \vec N} +\\ +(2,0) & (-1,-1) & -1 & y_1^{-1} y_2^{-1} & (3,1) + & 0 & y_1^3 y_1 + & 0 +\\ +(0,2) & (-1,-1) & -1 & y_1^{-1} y_2^{-1} & (1,3) + & 0 & y_1 y_2^3 + & 0 +\end{array} +$$ +Adding up all contributions in the final columns of these tables, +we obtain a grand total of +$$ +12. +$$ +\end{example} \subsection{Conversion to ``standard form''} \label{s:standard} diff --git a/doc/mydefs.sty b/doc/mydefs.sty index c761120..4940e5e 100644 --- a/doc/mydefs.sty +++ b/doc/mydefs.sty @@ -113,6 +113,7 @@ \DeclareMathOperator{\sgn}{sgn} \DeclareMathOperator{\vol}{vol} \DeclareMathOperator{\todd}{td} +\DeclareMathOperator{\Todd}{Todd} \DeclareMathOperator{\width}{width} \def\sm#1{ diff --git a/laurent.cc b/laurent.cc new file mode 100644 index 0000000..314778c --- /dev/null +++ b/laurent.cc @@ -0,0 +1,823 @@ +#include +#include +#include +#include +#include "bernoulli.h" +#include "binomial.h" +#include "conversion.h" +#include "decomposer.h" +#include "lattice_point.h" +#include "laurent.h" +#include "param_util.h" +#include "power.h" +#include "reduce_domain.h" +#include "config.h" + +using std::cerr; +using std::ostream; +using std::endl; +using std::vector; + +#ifdef HAVE_GNUCXX_HASHMAP + +#include + +#define HASH_MAP __gnu_cxx::hash_map + +namespace __gnu_cxx +{ + template<> struct hash< const std::vector > + { + size_t operator()( const std::vector& x ) const + { + unsigned long __h = 0; + for (int i = 0; i < x.size(); ++i) + __h = 5 * __h + x[i]; + return size_t(__h); + } + }; +} + +#else + +#warning "no hash_map available" +#include +#define HASH_MAP std::map + +#endif + +#define ALLOC(type) (type*)malloc(sizeof(type)) +#define ALLOCN(type,n) (type*)malloc((n) * sizeof(type)) + +static ostream & operator<< (ostream & os, const vector & v) +{ + os << "["; + for (int i = 0; i < v.size(); ++i) { + if (i) + os << ", "; + os << v[i]; + } + os << "]"; + return os; +} + +/* Compares the position of the first non-zero element + * in both vectors and the second non-zero element + * if the position of the first non-zero element is the same. + */ +static int pos_cmp(const void *a, const void *b) +{ + int pa, pb; + const Vector *va = *(const Vector **)a; + const Vector *vb = *(const Vector **)b; + + pa = First_Non_Zero(va->p, va->Size); + pb = First_Non_Zero(vb->p, vb->Size); + + if (pa == pb) { + pa = First_Non_Zero(va->p+pa+1, va->Size-pa-1); + pb = First_Non_Zero(vb->p+pa+1, vb->Size-pa-1); + } + + return pa - pb; +} + +/* Represents the vertex and the rays of a vertex cone */ +struct vertex_cone { + unsigned dim; + /* The coefficients of the rays, ordered according + * to the first non-zero coefficient. + */ + Vector **coeff; + Matrix Rays; + + /* The powers of the coefficients of the rays */ + struct power ***coeff_power; + + /* The coordinates of the integer point in the fundamental + * parallelepiped, in the basis formed by the rays of + * the vertex cone. + */ + evalue **E_vertex; + unsigned max_power; + /* Bernoulli polynomials corresponding to each E_vertex */ + evalue ***bernoulli_t; + + vertex_cone(unsigned dim); + void init(const mat_ZZ &rays, Param_Vertices *V, unsigned max_power); + void clear(); + + const evalue *get_bernoulli(int n, int i); + + ~vertex_cone() { + for (int i = 0; i < dim; ++i) + Vector_Free(coeff[i]); + free(coeff); + delete [] E_vertex; + free(Rays.p); + for (int i = 0; i < dim; ++i) + delete [] coeff_power[i]; + delete [] coeff_power; + delete [] bernoulli_t; + } +}; + +vertex_cone::vertex_cone(unsigned dim) : dim(dim) +{ + E_vertex = new evalue *[dim]; + bernoulli_t = new evalue **[dim]; + + coeff = ALLOCN(Vector *, dim); + for (int i = 0; i < dim; ++i) + coeff[i] = Vector_Alloc(dim); + + Rays.NbRows = Rays.NbColumns = dim; + Rays.p = ALLOCN(Value *, dim); + + coeff_power = new struct power **[dim]; + for (int i = 0; i < dim; ++i) + coeff_power[i] = new struct power *[dim]; +} + +void vertex_cone::init(const mat_ZZ &rays, Param_Vertices *V, + unsigned max_power) +{ + unsigned nparam = V->Vertex->NbColumns - 2; + this->max_power = max_power; + + for (int i = 0; i < dim; ++i) + zz2values(rays[i], coeff[i]->p); + qsort(coeff, dim, sizeof(Vector *), pos_cmp); + + for (int i = 0; i < dim; ++i) { + for (int j = 0; j < dim; ++j) { + if (value_notzero_p(coeff[i]->p[j])) + coeff_power[i][j] = new struct power(coeff[i]->p[j], max_power); + else + coeff_power[i][j] = NULL; + } + } + + for (int i = 0; i < dim; ++i) + Rays.p[i] = coeff[i]->p; + Matrix *T = Transpose(&Rays); + Matrix *L = relative_coordinates(V, T); + Matrix_Free(T); + + for (int i = 0; i < dim; ++i) + E_vertex[i] = ceiling(L->p[i], V->Vertex->p[0][nparam+1], nparam, NULL); + Matrix_Free(L); + + for (int j = 0; j < dim; ++j) { + bernoulli_t[j] = new evalue *[max_power]; + for (int k = 0; k < max_power; ++k) + bernoulli_t[j][k] = NULL; + } +} + +/* + * Returns b(n, E_vertex[i]) + */ +const evalue *vertex_cone::get_bernoulli(int n, int i) +{ + assert(n > 0); + if (!bernoulli_t[i][n-1]) { + struct poly_list *bernoulli = bernoulli_compute(n); + bernoulli_t[i][n-1] = evalue_polynomial(bernoulli->poly[n], + E_vertex[i]); + } + return bernoulli_t[i][n-1]; +} + +void vertex_cone::clear() +{ + for (int i = 0; i < dim; ++i) + if (E_vertex[i]) + evalue_free(E_vertex[i]); + + for (int i = 0; i < dim; ++i) { + for (int j = 0; j < max_power; ++j) { + if (bernoulli_t[i][j]) + evalue_free(bernoulli_t[i][j]); + } + delete [] bernoulli_t[i]; + } + + for (int i = 0; i < dim; ++i) { + for (int j = 0; j < dim; ++j) + if (coeff_power[i][j]) + delete coeff_power[i][j]; + } +} + +struct todd_product { + vertex_cone &vc; + unsigned dim; + /* The non-zero coefficients in the rays of the vertex cone */ + vector< vector > mask; + /* For each ray, the power of each variable it contributes */ + vector< vector > selected; + /* The powers of each variable that still need to be selected */ + vector left; + /* For each variable, the last ray that has a non-zero coefficient */ + vector last_level; + + HASH_MAP, const evalue *> cache; + + todd_product(vertex_cone &vc); + evalue *add(evalue *sum); + const evalue *get_coefficient(const vector &powers); + + ~todd_product() { + HASH_MAP, const evalue *>::iterator j; + for (j = cache.begin(); j != cache.end(); ++j) + if ((*j).second) + evalue_free(const_cast((*j).second)); + } +}; + +todd_product::todd_product(vertex_cone &vc) : vc(vc) +{ + dim = vc.dim; + for (int i = 0; i < dim; ++i) { + mask.push_back(vector(dim)); + selected.push_back(vector(dim)); + } + last_level = vector(dim); + for (int i = 0; i < dim; ++i) { + for (int j = 0; j < dim; ++j) { + if (vc.coeff_power[i][j]) { + mask[i][j] = 1; + last_level[j] = i; + } + } + } +} + +void multinomial(const vector &k, Value *m) +{ + int s = 0; + value_set_si(*m, 1); + for (int i = 0; i < k.size(); ++i) { + if (k[i] == 0) + continue; + s += k[i]; + value_multiply(*m, *m, *binomial(s, k[i])); + } +} + +/* Add coefficient of selected powers of variables to sum + * and return the result. + * The contribution of each ray j of the vertex cone is + * + * ( \sum k ) + * b(\sum k, \ceil{v}) ( k_1, \ldots, k_n ) c_1^{k_1} \cdots c_n^{k_n} + * + * with k_i the selected powers, c_i the coefficients of the ray + * and \ceil{v} the coordinate E_vertex[j] corresponding to the ray. + */ +evalue *todd_product::add(evalue *sum) +{ + evalue *c = NULL; + for (int i = 0; i < dim; ++i) { + int s = 0; + evalue *f = ALLOC(evalue); + value_init(f->d); + evalue_set_si(f, 1, 1); + for (int j = 0; j < dim; ++j) { + if (!selected[i][j]) + continue; + value_multiply(f->x.n, f->x.n, + *(*vc.coeff_power[i][j])[selected[i][j]]); + value_multiply(f->d, f->d, *factorial(selected[i][j])); + s += selected[i][j]; + } + if (s > 0) + emul(vc.get_bernoulli(s, i), f); + if (!c) + c = f; + else { + emul(f, c); + evalue_free(f); + } + } + if (!sum) + sum = c; + else { + eadd(c, sum); + evalue_free(c); + } + return sum; +} + +static int first_non_zero(const vector& row) +{ + for (int i = 0; i < row.size(); ++i) + if (row[i] != 0) + return i; + return -1; +} + +/* Return coefficient of the term with exponents powers by + * iterating over all combinations of exponents for each ray + * that sum up to the given exponents. + */ +const evalue *todd_product::get_coefficient(const vector &powers) +{ + evalue *c = NULL; + + HASH_MAP, const evalue *>::iterator i; + i = cache.find(powers); + if (i != cache.end()) + return (*i).second; + + for (int i = 0; i < vc.dim; ++i) + for (int j = 0; j < vc.dim; ++j) + selected[i][j] = 0; + + left = powers; + int nz = first_non_zero(left); + int level = last_level[nz]; + int p = nz; + while (level >= 0) { + if (mask[level][p] && left[p] > 0) { + selected[level][p]++; + left[p]--; + /* Fill out remaining powers and make sure we backtrack from + * the right position. + */ + for (int i = 0; i < vc.dim; ++i) { + if (left[i] == 0) + continue; + selected[last_level[i]][i] += left[i]; + left[i] = 0; + if (last_level[i] >= level) { + level = last_level[i]; + p = i; + } + } + + c = add(c); + } + if (selected[level][p]) { + left[p] += selected[level][p]; + selected[level][p] = 0; + } + if (--p < 0) { + --level; + p = dim-1; + } + } + cache[powers] = c; + return c; +} + +/* Maintains the coefficients of the reciprocals of the + * (negated) rays of the vertex cone vc. + */ +struct reciprocal { + vertex_cone &vc; + + /* can_borrow_from[i][j] = 1 if there is a ray + * with first non-zero coefficient i and a subsequent + * non-zero coefficient j. + */ + vector< vector > can_borrow_from; + /* Total exponent that a variable can borrow from subsequent vars */ + vector can_borrow; + /* has_borrowed_from[i][j] is the exponent borrowed by i from j */ + vector< vector > has_borrowed_from; + /* Total exponent borrowed from subsequent vars */ + vector has_borrowed; + /* The last variable that can borrow from subsequent vars */ + int last; + + /* Position of the first non-zero coefficient in each ray. */ + vector neg; + + /* Power without any "borrowing" */ + vector base_power; + /* Power after "borrowing" */ + vector power; + + /* The non-zero coefficients in the rays of the vertex cone, + * except the first. + */ + vector< vector > mask; + /* For each ray, the (positive) power of each variable it contributes */ + vector< vector > selected; + /* The powers of each variable that still need to be selected */ + vector left; + + HASH_MAP, const evalue *> cache; + + reciprocal(vertex_cone &vc); + void start(vector &power); + int next(); + + evalue *add(evalue *sum); + const evalue *get_coefficient(); + ~reciprocal() { + HASH_MAP, const evalue *>::iterator j; + for (j = cache.begin(); j != cache.end(); ++j) + if ((*j).second) + evalue_free(const_cast((*j).second)); + } +}; + +reciprocal::reciprocal(vertex_cone &vc) : vc(vc) +{ + for (int i = 0; i < vc.dim; ++i) { + can_borrow_from.push_back(vector(vc.dim)); + has_borrowed_from.push_back(vector(vc.dim)); + mask.push_back(vector(vc.dim)); + selected.push_back(vector(vc.dim)); + } + can_borrow = vector(vc.dim); + has_borrowed = vector(vc.dim); + neg = vector(vc.dim); + left = vector(vc.dim); + + for (int i = 0; i < vc.dim; ++i) { + int pos = First_Non_Zero(vc.coeff[i]->p, vc.coeff[i]->Size); + neg[i] = pos; + for (int j = pos+1; j < vc.dim; ++j) + if (value_notzero_p(vc.coeff[i]->p[j])) { + mask[i][j] = 1; + can_borrow_from[neg[i]][j] = 1; + } + } +} + +/* Initialize power to the exponent of the todd product + * required to compute the coefficient in the full product + * of the term with exponent req_power, without any + * "borrowing". + */ +void reciprocal::start(vector &req_power) +{ + power = req_power; + for (int j = 0; j < vc.dim; ++j) + power[neg[j]]++; + + base_power = power; + + last = -1; + for (int i = 0; i < vc.dim; ++i) { + can_borrow[i] = 0; + has_borrowed[i] = 0; + for (int j = i+1; j < vc.dim; ++j) { + has_borrowed_from[i][j] = 0; + if (can_borrow_from[i][j]) + can_borrow[i] += power[j]; + } + if (can_borrow[i]) + last = i; + } +} + +/* Set power to the next exponent of the todd product required + * and return 1 as long as there is any such exponent left. + */ +int reciprocal::next() +{ + int p = last; + while (p >= 0) { + if (has_borrowed[p] < can_borrow[p]) { + int j; + for (j = p+1; j < vc.dim; ++j) + if (can_borrow_from[p][j]) { + if (power[j] > 0) + break; + else if (has_borrowed_from[p][j]) { + power[j] += has_borrowed_from[p][j]; + power[p] -= has_borrowed_from[p][j]; + has_borrowed[p] -= has_borrowed_from[p][j]; + has_borrowed_from[p][j] = 0; + } + } + if (j < vc.dim) { + has_borrowed_from[p][j]++; + has_borrowed[p]++; + power[p]++; + power[j]--; + return 1; + } + } + if (has_borrowed[p]) { + for (int j = p+1; j < vc.dim; ++j) + if (has_borrowed_from[p][j]) { + power[j] += has_borrowed_from[p][j]; + has_borrowed_from[p][j] = 0; + } + power[p] -= has_borrowed[p]; + has_borrowed[p] = 0; + } + --p; + } + return 0; +} + +/* Add coefficient of selected powers of variables to sum + * and return the result. + * The contribution of each ray j of the vertex cone is + * + * ( K ) + * ( k_{f+1}, \ldots, k_n ) (-1)^{K+1} + * c_f^{-K-1} c_{f+1}^{k_{f+1}} \cdots c_n^{k_n} + * + * K = \sum_{i=f+1}^n k_i + */ +evalue *reciprocal::add(evalue *sum) +{ + evalue *t = NULL; + for (int i = 0; i < vc.dim; ++i) { + evalue *c = ALLOC(evalue); + value_init(c->d); + value_init(c->x.n); + value_set_si(c->d, 1); + multinomial(selected[i], &c->x.n); + int s = 0; + for (int j = 0; j < vc.dim; ++j) { + if (selected[i][j] == 0) + continue; + value_multiply(c->x.n, c->x.n, + *(*vc.coeff_power[i][j])[selected[i][j]]); + s += selected[i][j]; + } + value_multiply(c->d, c->d, *(*vc.coeff_power[i][neg[i]])[s+1]); + if (!(s % 2)) + value_oppose(c->x.n, c->x.n); + if (value_neg_p(c->d)) { + value_oppose(c->d, c->d); + value_oppose(c->x.n, c->x.n); + } + if (!t) + t = c; + else { + emul(c, t); + evalue_free(c); + } + } + if (!sum) + sum = t; + else { + eadd(t, sum); + evalue_free(t); + } + return sum; +} + +/* Return coefficient of the term with exponents powers by + * iterating over all combinations of exponents for each ray + * that sum up to the given exponents. + */ +const evalue *reciprocal::get_coefficient() +{ + evalue *c = NULL; + + for (int j = 0; j < vc.dim; ++j) + left[j] = base_power[j] - power[j]; + + HASH_MAP, const evalue *>::iterator i; + i = cache.find(left); + if (i != cache.end()) + return (*i).second; + + int nz = first_non_zero(left); + if (nz == -1) + return cache[power] = add(c); + if (left[nz] > 0) + return NULL; + + for (int i = 0; i < vc.dim; ++i) + for (int j = 0; j < vc.dim; ++j) + selected[i][j] = 0; + + int level = vc.dim-1; + int p = vc.dim-1; + while (level >= 0) { + int nz = first_non_zero(left); + if (nz < neg[level] || (nz == neg[level] && left[nz] > 0)) { + assert(p == vc.dim-1); + --level; + continue; + } + if (nz == neg[level] && mask[level][p]) { + selected[level][p]++; + left[p]--; + left[neg[level]]++; + int nz = first_non_zero(left); + if (nz == -1) + c = add(c); + if (left[nz] < 0) { + level = vc.dim-1; + p = vc.dim-1; + continue; + } + } + if (selected[level][p]) { + left[p] += selected[level][p]; + left[neg[level]] -= selected[level][p]; + selected[level][p] = 0; + } + if (--p < 0) { + --level; + p = vc.dim-1; + } + } + cache[left] = c; + + return c; +} + +/* A term in the input polynomial */ +struct term { + vector powers; + const evalue *coeff; +}; + +static void collect_terms_r(vector &terms, vector &powers, + const evalue *polynomial, unsigned nvar) +{ + if (EVALUE_IS_ZERO(*polynomial)) + return; + + if (value_zero_p(polynomial->d)) + assert(polynomial->x.p->type == ::polynomial); + if (value_notzero_p(polynomial->d) || polynomial->x.p->pos > nvar) { + struct term t; + t.powers = powers; + t.coeff = polynomial; + terms.push_back(t); + return; + } + + for (int i = polynomial->x.p->size-1; i >= 0; --i) { + powers[polynomial->x.p->pos-1] = i; + collect_terms_r(terms, powers, &polynomial->x.p->arr[i], nvar); + } +} + +/* Expand "polynomial" as a sum of powers of the "nvar" variables, + * collect the terms in "terms" and return the maximal total degree. + */ +static unsigned collect_terms(vector &terms, + const evalue *polynomial, unsigned nvar) +{ + vector powers(nvar); + for (int i = 0; i < nvar; ++i) + powers[i] = 0; + collect_terms_r(terms, powers, polynomial, nvar); + + unsigned max_degree = 0; + for (int i = 0; i < terms.size(); ++i) { + int sum = 0; + for (int j = 0; j < nvar; ++j) + sum += terms[i].powers[j]; + if (sum > max_degree) + max_degree = sum; + } + return max_degree; +} + +/* +static void dump_terms(const vector &terms) +{ + for (int i = 0; i < terms.size(); ++i) { + cerr << terms[i].powers; + print_evalue(stderr, terms[i].coeff, test); + } +} +*/ + +struct laurent_summator : public signed_cone_consumer, + public vertex_decomposer { + const evalue *polynomial; + unsigned dim; + vertex_cone vc; + vector terms; + evalue *result; + unsigned max_power; + + laurent_summator(const evalue *e, unsigned dim, Param_Polyhedron *PP) : + polynomial(e), dim(dim), vertex_decomposer(PP, *this), + vc(dim) { + max_power = dim + collect_terms(terms, polynomial, dim); + result = NULL; + } + ~laurent_summator() { + if (result) + evalue_free(result); + } + virtual void handle(const signed_cone& sc, barvinok_options *options); +}; + +static int first_non_zero(const vec_ZZ& row) +{ + for (int i = 0; i < row.length(); ++i) + if (row[i] != 0) + return i; + return -1; +} + +void laurent_summator::handle(const signed_cone& sc, barvinok_options *options) +{ + assert(sc.det == 1); + + vc.init(sc.rays, V, max_power); + reciprocal recip(vc); + todd_product tp(vc); + for (int i = 0; i < terms.size(); ++i) { + recip.start(terms[i].powers); + do { + const evalue *c = recip.get_coefficient(); + if (!c) + continue; + + const evalue *t = tp.get_coefficient(recip.power); + + evalue *f = evalue_dup(terms[i].coeff); + if (sc.sign < 0) + evalue_negate(f); + for (int j = 0; j < dim; ++j) + evalue_mul(f, *factorial(terms[i].powers[j])); + evalue_shift_variables(f, -dim); + emul(c, f); + emul(t, f); + if (!result) + result = f; + else { + eadd(f, result); + evalue_free(f); + } + } while (recip.next()); + } + vc.clear(); +} + +static evalue *summate_over_domain(const evalue *e, unsigned nvar, + Polyhedron *D, + struct barvinok_options *options) +{ + Polyhedron *U, *TC; + Param_Polyhedron *PP; + struct evalue_section *s; + int nd = -1; + Param_Domain *PD; + evalue *sum; + + U = Universe_Polyhedron(D->Dimension - nvar); + PP = Polyhedron2Param_Polyhedron(D, U, options); + + for (nd = 0, PD = PP->D; PD; ++nd, PD = PD->next); + s = ALLOCN(struct evalue_section, nd); + + TC = true_context(D, U, options->MaxRays); + FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, PD, rVD) + Param_Vertices *V; + laurent_summator ls(e, nvar, PP); + + FORALL_PVertex_in_ParamPolyhedron(V, PD, PP) // _i internal counter + ls.decompose_at_vertex(V, _i, options); + END_FORALL_PVertex_in_ParamPolyhedron; + + s[i].D = rVD; + s[i].E = ls.result; + ls.result = NULL; + END_FORALL_REDUCED_DOMAIN + Polyhedron_Free(TC); + Polyhedron_Free(U); + Param_Polyhedron_Free(PP); + + sum = evalue_from_section_array(s, nd); + free(s); + + return sum; +} + +evalue *laurent_summate(evalue *e, unsigned nvar, + struct barvinok_options *options) +{ + int i; + evalue *res; + + assert(nvar >= 0); + if (nvar == 0 || EVALUE_IS_ZERO(*e)) + return evalue_dup(e); + + assert(value_zero_p(e->d)); + assert(e->x.p->type == partition); + + res = evalue_zero(); + + for (i = 0; i < e->x.p->size/2; ++i) { + evalue *t; + t = summate_over_domain(&e->x.p->arr[2*i+1], nvar, + EVALUE_DOMAIN(e->x.p->arr[2*i]), options); + eadd(t, res); + evalue_free(t); + } + + return res; +} diff --git a/laurent.h b/laurent.h new file mode 100644 index 0000000..b0b0404 --- /dev/null +++ b/laurent.h @@ -0,0 +1,14 @@ +#include + +#if defined(__cplusplus) +extern "C" { +#endif + +struct barvinok_options; + +evalue *laurent_summate(evalue *e, unsigned nvar, + struct barvinok_options *options); + +#if defined(__cplusplus) +} +#endif diff --git a/m4/ac_cxx_gnucxx_hashmap.m4 b/m4/ac_cxx_gnucxx_hashmap.m4 new file mode 100644 index 0000000..e503956 --- /dev/null +++ b/m4/ac_cxx_gnucxx_hashmap.m4 @@ -0,0 +1,39 @@ +# =========================================================================== +# http://autoconf-archive.cryp.to/ac_cxx_gnucxx_hashmap.html +# =========================================================================== +# +# SYNOPSIS +# +# AC_CXX_GNUCXX_HASHMAP +# +# DESCRIPTION +# +# Test for the presence of GCC's hashmap STL extension. +# +# LAST MODIFICATION +# +# 2008-04-12 +# +# COPYLEFT +# +# Copyright (c) 2008 Patrick Mauritz +# +# Copying and distribution of this file, with or without modification, are +# permitted in any medium without royalty provided the copyright notice +# and this notice are preserved. + +AC_DEFUN([AC_CXX_GNUCXX_HASHMAP],[ +AC_CACHE_CHECK(whether the compiler supports __gnu_cxx::hash_map, +ac_cv_cxx_gnucxx_hashmap, +[AC_LANG_SAVE + AC_LANG_CPLUSPLUS + AC_TRY_COMPILE([#include +using __gnu_cxx::hash_map;], + [], + ac_cv_cxx_gnucxx_hashmap=yes, ac_cv_cxx_gnucxx_hashmap=no) + AC_LANG_RESTORE +]) +if test "$ac_cv_cxx_gnucxx_hashmap" = yes; then + AC_DEFINE(HAVE_GNUCXX_HASHMAP,,[define if the compiler supports __gnu_cxx::hash_map]) +fi +]) diff --git a/options.c b/options.c index 22e52b9..6a4a5b1 100644 --- a/options.c +++ b/options.c @@ -190,7 +190,7 @@ static struct argp_option barvinok_argp_options[] = { "[default: pip]", #endif }, - { "summation", BV_OPT_SUM, "barvinok|bernoulli|euler", 0, + { "summation", BV_OPT_SUM, "barvinok|bernoulli|euler|laurent", 0, "[default: barvinok]" }, { "bernstein-recurse", BV_OPT_RECURSE, "none|factors|intervals|full", 0, "[default: factors]" }, @@ -355,10 +355,14 @@ static error_t barvinok_parse_opt(int key, char *arg, struct argp_state *state) case BV_OPT_SUM: if (!strcmp(arg, "barvinok")) options->summation = BV_SUM_BARVINOK; - if (!strcmp(arg, "euler")) + else if (!strcmp(arg, "euler")) options->summation = BV_SUM_EULER; - if (!strcmp(arg, "bernoulli")) + else if (!strcmp(arg, "bernoulli")) options->summation = BV_SUM_BERNOULLI; + else if (!strcmp(arg, "laurent")) + options->summation = BV_SUM_LAURENT; + else + argp_error(state, "unknown summation method '%s'\n", arg); break; case BV_OPT_RECURSE: if (!strcmp(arg, "none")) diff --git a/testlib.cc b/testlib.cc index 4818924..586a2db 100644 --- a/testlib.cc +++ b/testlib.cc @@ -14,6 +14,7 @@ #include "hilbert.h" #include "hull.h" #include "ilp.h" +#include "laurent.h" #include "matrix_read.h" #include "remove_equalities.h" #include "config.h" @@ -716,6 +717,39 @@ int test_hull(struct barvinok_options *options) Matrix_Free(hull); } +static int test_laurent(struct barvinok_options *options) +{ + unsigned nvar, nparam; + char **all_vars; + evalue *e, *sum, *res; + + e = evalue_read_from_str(" x1 >= 0\n" + " x2 >= 0\n" + " -x1 -x2 + 2 >= 0\n" + "\n" + "(N + M * x2)\n", + "x1,x2", &all_vars, &nvar, &nparam, + options->MaxRays); + Free_ParamNames(all_vars, nvar+nparam); + + sum = laurent_summate(e, nvar, options); + + res = evalue_read_from_str("(6 * N + 4 * M)", + "", &all_vars, &nvar, &nparam, + options->MaxRays); + Free_ParamNames(all_vars, nvar+nparam); + + assert(value_zero_p(sum->d)); + assert(sum->x.p->type == partition); + assert(sum->x.p->size == 2); + + assert(eequal(res, &sum->x.p->arr[1])); + + evalue_free(e); + evalue_free(sum); + evalue_free(res); +} + int main(int argc, char **argv) { struct barvinok_options *options = barvinok_options_new_with_defaults(); @@ -736,5 +770,6 @@ int main(int argc, char **argv) test_hilbert(options); test_ilp(options); test_hull(options); + test_laurent(options); barvinok_options_free(options); } -- 2.11.4.GIT