[Mlir-commits] [mlir] [MLIR][Presburger] Implement function to evaluate the number of terms in a generating function. (PR #78078)
Arjun P
llvmlistbot at llvm.org
Sun Jan 21 12:21:22 PST 2024
================
@@ -245,3 +246,235 @@ QuasiPolynomial mlir::presburger::detail::getCoefficientInRationalFunction(
}
return coefficients[power].simplify();
}
+
+/// Substitute x_i = t^μ_i in one term of a generating function, returning
+/// a quasipolynomial which represents the exponent of the numerator
+/// of the result, and a vector which represents the exponents of the
+/// denominator of the result.
+/// v represents the affine functions whose floors are multiplied by the
+/// generators, and ds represents the list of generators.
+std::pair<QuasiPolynomial, std::vector<Fraction>>
+substituteMuInTerm(unsigned numParams, ParamPoint v, std::vector<Point> ds,
+ Point mu) {
+ unsigned numDims = mu.size();
+ for (const Point &d : ds)
+ assert(d.size() == numDims &&
+ "μ has to have the same number of dimensions as the generators!");
+
+ // First, the exponent in the numerator becomes
+ // - (μ • u_1) * (floor(first col of v))
+ // - (μ • u_2) * (floor(second col of v)) - ...
+ // - (μ • u_d) * (floor(d'th col of v))
+ // So we store the negation of the dot products.
+
+ // We have d terms, each of whose coefficient is the negative dot product,
+ SmallVector<Fraction> coefficients;
+ coefficients.reserve(numDims);
+ for (const Point &d : ds)
+ coefficients.push_back(-dotProduct(mu, d));
+
+ // Then, the affine function is a single floor expression, given by the
+ // corresponding column of v.
+ ParamPoint vTranspose = v.transpose();
+ std::vector<std::vector<SmallVector<Fraction>>> affine;
+ affine.reserve(numDims);
+ for (unsigned j = 0; j < numDims; ++j)
+ affine.push_back({SmallVector<Fraction>(vTranspose.getRow(j))});
+
+ QuasiPolynomial num(numParams, coefficients, affine);
+ num = num.simplify();
+
+ std::vector<Fraction> dens;
+ dens.reserve(ds.size());
+ // Similarly, each term in the denominator has exponent
+ // given by the dot product of μ with u_i.
+ for (const Point &d : ds) {
+ // This term in the denominator is
+ // (1 - t^dens.back())
+ dens.push_back(dotProduct(d, mu));
+ }
+
+ return {num, dens};
+}
+
+/// Normalize all denominator exponents `dens` to their absolute values
+/// by multiplying and dividing by the inverses.
+/// Also, take the factors (-s) out of each term of the product.
+void normalizeDenominatorExponents(int &sign, QuasiPolynomial &num,
+ std::vector<Fraction> &dens) {
+ // We track the number of exponents that are negative in the
+ // denominator, and convert them to their absolute values.
+ unsigned numNegExps = 0;
+ Fraction sumNegExps(0, 1);
+ for (unsigned j = 0, e = dens.size(); j < e; ++j) {
+ if (dens[j] < 0) {
+ numNegExps += 1;
+ sumNegExps += dens[j];
+ }
+ }
+
+ // If we have (1 - (s+1)^-c) in the denominator,
+ // multiply and divide by (s+1)^c.
+ // We convert all negative-exponent terms at once; therefore
+ // we multiply and divide by (s+1)^sumNegExps.
+ // Then we get
+ // -(1 - (s+1)^c) in the denominator,
+ // increase the numerator by c, and
+ // flip the sign.
+ if (numNegExps % 2 == 1)
+ sign = -sign;
+ num = num - QuasiPolynomial(num.getNumInputs(), sumNegExps);
+}
+
+/// Compute the binomial coefficients nCi for 0 ≤ i ≤ r,
+/// where n is a QuasiPolynomial.
+std::vector<QuasiPolynomial> getBinomialCoefficients(QuasiPolynomial n,
+ unsigned r) {
+ unsigned numParams = n.getNumInputs();
+ std::vector<QuasiPolynomial> coefficients;
+ coefficients.reserve(r + 1);
+ coefficients.push_back(QuasiPolynomial(numParams, 1));
+ for (unsigned j = 1; j <= r; ++j)
+ // We use the recursive formula for binomial coefficients here and below.
+ coefficients.push_back(
+ (coefficients[j - 1] * (n - QuasiPolynomial(numParams, j - 1)) /
+ Fraction(j, 1))
+ .simplify());
+ return coefficients;
+}
+
+/// Compute the binomial coefficients nCi for 0 ≤ i ≤ r,
+/// where n is a QuasiPolynomial.
+std::vector<Fraction> getBinomialCoefficients(Fraction n, Fraction r) {
+ std::vector<Fraction> coefficients;
+ coefficients.reserve((int64_t)floor(r));
+ coefficients.push_back(1);
+ for (unsigned j = 1; j <= r; ++j)
+ coefficients.push_back(coefficients[j - 1] * (n - (j - 1)) / (j));
+ return coefficients;
+}
+
+/// We have a generating function of the form
+/// f_p(x) = \sum_i sign_i * (x^n_i(p)) / (\prod_j (1 - x^d_{ij})
+///
+/// where sign_i is ±1,
+/// n_i \in Q^p -> Q^d is the sum of the vectors d_{ij}, weighted by the
+/// floors of d affine functions on p parameters.
+/// d_{ij} \in Q^d are vectors.
+///
+/// We need to find the number of terms of the form x^t in the expansion of
+/// this function.
+/// However, direct substitution (x = (1, ..., 1)) causes the denominator
+/// to become zero.
+///
+/// We therefore use the following procedure instead:
+/// 1. Substitute x_i = (s+1)^μ_i for some vector μ. This makes the generating
+/// function a function of a scalar s.
+/// 2. Write each term in this function as P(s)/Q(s), where P and Q are
+/// polynomials. P has coefficients as quasipolynomials in d parameters, while
+/// Q has coefficients as scalars.
+/// 3. Find the constant term in the expansion of each term P(s)/Q(s). This is
+/// equivalent to substituting s = 0.
+///
+/// Verdoolaege, Sven, et al. "Counting integer points in parametric
+/// polytopes using Barvinok's rational functions." Algorithmica 48 (2007):
+/// 37-66.
+QuasiPolynomial
+mlir::presburger::detail::computeNumTerms(const GeneratingFunction &gf) {
+ // Step (1) We need to find a μ such that we can substitute x_i =
+ // (s+1)^μ_i. After this substitution, the exponent of (s+1) in the
+ // denominator is (μ_i • d_{ij}) in each term. Clearly, this cannot become
+ // zero. Hence we find a vector μ that is not orthogonal to any of the
+ // d_{ij} and substitute x accordingly.
+ std::vector<Point> allDenominators;
+ for (ArrayRef<Point> den : gf.getDenominators())
+ allDenominators.insert(allDenominators.end(), den.begin(), den.end());
+ Point mu = getNonOrthogonalVector(allDenominators);
+
+ unsigned numParams = gf.getNumParams();
+ const std::vector<std::vector<Point>> &ds = gf.getDenominators();
+ QuasiPolynomial totalTerm(numParams, 0);
+ for (unsigned i = 0, e = ds.size(); i < e; ++i) {
+ int sign = gf.getSigns()[i];
+
+ // Compute the new exponents of (s+1) for the numerator and the
+ // denominator after substituting μ.
+ auto [numExp, dens] =
+ substituteMuInTerm(numParams, gf.getNumerators()[i], ds[i], mu);
+ // Now the numerator is (s+1)^numExp
+ // and the denominator is \prod_j (1 - (s+1)^dens[j]).
+
+ // Step (2) We need to express the terms in the function as quotients of
+ // polynomials. Each term is now of the form
+ // sign_i * (s+1)^numExp / (\prod_j (1 - (s+1)^dens[j]))
+ // For the i'th term, we first normalize the denominator to have only
+ // positive exponents. We convert all the dens[j] to their
+ // absolute values and change the sign and exponent in the numerator.
+ normalizeDenominatorExponents(sign, numExp, dens);
+
+ // Then, using the formula for geometric series, we replace each (1 -
+ // (s+1)^(dens[j])) with
+ // (-s)(\sum_{0 ≤ k < dens[j]} (s+1)^k).
+ for (unsigned j = 0, e = dens.size(); j < e; ++j)
+ dens[j] = abs(dens[j]) - 1;
+ // Note that at this point, the semantics of `dens[j]` changes to mean
+ // the limit of summation of (s+1)^t terms, and not the actual exponent.
+
+ // Since the -s are taken out, the sign changes if there is an odd number
+ // of such terms.
+ unsigned r = dens.size();
+ if (dens.size() % 2 == 1)
+ sign = -sign;
+
+ // Thus the term overall has now the form
+ // sign'_i * (s+1)^numExp /
+ // (s^r * \prod_j (\sum_{0 ≤ k < dens[j]} (s+1)^k)).
+ // This means that
+ // the numerator is a polynomial in s, with coefficients as
+ // quasipolynomials (given by binomial coefficients), and the denominator
+ // is polynomial in s, with integral coefficients (given by taking the
----------------
Superty wrote:
is a polynomial
https://github.com/llvm/llvm-project/pull/78078
More information about the Mlir-commits
mailing list