[Mlir-commits] [mlir] cfd51fb - [MLIR][Presburger] Add LLL basis reduction (#75565)
llvmlistbot at llvm.org
llvmlistbot at llvm.org
Tue Dec 19 08:31:43 PST 2023
Author: Abhinav271828
Date: 2023-12-19T17:31:38+01:00
New Revision: cfd51fbadd19a7b83805a84cb5745f1350c798c4
URL: https://github.com/llvm/llvm-project/commit/cfd51fbadd19a7b83805a84cb5745f1350c798c4
DIFF: https://github.com/llvm/llvm-project/commit/cfd51fbadd19a7b83805a84cb5745f1350c798c4.diff
LOG: [MLIR][Presburger] Add LLL basis reduction (#75565)
Add a method for LLL basis reduction to the FracMatrix class.
This needs an abs() method for Fractions, which is added to Fraction.h.
Added:
Modified:
mlir/include/mlir/Analysis/Presburger/Fraction.h
mlir/include/mlir/Analysis/Presburger/Matrix.h
mlir/lib/Analysis/Presburger/Matrix.cpp
mlir/unittests/Analysis/Presburger/MatrixTest.cpp
Removed:
################################################################################
diff --git a/mlir/include/mlir/Analysis/Presburger/Fraction.h b/mlir/include/mlir/Analysis/Presburger/Fraction.h
index afcbed84c66bc3..c07bb767f50bfc 100644
--- a/mlir/include/mlir/Analysis/Presburger/Fraction.h
+++ b/mlir/include/mlir/Analysis/Presburger/Fraction.h
@@ -101,6 +101,11 @@ inline bool operator>=(const Fraction &x, const Fraction &y) {
return compare(x, y) >= 0;
}
+inline Fraction abs(const Fraction &f) {
+ assert(f.den > 0 && "denominator of fraction must be positive!");
+ return Fraction(abs(f.num), f.den);
+}
+
inline Fraction reduce(const Fraction &f) {
if (f == Fraction(0))
return Fraction(0, 1);
@@ -124,6 +129,12 @@ inline Fraction operator-(const Fraction &x, const Fraction &y) {
return reduce(Fraction(x.num * y.den - x.den * y.num, x.den * y.den));
}
+// Find the integer nearest to a given fraction.
+inline MPInt round(const Fraction &f) {
+ MPInt rem = f.num % f.den;
+ return (f.num / f.den) + (rem > f.den / 2);
+}
+
inline Fraction &operator+=(Fraction &x, const Fraction &y) {
x = x + y;
return x;
diff --git a/mlir/include/mlir/Analysis/Presburger/Matrix.h b/mlir/include/mlir/Analysis/Presburger/Matrix.h
index 89fad85c0c3374..347e2e0489786f 100644
--- a/mlir/include/mlir/Analysis/Presburger/Matrix.h
+++ b/mlir/include/mlir/Analysis/Presburger/Matrix.h
@@ -270,6 +270,12 @@ class FracMatrix : public Matrix<Fraction> {
// of the rows of matrix (cubic time).
// The rows of the matrix must be linearly independent.
FracMatrix gramSchmidt() const;
+
+ // Run LLL basis reduction on the matrix, modifying it in-place.
+ // The parameter is what [the original
+ // paper](https://www.cs.cmu.edu/~avrim/451f11/lectures/lect1129_LLL.pdf)
+ // calls `y`, usually 3/4.
+ void LLL(Fraction delta);
};
} // namespace presburger
diff --git a/mlir/lib/Analysis/Presburger/Matrix.cpp b/mlir/lib/Analysis/Presburger/Matrix.cpp
index 1fcc6d072b44b7..25300f84cfc047 100644
--- a/mlir/lib/Analysis/Presburger/Matrix.cpp
+++ b/mlir/lib/Analysis/Presburger/Matrix.cpp
@@ -576,4 +576,71 @@ FracMatrix FracMatrix::gramSchmidt() const {
}
}
return orth;
-}
\ No newline at end of file
+}
+
+// Convert the matrix, interpreted (row-wise) as a basis
+// to an LLL-reduced basis.
+//
+// This is an implementation of the algorithm described in
+// "Factoring polynomials with rational coefficients" by
+// A. K. Lenstra, H. W. Lenstra Jr., L. Lovasz.
+//
+// Let {b_1, ..., b_n} be the current basis and
+// {b_1*, ..., b_n*} be the Gram-Schmidt orthogonalised
+// basis (unnormalized).
+// Define the Gram-Schmidt coefficients μ_ij as
+// (b_i • b_j*) / (b_j* • b_j*), where (•) represents the inner product.
+//
+// We iterate starting from the second row to the last row.
+//
+// For the kth row, we first check μ_kj for all rows j < k.
+// We subtract b_j (scaled by the integer nearest to μ_kj)
+// from b_k.
+//
+// Now, we update k.
+// If b_k and b_{k-1} satisfy the Lovasz condition
+// |b_k|^2 ≥ (δ - μ_k{k-1}^2) |b_{k-1}|^2,
+// we are done and we increment k.
+// Otherwise, we swap b_k and b_{k-1} and decrement k.
+//
+// We repeat this until k = n and return.
+void FracMatrix::LLL(Fraction delta) {
+ MPInt nearest;
+ Fraction mu;
+
+ // `gsOrth` holds the Gram-Schmidt orthogonalisation
+ // of the matrix at all times. It is recomputed every
+ // time the matrix is modified during the algorithm.
+ // This is naive and can be optimised.
+ FracMatrix gsOrth = gramSchmidt();
+
+ // We start from the second row.
+ unsigned k = 1;
+ while (k < getNumRows()) {
+ for (unsigned j = k - 1; j < k; j--) {
+ // Compute the Gram-Schmidt coefficient μ_jk.
+ mu = dotProduct(getRow(k), gsOrth.getRow(j)) /
+ dotProduct(gsOrth.getRow(j), gsOrth.getRow(j));
+ nearest = round(mu);
+ // Subtract b_j scaled by the integer nearest to μ_jk from b_k.
+ addToRow(k, getRow(j), -Fraction(nearest, 1));
+ gsOrth = gramSchmidt(); // Update orthogonalization.
+ }
+ mu = dotProduct(getRow(k), gsOrth.getRow(k - 1)) /
+ dotProduct(gsOrth.getRow(k - 1), gsOrth.getRow(k - 1));
+ // Check the Lovasz condition for b_k and b_{k-1}.
+ if (dotProduct(gsOrth.getRow(k), gsOrth.getRow(k)) >
+ (delta - mu * mu) *
+ dotProduct(gsOrth.getRow(k - 1), gsOrth.getRow(k - 1))) {
+ // If it is satisfied, proceed to the next k.
+ k += 1;
+ } else {
+ // If it is not satisfied, decrement k (without
+ // going beyond the second row).
+ swapRows(k, k - 1);
+ gsOrth = gramSchmidt(); // Update orthogonalization.
+ k = k > 1 ? k - 1 : 1;
+ }
+ }
+ return;
+}
diff --git a/mlir/unittests/Analysis/Presburger/MatrixTest.cpp b/mlir/unittests/Analysis/Presburger/MatrixTest.cpp
index 508d4fa369c14c..e6e452790f82d1 100644
--- a/mlir/unittests/Analysis/Presburger/MatrixTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/MatrixTest.cpp
@@ -377,4 +377,52 @@ TEST(MatrixTest, gramSchmidt) {
gs = mat.gramSchmidt();
EXPECT_EQ_FRAC_MATRIX(gs, FracMatrix::identity(10));
+}
+
+void checkReducedBasis(FracMatrix mat, Fraction delta) {
+ FracMatrix gsOrth = mat.gramSchmidt();
+
+ // Size-reduced check.
+ for (unsigned i = 0, e = mat.getNumRows(); i < e; i++) {
+ for (unsigned j = 0; j < i; j++) {
+ Fraction mu = dotProduct(mat.getRow(i), gsOrth.getRow(j)) /
+ dotProduct(gsOrth.getRow(j), gsOrth.getRow(j));
+ EXPECT_TRUE(abs(mu) <= Fraction(1, 2));
+ }
+ }
+
+ // Lovasz condition check.
+ for (unsigned i = 1, e = mat.getNumRows(); i < e; i++) {
+ Fraction mu = dotProduct(mat.getRow(i), gsOrth.getRow(i - 1)) /
+ dotProduct(gsOrth.getRow(i - 1), gsOrth.getRow(i - 1));
+ EXPECT_TRUE(dotProduct(mat.getRow(i), mat.getRow(i)) >
+ (delta - mu * mu) *
+ dotProduct(gsOrth.getRow(i - 1), gsOrth.getRow(i - 1)));
+ }
+}
+
+TEST(MatrixTest, LLL) {
+ FracMatrix mat =
+ makeFracMatrix(3, 3,
+ {{Fraction(1, 1), Fraction(1, 1), Fraction(1, 1)},
+ {Fraction(-1, 1), Fraction(0, 1), Fraction(2, 1)},
+ {Fraction(3, 1), Fraction(5, 1), Fraction(6, 1)}});
+ mat.LLL(Fraction(3, 4));
+
+ checkReducedBasis(mat, Fraction(3, 4));
+
+ mat = makeFracMatrix(
+ 2, 2,
+ {{Fraction(12, 1), Fraction(2, 1)}, {Fraction(13, 1), Fraction(4, 1)}});
+ mat.LLL(Fraction(3, 4));
+
+ checkReducedBasis(mat, Fraction(3, 4));
+
+ mat = makeFracMatrix(3, 3,
+ {{Fraction(1, 1), Fraction(0, 1), Fraction(2, 1)},
+ {Fraction(0, 1), Fraction(1, 3), -Fraction(5, 3)},
+ {Fraction(0, 1), Fraction(0, 1), Fraction(1, 1)}});
+ mat.LLL(Fraction(3, 4));
+
+ checkReducedBasis(mat, Fraction(3, 4));
}
\ No newline at end of file
More information about the Mlir-commits
mailing list