[Mlir-commits] [mlir] [MLIR][Presburger] Helper functions to compute the constant term of a generating function (PR #77819)

llvmlistbot at llvm.org llvmlistbot at llvm.org
Thu Jan 11 11:15:47 PST 2024


llvmbot wrote:


<!--LLVM PR SUMMARY COMMENT-->

@llvm/pr-subscribers-mlir-presburger

Author: None (Abhinav271828)

<details>
<summary>Changes</summary>

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).

---
Full diff: https://github.com/llvm/llvm-project/pull/77819.diff


6 Files Affected:

- (modified) mlir/include/mlir/Analysis/Presburger/Barvinok.h (+14) 
- (modified) mlir/include/mlir/Analysis/Presburger/QuasiPolynomial.h (+2) 
- (modified) mlir/lib/Analysis/Presburger/Barvinok.cpp (+99) 
- (modified) mlir/lib/Analysis/Presburger/QuasiPolynomial.cpp (+6) 
- (modified) mlir/unittests/Analysis/Presburger/BarvinokTest.cpp (+56) 
- (modified) mlir/unittests/Analysis/Presburger/Utils.h (+14) 


``````````diff
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index 213af636e5964d..81259e0b9f87b9 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,19 @@ 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(std::vector<Point> vectors);
+
+/// Find the coefficient of a given power of s in a rational function
+/// given by P(s)/Q(s), where the coefficients in P are QuasiPolynomials,
+/// and those in Q are Fractions.
+QuasiPolynomial getCoefficientInRationalFunction(unsigned power,
+                                                 std::vector<QuasiPolynomial>,
+                                                 std::vector<Fraction>);
+
 } // 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..e656c7e3fb00c1 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 {
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index 0bdc9015c3d647..7572267c87a8ee 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -144,3 +144,102 @@ GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
                             std::vector({numerator}),
                             std::vector({denominator}));
 }
+
+/// We use a recursive procedure to find a vector not orthogonal
+/// to a given set. 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 v and u,
+/// vs ++ [v] means the vector vs with the new element v appended to it.
+///
+/// Suppose we have a vector vs which is not orthogonal to
+/// any of {x_1[:n-1], ..., x_k[:n-1]}.
+/// Then we need v s.t. <x_i, vs++[v]> != 0 for all i.
+/// => <x_i[:n-1], vs> + x_i[-1]*v != 0
+/// => v != - <x_i[:n-1], vs> / x_i[-1]
+/// We compute this value for all i, and then
+/// set v to be the maximum element of this set + 1. Thus
+/// v is outside the set as desired, and we append it to vs.
+///
+/// The base case is given in one dimension,
+/// where the vector [1] is not orthogonal to any
+/// of the input vectors (since they are all nonzero).
+Point mlir::presburger::detail::getNonOrthogonalVector(
+    std::vector<Point> vectors) {
+  unsigned dim = vectors[0].size();
+
+  SmallVector<Fraction> newPoint = {Fraction(1, 1)};
+  std::vector<Fraction> lowerDimDotProducts;
+  Fraction dotProduct;
+  Fraction maxDisallowedValue = Fraction(-1, 0),
+           disallowedValue = Fraction(0, 1);
+  Fraction newValue;
+
+  for (unsigned d = 2; d <= dim; d++) {
+    lowerDimDotProducts.clear();
+
+    // Compute the set of dot products <x_i[:d-1], vs> for each i.
+    for (const Point &vector : vectors) {
+      dotProduct = Fraction(0, 1);
+      for (unsigned i = 0; i < d - 1; i++)
+        dotProduct = dotProduct + vector[i] * newPoint[i];
+      lowerDimDotProducts.push_back(dotProduct);
+    }
+
+    // Compute - <x_i[:n-1], vs> / x_i[-1] for each i,
+    // and find the biggest such value.
+    for (unsigned i = 0, e = vectors.size(); i < e; ++i) {
+      if (vectors[i][d - 1] == 0)
+        continue;
+      disallowedValue = -lowerDimDotProducts[i] / vectors[i][d - 1];
+      if (maxDisallowedValue < disallowedValue)
+        maxDisallowedValue = disallowedValue;
+    }
+
+    newValue = Fraction(ceil(maxDisallowedValue + Fraction(1, 1)), 1);
+    newPoint.append(1, newValue);
+  }
+  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 `getCoefficientInRationalFuncion` on
+/// all i \in [0, power).
+///
+/// https://math.ucdavis.edu/~deloera/researchsummary/
+/// barvinokalgorithm-latte1.pdf, p. 1285
+QuasiPolynomial mlir::presburger::detail::getCoefficientInRationalFunction(
+    unsigned power, std::vector<QuasiPolynomial> num,
+    std::vector<Fraction> den) {
+
+  unsigned numParam = num[0].getNumInputs();
+  unsigned limit;
+
+  std::vector<QuasiPolynomial> coefficients(power + 1,
+                                            QuasiPolynomial(numParam, 0));
+
+  coefficients[0] = num[0] / den[0];
+
+  for (unsigned i = 1; i <= power; i++) {
+    if (i < num.size())
+      coefficients[i] = num[i];
+
+    limit = i + 1 < den.size() ? i + 1 : den.size();
+
+    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..24878dd47a2acf 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 different numbers of symbols cannot "
diff --git a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
index 2936d95c802e9c..072a24262604cb 100644
--- a/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/BarvinokTest.cpp
@@ -82,3 +82,59 @@ 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 nonzero
+// 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_FALSE(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_FALSE(dotProduct(nonOrth, vector) == 0);
+}
+
+Fraction getC(QuasiPolynomial q) {
+  Fraction t = 0;
+  for (const Fraction &c : q.getCoefficients())
+    t += c;
+
+  return t;
+}
+
+// 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_CONSTANT_TERM_QUASIPOLYNOMIAL(coeff, 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_CONSTANT_TERM_QUASIPOLYNOMIAL(coeff, Fraction(55, 64));
+}
diff --git a/mlir/unittests/Analysis/Presburger/Utils.h b/mlir/unittests/Analysis/Presburger/Utils.h
index 6b00898a7e2749..2c05fe7e2057f9 100644
--- a/mlir/unittests/Analysis/Presburger/Utils.h
+++ b/mlir/unittests/Analysis/Presburger/Utils.h
@@ -128,6 +128,20 @@ inline void EXPECT_EQ_REPR_QUASIPOLYNOMIAL(QuasiPolynomial a,
   }
 }
 
+// Check the constant term of the quasipolynomial against a constant.
+inline void EXPECT_CONSTANT_TERM_QUASIPOLYNOMIAL(QuasiPolynomial qp,
+                                                 Fraction f) {
+  SmallVector<Fraction> coeffs = qp.getCoefficients();
+  std::vector<std::vector<SmallVector<Fraction>>> aff = qp.getAffine();
+
+  Fraction t = 0;
+  for (unsigned i = 0, e = coeffs.size(); i < e; ++i) {
+    if (aff[i].size() == 0)
+      t += coeffs[i];
+  }
+  EXPECT_EQ(t, f);
+}
+
 /// lhs and rhs represent non-negative integers or positive infinity. The
 /// infinity case corresponds to when the Optional is empty.
 inline bool infinityOrUInt64LE(std::optional<MPInt> lhs,

``````````

</details>


https://github.com/llvm/llvm-project/pull/77819


More information about the Mlir-commits mailing list