[Mlir-commits] [mlir] 850f713 - [MLIR][Presburger] Helper functions to compute the constant term of a generating function (#77819)
llvmlistbot at llvm.org
llvmlistbot at llvm.org
Sat Jan 13 08:00:12 PST 2024
Author: Abhinav271828
Date: 2024-01-13T21:30:06+05:30
New Revision: 850f713e80426f1706c0d3dad143c330ca872d5d
URL: https://github.com/llvm/llvm-project/commit/850f713e80426f1706c0d3dad143c330ca872d5d
DIFF: https://github.com/llvm/llvm-project/commit/850f713e80426f1706c0d3dad143c330ca872d5d.diff
LOG: [MLIR][Presburger] Helper functions to compute the constant term of a generating function (#77819)
We implement two functions that are needed to compute the constant term
of a GF.
One finds a vector not orthogonal to all the non-null vectors in a given
set.
One computes the coefficient of any term in an arbitrary rational
function (quotient of two polynomials).
Added:
Modified:
mlir/include/mlir/Analysis/Presburger/Barvinok.h
mlir/include/mlir/Analysis/Presburger/QuasiPolynomial.h
mlir/lib/Analysis/Presburger/Barvinok.cpp
mlir/lib/Analysis/Presburger/QuasiPolynomial.cpp
mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
Removed:
################################################################################
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 213af636e5964d..edee19f0e1a535 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -27,6 +27,7 @@
#include "mlir/Analysis/Presburger/GeneratingFunction.h"
#include "mlir/Analysis/Presburger/IntegerRelation.h"
#include "mlir/Analysis/Presburger/Matrix.h"
+#include "mlir/Analysis/Presburger/QuasiPolynomial.h"
#include <optional>
namespace mlir {
@@ -83,6 +84,21 @@ ConeH getDual(ConeV cone);
GeneratingFunction unimodularConeGeneratingFunction(ParamPoint vertex, int sign,
ConeH cone);
+/// Find a vector that is not orthogonal to any of the given vectors,
+/// i.e., has nonzero dot product with those of the given vectors
+/// that are not null.
+/// If any of the vectors is null, it is ignored.
+Point getNonOrthogonalVector(ArrayRef<Point> vectors);
+
+/// Find the coefficient of a given power of s in a rational function
+/// given by P(s)/Q(s), where
+/// P is a polynomial, in which the coefficients are QuasiPolynomials
+/// over d parameters (distinct from s), and
+/// and Q is a polynomial with Fraction coefficients.
+QuasiPolynomial getCoefficientInRationalFunction(unsigned power,
+ ArrayRef<QuasiPolynomial> num,
+ ArrayRef<Fraction> den);
+
} // namespace detail
} // namespace presburger
} // namespace mlir
diff --git a/mlir/include/mlir/Analysis/Presburger/QuasiPolynomial.h b/mlir/include/mlir/Analysis/Presburger/QuasiPolynomial.h
index f8ce8524e41b21..d03446f50264fa 100644
--- a/mlir/include/mlir/Analysis/Presburger/QuasiPolynomial.h
+++ b/mlir/include/mlir/Analysis/Presburger/QuasiPolynomial.h
@@ -39,6 +39,8 @@ class QuasiPolynomial : public PresburgerSpace {
QuasiPolynomial(unsigned numVars, SmallVector<Fraction> coeffs = {},
std::vector<std::vector<SmallVector<Fraction>>> aff = {});
+ QuasiPolynomial(unsigned numVars, Fraction constant);
+
// Find the number of inputs (numDomain) to the polynomial.
// numSymbols is set to zero.
unsigned getNumInputs() const {
@@ -60,6 +62,8 @@ class QuasiPolynomial : public PresburgerSpace {
// Removes terms which evaluate to zero from the expression.
QuasiPolynomial simplify();
+ Fraction getConstantTerm();
+
private:
SmallVector<Fraction> coefficients;
std::vector<std::vector<SmallVector<Fraction>>> affine;
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 0bdc9015c3d647..a3d210b3b8fd16 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -8,6 +8,7 @@
#include "mlir/Analysis/Presburger/Barvinok.h"
#include "llvm/ADT/Sequence.h"
+#include <algorithm>
using namespace mlir;
using namespace presburger;
@@ -144,3 +145,100 @@ GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
std::vector({numerator}),
std::vector({denominator}));
}
+
+/// We use an iterative procedure to find a vector not orthogonal
+/// to a given set, ignoring the null vectors.
+/// Let the inputs be {x_1, ..., x_k}, all vectors of length n.
+///
+/// In the following,
+/// vs[:i] means the elements of vs up to and including the i'th one,
+/// <vs, us> means the dot product of vs and us,
+/// vs ++ [v] means the vector vs with the new element v appended to it.
+///
+/// We proceed iteratively; for steps d = 0, ... n-1, we construct a vector
+/// which is not orthogonal to any of {x_1[:d], ..., x_n[:d]}, ignoring
+/// the null vectors.
+/// At step d = 0, we let vs = [1]. Clearly this is not orthogonal to
+/// any vector in the set {x_1[0], ..., x_n[0]}, except the null ones,
+/// which we ignore.
+/// At step d > 0 , we need a number v
+/// s.t. <x_i[:d], vs++[v]> != 0 for all i.
+/// => <x_i[:d-1], vs> + x_i[d]*v != 0
+/// => v != - <x_i[:d-1], vs> / x_i[d]
+/// We compute this value for all x_i, and then
+/// set v to be the maximum element of this set plus one. Thus
+/// v is outside the set as desired, and we append it to vs
+/// to obtain the result of the d'th step.
+Point mlir::presburger::detail::getNonOrthogonalVector(
+ ArrayRef<Point> vectors) {
+ unsigned dim = vectors[0].size();
+ for (const Point &vector : vectors)
+ assert(vector.size() == dim && "all vectors need to be the same size!");
+
+ SmallVector<Fraction> newPoint = {Fraction(1, 1)};
+ Fraction maxDisallowedValue = -Fraction(1, 0),
+ disallowedValue = Fraction(0, 1);
+
+ for (unsigned d = 1; d < dim; ++d) {
+ // Compute the disallowed values - <x_i[:d-1], vs> / x_i[d] for each i.
+ maxDisallowedValue = -Fraction(1, 0);
+ for (const Point &vector : vectors) {
+ if (vector[d] == 0)
+ continue;
+ disallowedValue =
+ -dotProduct(ArrayRef(vector).slice(0, d), newPoint) / vector[d];
+
+ // Find the biggest such value
+ maxDisallowedValue = std::max(maxDisallowedValue, disallowedValue);
+ }
+ newPoint.push_back(maxDisallowedValue + 1);
+ }
+ return newPoint;
+}
+
+/// We use the following recursive formula to find the coefficient of
+/// s^power in the rational function given by P(s)/Q(s).
+///
+/// Let P[i] denote the coefficient of s^i in the polynomial P(s).
+/// (P/Q)[r] =
+/// if (r == 0) then
+/// P[0]/Q[0]
+/// else
+/// (P[r] - {Σ_{i=1}^r (P/Q)[r-i] * Q[i])}/(Q[0])
+/// We therefore recursively call `getCoefficientInRationalFunction` on
+/// all i \in [0, power).
+///
+/// https://math.ucdavis.edu/~deloera/researchsummary/
+/// barvinokalgorithm-latte1.pdf, p. 1285
+QuasiPolynomial mlir::presburger::detail::getCoefficientInRationalFunction(
+ unsigned power, ArrayRef<QuasiPolynomial> num, ArrayRef<Fraction> den) {
+ assert(den.size() != 0 &&
+ "division by empty denominator in rational function!");
+
+ unsigned numParam = num[0].getNumInputs();
+ for (const QuasiPolynomial &qp : num)
+ // We use the `isEqual` method of PresburgerSpace, which QuasiPolynomial
+ // inherits from.
+ assert(num[0].isEqual(qp) &&
+ "the quasipolynomials should all belong to the same space!");
+
+ std::vector<QuasiPolynomial> coefficients;
+ coefficients.reserve(power + 1);
+
+ coefficients.push_back(num[0] / den[0]);
+ for (unsigned i = 1; i <= power; ++i) {
+ // If the power is not there in the numerator, the coefficient is zero.
+ coefficients.push_back(i < num.size() ? num[i]
+ : QuasiPolynomial(numParam, 0));
+
+ // After den.size(), the coefficients are zero, so we stop
+ // subtracting at that point (if it is less than i).
+ unsigned limit = std::min<unsigned long>(i, den.size() - 1);
+ for (unsigned j = 1; j <= limit; ++j)
+ coefficients[i] = coefficients[i] -
+ coefficients[i - j] * QuasiPolynomial(numParam, den[j]);
+
+ coefficients[i] = coefficients[i] / den[0];
+ }
+ return coefficients[power].simplify();
+}
\ No newline at end of file
diff --git a/mlir/lib/Analysis/Presburger/QuasiPolynomial.cpp b/mlir/lib/Analysis/Presburger/QuasiPolynomial.cpp
index feed683a203cff..385d4c354be18a 100644
--- a/mlir/lib/Analysis/Presburger/QuasiPolynomial.cpp
+++ b/mlir/lib/Analysis/Presburger/QuasiPolynomial.cpp
@@ -36,6 +36,12 @@ QuasiPolynomial::QuasiPolynomial(
#endif // NDEBUG
}
+/// Define a quasipolynomial which is a single constant.
+QuasiPolynomial::QuasiPolynomial(unsigned numVars, Fraction constant)
+ : PresburgerSpace(/*numDomain=*/numVars, /*numRange=*/1, /*numSymbols=*/0,
+ /*numLocals=*/0),
+ coefficients({constant}), affine({{}}) {}
+
QuasiPolynomial QuasiPolynomial::operator+(const QuasiPolynomial &x) const {
assert(getNumInputs() == x.getNumInputs() &&
"two quasi-polynomials with
diff erent numbers of symbols cannot "
@@ -113,3 +119,11 @@ QuasiPolynomial QuasiPolynomial::simplify() {
}
return QuasiPolynomial(getNumInputs(), newCoeffs, newAffine);
}
+
+Fraction QuasiPolynomial::getConstantTerm() {
+ Fraction constTerm = 0;
+ for (unsigned i = 0, e = coefficients.size(); i < e; ++i)
+ if (affine[i].size() == 0)
+ constTerm += coefficients[i];
+ return constTerm;
+}
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 2936d95c802e9c..e304f81de21f0b 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -82,3 +82,45 @@ TEST(BarvinokTest, unimodularConeGeneratingFunction) {
1, {1}, {makeFracMatrix(2, 3, {{-83, -100, -41}, {-22, -27, -15}})},
{{{8, 47, -17}, {-7, -41, 15}, {1, 5, -2}}}));
}
+
+// The following vectors are randomly generated.
+// We then check that the output of the function has non-zero
+// dot product with all non-null vectors.
+TEST(BarvinokTest, getNonOrthogonalVector) {
+ std::vector<Point> vectors = {Point({1, 2, 3, 4}), Point({-1, 0, 1, 1}),
+ Point({2, 7, 0, 0}), Point({0, 0, 0, 0})};
+ Point nonOrth = getNonOrthogonalVector(vectors);
+
+ for (unsigned i = 0; i < 3; i++)
+ EXPECT_NE(dotProduct(nonOrth, vectors[i]), 0);
+
+ vectors = {Point({0, 1, 3}), Point({-2, -1, 1}), Point({6, 3, 0}),
+ Point({0, 0, -3}), Point({5, 0, -1})};
+ nonOrth = getNonOrthogonalVector(vectors);
+
+ for (const Point &vector : vectors)
+ EXPECT_NE(dotProduct(nonOrth, vector), 0);
+}
+
+// The following polynomials are randomly generated and the
+// coefficients are computed by hand.
+// Although the function allows the coefficients of the numerator
+// to be arbitrary quasipolynomials, we stick to constants for simplicity,
+// as the relevant arithmetic operations on quasipolynomials
+// are tested separately.
+TEST(BarvinokTest, getCoefficientInRationalFunction) {
+ std::vector<QuasiPolynomial> numerator = {
+ QuasiPolynomial(0, 2), QuasiPolynomial(0, 3), QuasiPolynomial(0, 5)};
+ std::vector<Fraction> denominator = {Fraction(1), Fraction(0), Fraction(4),
+ Fraction(3)};
+ QuasiPolynomial coeff =
+ getCoefficientInRationalFunction(1, numerator, denominator);
+ EXPECT_EQ(coeff.getConstantTerm(), 3);
+
+ numerator = {QuasiPolynomial(0, -1), QuasiPolynomial(0, 4),
+ QuasiPolynomial(0, -2), QuasiPolynomial(0, 5),
+ QuasiPolynomial(0, 6)};
+ denominator = {Fraction(8), Fraction(4), Fraction(0), Fraction(-2)};
+ coeff = getCoefficientInRationalFunction(3, numerator, denominator);
+ EXPECT_EQ(coeff.getConstantTerm(), Fraction(55, 64));
+}
More information about the Mlir-commits
mailing list