[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
Sat Jan 20 08:56:12 PST 2024
================
@@ -245,3 +245,270 @@ QuasiPolynomial mlir::presburger::detail::getCoefficientInRationalFunction(
}
return coefficients[power].simplify();
}
+
+static std::vector<Fraction> convolution(ArrayRef<Fraction> a,
+ ArrayRef<Fraction> b) {
+ // The length of the convolution is the sum of the lengths
+ // of the two sequences. We pad the shorter one with zeroes.
+ unsigned len = a.size() + b.size();
+
+ // We define accessors to avoid out-of-bounds errors.
+ std::function<Fraction(unsigned i)> aGetItem = [a](unsigned i) -> Fraction {
+ if (i < a.size())
+ return a[i];
+ else
+ return 0;
+ };
+ std::function<Fraction(unsigned i)> bGetItem = [b](unsigned i) -> Fraction {
+ if (i < b.size())
+ return b[i];
+ else
+ return 0;
+ };
+
+ std::vector<Fraction> convolution;
+ convolution.reserve(len);
+ for (unsigned k = 0; k < len; ++k) {
+ Fraction sum(0, 1);
+ for (unsigned l = 0; l <= k; ++l)
+ sum += aGetItem(l) * bGetItem(k - l);
+ convolution.push_back(sum);
+ }
+ return convolution;
+}
+
+/// 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
+/// 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];
+ }
+ // All exponents will be made positive; then reduce
+ // (1 - (s+1)^x)
+ // to
+ // (-s)*(Σ_{x-1} (s+1)^j) because x > 0
+ dens[j] = abs(dens[j]) - 1;
+ }
+
+ // 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);
+
+ // Take all the (-s) out.
+ unsigned r = dens.size();
+ if (r % 2 == 1)
+ sign = -sign;
+}
+
+/// 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.
+/// 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 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 numerator and 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 by multiplying and dividing by (s+1)^(-dens[j]) if it is
+ // negative. We change the sign accordingly to keep the denominator in the
----------------
Superty wrote:
these are some implementation details that have been abstracted away into the function below, hence the user doesn't need to read them here. they can go to the function for that
https://github.com/llvm/llvm-project/pull/78078
More information about the Mlir-commits
mailing list