[Mlir-commits] [mlir] 49b754b - [MLIR][NFC] Move presburger functionality from FlatAffineConstraints to IntegerPolyhedron

llvmlistbot at llvm.org llvmlistbot at llvm.org
Fri Jan 7 11:59:21 PST 2022


Author: Groverkss
Date: 2022-01-08T01:22:49+05:30
New Revision: 49b754b5c688e0b819c2422abb10f77775042e91

URL: https://github.com/llvm/llvm-project/commit/49b754b5c688e0b819c2422abb10f77775042e91
DIFF: https://github.com/llvm/llvm-project/commit/49b754b5c688e0b819c2422abb10f77775042e91.diff

LOG: [MLIR][NFC] Move presburger functionality from FlatAffineConstraints to IntegerPolyhedron

This patch moves all presburger functionality from FlatAffineConstraints to
IntegerPolyhedron. This patch is purely mechanical, it only moves and renames
functionality and tests.

This patch is part of a series of patches to move presburger functionality to
Presburger/ directory.

Reviewed By: bondhugula

Differential Revision: https://reviews.llvm.org/D116681

Added: 
    

Modified: 
    mlir/include/mlir/Analysis/AffineStructures.h
    mlir/include/mlir/Analysis/Presburger/IntegerPolyhedron.h
    mlir/lib/Analysis/AffineStructures.cpp
    mlir/lib/Analysis/Presburger/IntegerPolyhedron.cpp
    mlir/unittests/Analysis/CMakeLists.txt
    mlir/unittests/Analysis/Presburger/IntegerPolyhedronTest.cpp

Removed: 
    mlir/unittests/Analysis/AffineStructuresTest.cpp


################################################################################
diff  --git a/mlir/include/mlir/Analysis/AffineStructures.h b/mlir/include/mlir/Analysis/AffineStructures.h
index 3991628a6100..ed4ca4cea96a 100644
--- a/mlir/include/mlir/Analysis/AffineStructures.h
+++ b/mlir/include/mlir/Analysis/AffineStructures.h
@@ -102,66 +102,9 @@ class FlatAffineConstraints : public IntegerPolyhedron {
     return cst->getKind() == Kind::FlatAffineConstraints;
   }
 
-  /// Checks for emptiness by performing variable elimination on all
-  /// identifiers, running the GCD test on each equality constraint, and
-  /// checking for invalid constraints. Returns true if the GCD test fails for
-  /// any equality, or if any invalid constraints are discovered on any row.
-  /// Returns false otherwise.
-  bool isEmpty() const;
-
-  /// Runs the GCD test on all equality constraints. Returns true if this test
-  /// fails on any equality. Returns false otherwise.
-  /// This test can be used to disprove the existence of a solution. If it
-  /// returns true, no integer solution to the equality constraints can exist.
-  bool isEmptyByGCDTest() const;
-
-  /// Returns true if the set of constraints is found to have no solution,
-  /// false if a solution exists. Uses the same algorithm as
-  /// `findIntegerSample`.
-  bool isIntegerEmpty() const;
-
-  /// Returns a matrix where each row is a vector along which the polytope is
-  /// bounded. The span of the returned vectors is guaranteed to contain all
-  /// such vectors. The returned vectors are NOT guaranteed to be linearly
-  /// independent. This function should not be called on empty sets.
-  Matrix getBoundedDirections() const;
-
-  /// Find an integer sample point satisfying the constraints using a
-  /// branch and bound algorithm with generalized basis reduction, with some
-  /// additional processing using Simplex for unbounded sets.
-  ///
-  /// Returns an integer sample point if one exists, or an empty Optional
-  /// otherwise.
-  Optional<SmallVector<int64_t, 8>> findIntegerSample() const;
-
-  /// Returns true if the given point satisfies the constraints, or false
-  /// otherwise.
-  bool containsPoint(ArrayRef<int64_t> point) const;
-
-  /// Find pairs of inequalities identified by their position indices, using
-  /// which an explicit representation for each local variable can be computed.
-  /// The pairs are stored as indices of upperbound, lowerbound inequalities. If
-  /// no such pair can be found, it is stored as llvm::None.
-  ///
-  /// The dividends of the explicit representations are stored in `dividends`
-  /// and the denominators in `denominators`. If no explicit representation
-  /// could be found for the `i^th` local identifier, `denominators[i]` is set
-  /// to 0.
-  void getLocalReprs(
-      std::vector<SmallVector<int64_t, 8>> &dividends,
-      SmallVector<unsigned, 4> &denominators,
-      std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const;
-  void getLocalReprs(
-      std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const;
-  void getLocalReprs(std::vector<SmallVector<int64_t, 8>> &dividends,
-                     SmallVector<unsigned, 4> &denominators) const;
-
   // Clones this object.
   std::unique_ptr<FlatAffineConstraints> clone() const;
 
-  /// The type of bound: equal, lower bound or upper bound.
-  enum BoundType { EQ, LB, UB };
-
   /// Adds a bound for the identifier at the specified position with constraints
   /// being drawn from the specified bound map. In case of an EQ bound, the
   /// bound map is expected to have exactly one result. In case of a LB/UB, the
@@ -171,11 +114,9 @@ class FlatAffineConstraints : public IntegerPolyhedron {
   /// dimensions/symbols of the affine map.
   LogicalResult addBound(BoundType type, unsigned pos, AffineMap boundMap);
 
-  /// Adds a constant bound for the specified identifier.
-  void addBound(BoundType type, unsigned pos, int64_t value);
-
-  /// Adds a constant bound for the specified expression.
-  void addBound(BoundType type, ArrayRef<int64_t> expr, int64_t value);
+  /// The `addBound` overload above hides the inherited overloads by default, so
+  /// we explicitly introduce them here.
+  using IntegerPolyhedron::addBound;
 
   /// Returns the constraint system as an integer set. Returns a null integer
   /// set if the system has no constraints, or if an integer set couldn't be
@@ -193,91 +134,15 @@ class FlatAffineConstraints : public IntegerPolyhedron {
                       SmallVectorImpl<AffineMap> *lbMaps,
                       SmallVectorImpl<AffineMap> *ubMaps);
 
-  /// Adds a new local identifier as the floordiv of an affine function of other
-  /// identifiers, the coefficients of which are provided in `dividend` and with
-  /// respect to a positive constant `divisor`. Two constraints are added to the
-  /// system to capture equivalence with the floordiv:
-  /// q = dividend floordiv c    <=>   c*q <= dividend <= c*q + c - 1.
-  void addLocalFloorDiv(ArrayRef<int64_t> dividend, int64_t divisor);
-
   /// Composes an affine map whose dimensions and symbols match one to one with
   /// the dimensions and symbols of this FlatAffineConstraints. The results of
   /// the map `other` are added as the leading dimensions of this constraint
   /// system. Returns failure if `other` is a semi-affine map.
   LogicalResult composeMatchingMap(AffineMap other);
 
-  /// Projects out (aka eliminates) `num` identifiers starting at position
-  /// `pos`. The resulting constraint system is the shadow along the dimensions
-  /// that still exist. This method may not always be integer exact.
-  // TODO: deal with integer exactness when necessary - can return a value to
-  // mark exactness for example.
-  void projectOut(unsigned pos, unsigned num);
-  inline void projectOut(unsigned pos) { return projectOut(pos, 1); }
-
-  /// Changes the partition between dimensions and symbols. Depending on the new
-  /// symbol count, either a chunk of trailing dimensional identifiers becomes
-  /// symbols, or some of the leading symbols become dimensions.
-  void setDimSymbolSeparation(unsigned newSymbolCount);
-
-  /// Tries to fold the specified identifier to a constant using a trivial
-  /// equality detection; if successful, the constant is substituted for the
-  /// identifier everywhere in the constraint system and then removed from the
-  /// system.
-  LogicalResult constantFoldId(unsigned pos);
-
-  /// This method calls `constantFoldId` for the specified range of identifiers,
-  /// `num` identifiers starting at position `pos`.
-  void constantFoldIdRange(unsigned pos, unsigned num);
-
-  /// Updates the constraints to be the smallest bounding (enclosing) box that
-  /// contains the points of `this` set and that of `other`, with the symbols
-  /// being treated specially. For each of the dimensions, the min of the lower
-  /// bounds (symbolic) and the max of the upper bounds (symbolic) is computed
-  /// to determine such a bounding box. `other` is expected to have the same
-  /// dimensional identifiers as this constraint system (in the same order).
-  ///
-  /// E.g.:
-  /// 1) this   = {0 <= d0 <= 127},
-  ///    other  = {16 <= d0 <= 192},
-  ///    output = {0 <= d0 <= 192}
-  /// 2) this   = {s0 + 5 <= d0 <= s0 + 20},
-  ///    other  = {s0 + 1 <= d0 <= s0 + 9},
-  ///    output = {s0 + 1 <= d0 <= s0 + 20}
-  /// 3) this   = {0 <= d0 <= 5, 1 <= d1 <= 9}
-  ///    other  = {2 <= d0 <= 6, 5 <= d1 <= 15},
-  ///    output = {0 <= d0 <= 6, 1 <= d1 <= 15}
-  LogicalResult unionBoundingBox(const FlatAffineConstraints &other);
-
   /// Replaces the contents of this FlatAffineConstraints with `other`.
   void clearAndCopyFrom(const IntegerPolyhedron &other) override;
 
-  /// Returns the smallest known constant bound for the extent of the specified
-  /// identifier (pos^th), i.e., the smallest known constant that is greater
-  /// than or equal to 'exclusive upper bound' - 'lower bound' of the
-  /// identifier. This constant bound is guaranteed to be non-negative. Returns
-  /// None if it's not a constant. This method employs trivial (low complexity /
-  /// cost) checks and detection. Symbolic identifiers are treated specially,
-  /// i.e., it looks for constant 
diff erences between affine expressions
-  /// involving only the symbolic identifiers. `lb` and `ub` (along with the
-  /// `boundFloorDivisor`) are set to represent the lower and upper bound
-  /// associated with the constant 
diff erence: `lb`, `ub` have the coefficients,
-  /// and `boundFloorDivisor`, their divisor. `minLbPos` and `minUbPos` if
-  /// non-null are set to the position of the constant lower bound and upper
-  /// bound respectively (to the same if they are from an equality). Ex: if the
-  /// lower bound is [(s0 + s2 - 1) floordiv 32] for a system with three
-  /// symbolic identifiers, *lb = [1, 0, 1], lbDivisor = 32. See comments at
-  /// function definition for examples.
-  Optional<int64_t> getConstantBoundOnDimSize(
-      unsigned pos, SmallVectorImpl<int64_t> *lb = nullptr,
-      int64_t *boundFloorDivisor = nullptr,
-      SmallVectorImpl<int64_t> *ub = nullptr, unsigned *minLbPos = nullptr,
-      unsigned *minUbPos = nullptr) const;
-
-  /// Returns the constant bound for the pos^th identifier if there is one;
-  /// None otherwise.
-  // TODO: Support EQ bounds.
-  Optional<int64_t> getConstantBound(BoundType type, unsigned pos) const;
-
   /// Gets the lower and upper bound of the `offset` + `pos`th identifier
   /// treating [0, offset) U [offset + num, symStartPos) as dimensions and
   /// [symStartPos, getNumDimAndSymbolIds) as symbols, and `pos` lies in
@@ -290,61 +155,7 @@ class FlatAffineConstraints : public IntegerPolyhedron {
                         unsigned symStartPos, ArrayRef<AffineExpr> localExprs,
                         MLIRContext *context) const;
 
-  /// Removes constraints that are independent of (i.e., do not have a
-  /// coefficient) identifiers in the range [pos, pos + num).
-  void removeIndependentConstraints(unsigned pos, unsigned num);
-
-  /// Returns true if the set can be trivially detected as being
-  /// hyper-rectangular on the specified contiguous set of identifiers.
-  bool isHyperRectangular(unsigned pos, unsigned num) const;
-
-  /// Removes duplicate constraints, trivially true constraints, and constraints
-  /// that can be detected as redundant as a result of 
diff ering only in their
-  /// constant term part. A constraint of the form <non-negative constant> >= 0
-  /// is considered trivially true. This method is a linear time method on the
-  /// constraints, does a single scan, and updates in place. It also normalizes
-  /// constraints by their GCD and performs GCD tightening on inequalities.
-  void removeTrivialRedundancy();
-
-  /// A more expensive check than `removeTrivialRedundancy` to detect redundant
-  /// inequalities.
-  void removeRedundantInequalities();
-
-  /// Removes redundant constraints using Simplex. Although the algorithm can
-  /// theoretically take exponential time in the worst case (rare), it is known
-  /// to perform much better in the average case. If V is the number of vertices
-  /// in the polytope and C is the number of constraints, the algorithm takes
-  /// O(VC) time.
-  void removeRedundantConstraints();
-
-  /// Converts identifiers in the column range [idStart, idLimit) to local
-  /// variables.
-  void convertDimToLocal(unsigned dimStart, unsigned dimLimit);
-
-  /// Adds additional local ids to the sets such that they both have the union
-  /// of the local ids in each set, without changing the set of points that
-  /// lie in `this` and `other`. The ordering of the local ids in the
-  /// sets may also be changed. After merging, if the `i^th` local variable in
-  /// one set has a known division representation, then the `i^th` local
-  /// variable in the other set either has the same division representation or
-  /// no known division representation.
-  ///
-  /// The number of dimensions and symbol ids in `this` and `other` should
-  /// match.
-  void mergeLocalIds(FlatAffineConstraints &other);
-
 protected:
-  /// Checks all rows of equality/inequality constraints for trivial
-  /// contradictions (for example: 1 == 0, 0 >= 1), which may have surfaced
-  /// after elimination. Returns true if an invalid constraint is found;
-  /// false otherwise.
-  bool hasInvalidConstraint() const;
-
-  /// Returns the constant lower bound bound if isLower is true, and the upper
-  /// bound if isLower is false.
-  template <bool isLower>
-  Optional<int64_t> computeConstantLowerOrUpperBound(unsigned pos);
-
   /// Given an affine map that is aligned with this constraint system:
   /// * Flatten the map.
   /// * Add newly introduced local columns at the beginning of this constraint
@@ -358,73 +169,10 @@ class FlatAffineConstraints : public IntegerPolyhedron {
   LogicalResult flattenAlignedMapAndMergeLocals(
       AffineMap map, std::vector<SmallVector<int64_t, 8>> *flattenedExprs);
 
-  /// Eliminates a single identifier at `position` from equality and inequality
-  /// constraints. Returns `success` if the identifier was eliminated, and
-  /// `failure` otherwise.
-  inline LogicalResult gaussianEliminateId(unsigned position) {
-    return success(gaussianEliminateIds(position, position + 1) == 1);
-  }
-
-  /// Removes local variables using equalities. Each equality is checked if it
-  /// can be reduced to the form: `e = affine-expr`, where `e` is a local
-  /// variable and `affine-expr` is an affine expression not containing `e`.
-  /// If an equality satisfies this form, the local variable is replaced in
-  /// each constraint and then removed. The equality used to replace this local
-  /// variable is also removed.
-  void removeRedundantLocalVars();
-
-  /// Eliminates identifiers from equality and inequality constraints
-  /// in column range [posStart, posLimit).
-  /// Returns the number of variables eliminated.
-  unsigned gaussianEliminateIds(unsigned posStart, unsigned posLimit);
-
-  /// Eliminates the identifier at the specified position using Fourier-Motzkin
-  /// variable elimination, but uses Gaussian elimination if there is an
-  /// equality involving that identifier. If the result of the elimination is
-  /// integer exact, `*isResultIntegerExact` is set to true. If `darkShadow` is
-  /// set to true, a potential under approximation (subset) of the rational
-  /// shadow / exact integer shadow is computed.
-  // See implementation comments for more details.
-  virtual void fourierMotzkinEliminate(unsigned pos, bool darkShadow = false,
-                                       bool *isResultIntegerExact = nullptr);
-
-  /// Tightens inequalities given that we are dealing with integer spaces. This
-  /// is similar to the GCD test but applied to inequalities. The constant term
-  /// can be reduced to the preceding multiple of the GCD of the coefficients,
-  /// i.e.,
-  ///  64*i - 100 >= 0  =>  64*i - 128 >= 0 (since 'i' is an integer). This is a
-  /// fast method (linear in the number of coefficients).
-  void gcdTightenInequalities();
-
-  /// Normalized each constraints by the GCD of its coefficients.
-  void normalizeConstraintsByGCD();
-
-  /// Searches for a constraint with a non-zero coefficient at `colIdx` in
-  /// equality (isEq=true) or inequality (isEq=false) constraints.
-  /// Returns true and sets row found in search in `rowIdx`, false otherwise.
-  bool findConstraintWithNonZeroAt(unsigned colIdx, bool isEq,
-                                   unsigned *rowIdx) const;
-
-  /// Returns true if the pos^th column is all zero for both inequalities and
-  /// equalities.
-  bool isColZero(unsigned pos) const;
-
   /// Prints the number of constraints, dimensions, symbols and locals in the
   /// FlatAffineConstraints. Also, prints for each identifier whether there is
   /// an SSA Value attached to it.
   void printSpace(raw_ostream &os) const override;
-
-  /// A parameter that controls detection of an unrealistic number of
-  /// constraints. If the number of constraints is this many times the number of
-  /// variables, we consider such a system out of line with the intended use
-  /// case of FlatAffineConstraints.
-  // The rationale for 32 is that in the typical simplest of cases, an
-  // identifier is expected to have one lower bound and one upper bound
-  // constraint. With a level of tiling or a connection to another identifier
-  // through a div or mod, an extra pair of bounds gets added. As a limit, we
-  // don't expect an identifier to have more than 32 lower/upper/equality
-  // constraints. This is conservatively set low and can be raised if needed.
-  constexpr static unsigned kExplosionFactor = 32;
 };
 
 /// An extension of FlatAffineConstraints in which dimensions and symbols can

diff  --git a/mlir/include/mlir/Analysis/Presburger/IntegerPolyhedron.h b/mlir/include/mlir/Analysis/Presburger/IntegerPolyhedron.h
index eae56e3acd0a..21e31cee1abe 100644
--- a/mlir/include/mlir/Analysis/Presburger/IntegerPolyhedron.h
+++ b/mlir/include/mlir/Analysis/Presburger/IntegerPolyhedron.h
@@ -14,6 +14,7 @@
 #define MLIR_ANALYSIS_PRESBURGER_INTEGERPOLYHEDRON_H
 
 #include "mlir/Analysis/Presburger/Matrix.h"
+#include "mlir/Support/LogicalResult.h"
 
 namespace mlir {
 
@@ -214,10 +215,254 @@ class IntegerPolyhedron {
                                SmallVectorImpl<unsigned> *eqIndices = nullptr,
                                unsigned offset = 0, unsigned num = 0) const;
 
+  /// Checks for emptiness by performing variable elimination on all
+  /// identifiers, running the GCD test on each equality constraint, and
+  /// checking for invalid constraints. Returns true if the GCD test fails for
+  /// any equality, or if any invalid constraints are discovered on any row.
+  /// Returns false otherwise.
+  bool isEmpty() const;
+
+  /// Runs the GCD test on all equality constraints. Returns true if this test
+  /// fails on any equality. Returns false otherwise.
+  /// This test can be used to disprove the existence of a solution. If it
+  /// returns true, no integer solution to the equality constraints can exist.
+  bool isEmptyByGCDTest() const;
+
+  /// Returns true if the set of constraints is found to have no solution,
+  /// false if a solution exists. Uses the same algorithm as
+  /// `findIntegerSample`.
+  bool isIntegerEmpty() const;
+
+  /// Returns a matrix where each row is a vector along which the polytope is
+  /// bounded. The span of the returned vectors is guaranteed to contain all
+  /// such vectors. The returned vectors are NOT guaranteed to be linearly
+  /// independent. This function should not be called on empty sets.
+  Matrix getBoundedDirections() const;
+
+  /// Find an integer sample point satisfying the constraints using a
+  /// branch and bound algorithm with generalized basis reduction, with some
+  /// additional processing using Simplex for unbounded sets.
+  ///
+  /// Returns an integer sample point if one exists, or an empty Optional
+  /// otherwise.
+  Optional<SmallVector<int64_t, 8>> findIntegerSample() const;
+
+  /// Returns true if the given point satisfies the constraints, or false
+  /// otherwise.
+  bool containsPoint(ArrayRef<int64_t> point) const;
+
+  /// Find pairs of inequalities identified by their position indices, using
+  /// which an explicit representation for each local variable can be computed.
+  /// The pairs are stored as indices of upperbound, lowerbound inequalities. If
+  /// no such pair can be found, it is stored as llvm::None.
+  ///
+  /// The dividends of the explicit representations are stored in `dividends`
+  /// and the denominators in `denominators`. If no explicit representation
+  /// could be found for the `i^th` local identifier, `denominators[i]` is set
+  /// to 0.
+  void getLocalReprs(
+      std::vector<SmallVector<int64_t, 8>> &dividends,
+      SmallVector<unsigned, 4> &denominators,
+      std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const;
+  void getLocalReprs(
+      std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const;
+  void getLocalReprs(std::vector<SmallVector<int64_t, 8>> &dividends,
+                     SmallVector<unsigned, 4> &denominators) const;
+
+  /// The type of bound: equal, lower bound or upper bound.
+  enum BoundType { EQ, LB, UB };
+
+  /// Adds a constant bound for the specified identifier.
+  void addBound(BoundType type, unsigned pos, int64_t value);
+
+  /// Adds a constant bound for the specified expression.
+  void addBound(BoundType type, ArrayRef<int64_t> expr, int64_t value);
+
+  /// Adds a new local identifier as the floordiv of an affine function of other
+  /// identifiers, the coefficients of which are provided in `dividend` and with
+  /// respect to a positive constant `divisor`. Two constraints are added to the
+  /// system to capture equivalence with the floordiv:
+  /// q = dividend floordiv c    <=>   c*q <= dividend <= c*q + c - 1.
+  void addLocalFloorDiv(ArrayRef<int64_t> dividend, int64_t divisor);
+
+  /// Projects out (aka eliminates) `num` identifiers starting at position
+  /// `pos`. The resulting constraint system is the shadow along the dimensions
+  /// that still exist. This method may not always be integer exact.
+  // TODO: deal with integer exactness when necessary - can return a value to
+  // mark exactness for example.
+  void projectOut(unsigned pos, unsigned num);
+  inline void projectOut(unsigned pos) { return projectOut(pos, 1); }
+
+  /// Changes the partition between dimensions and symbols. Depending on the new
+  /// symbol count, either a chunk of trailing dimensional identifiers becomes
+  /// symbols, or some of the leading symbols become dimensions.
+  void setDimSymbolSeparation(unsigned newSymbolCount);
+
+  /// Tries to fold the specified identifier to a constant using a trivial
+  /// equality detection; if successful, the constant is substituted for the
+  /// identifier everywhere in the constraint system and then removed from the
+  /// system.
+  LogicalResult constantFoldId(unsigned pos);
+
+  /// This method calls `constantFoldId` for the specified range of identifiers,
+  /// `num` identifiers starting at position `pos`.
+  void constantFoldIdRange(unsigned pos, unsigned num);
+
+  /// Updates the constraints to be the smallest bounding (enclosing) box that
+  /// contains the points of `this` set and that of `other`, with the symbols
+  /// being treated specially. For each of the dimensions, the min of the lower
+  /// bounds (symbolic) and the max of the upper bounds (symbolic) is computed
+  /// to determine such a bounding box. `other` is expected to have the same
+  /// dimensional identifiers as this constraint system (in the same order).
+  ///
+  /// E.g.:
+  /// 1) this   = {0 <= d0 <= 127},
+  ///    other  = {16 <= d0 <= 192},
+  ///    output = {0 <= d0 <= 192}
+  /// 2) this   = {s0 + 5 <= d0 <= s0 + 20},
+  ///    other  = {s0 + 1 <= d0 <= s0 + 9},
+  ///    output = {s0 + 1 <= d0 <= s0 + 20}
+  /// 3) this   = {0 <= d0 <= 5, 1 <= d1 <= 9}
+  ///    other  = {2 <= d0 <= 6, 5 <= d1 <= 15},
+  ///    output = {0 <= d0 <= 6, 1 <= d1 <= 15}
+  LogicalResult unionBoundingBox(const IntegerPolyhedron &other);
+
+  /// Returns the smallest known constant bound for the extent of the specified
+  /// identifier (pos^th), i.e., the smallest known constant that is greater
+  /// than or equal to 'exclusive upper bound' - 'lower bound' of the
+  /// identifier. This constant bound is guaranteed to be non-negative. Returns
+  /// None if it's not a constant. This method employs trivial (low complexity /
+  /// cost) checks and detection. Symbolic identifiers are treated specially,
+  /// i.e., it looks for constant 
diff erences between affine expressions
+  /// involving only the symbolic identifiers. `lb` and `ub` (along with the
+  /// `boundFloorDivisor`) are set to represent the lower and upper bound
+  /// associated with the constant 
diff erence: `lb`, `ub` have the coefficients,
+  /// and `boundFloorDivisor`, their divisor. `minLbPos` and `minUbPos` if
+  /// non-null are set to the position of the constant lower bound and upper
+  /// bound respectively (to the same if they are from an equality). Ex: if the
+  /// lower bound is [(s0 + s2 - 1) floordiv 32] for a system with three
+  /// symbolic identifiers, *lb = [1, 0, 1], lbDivisor = 32. See comments at
+  /// function definition for examples.
+  Optional<int64_t> getConstantBoundOnDimSize(
+      unsigned pos, SmallVectorImpl<int64_t> *lb = nullptr,
+      int64_t *boundFloorDivisor = nullptr,
+      SmallVectorImpl<int64_t> *ub = nullptr, unsigned *minLbPos = nullptr,
+      unsigned *minUbPos = nullptr) const;
+
+  /// Returns the constant bound for the pos^th identifier if there is one;
+  /// None otherwise.
+  // TODO: Support EQ bounds.
+  Optional<int64_t> getConstantBound(BoundType type, unsigned pos) const;
+
+  /// Removes constraints that are independent of (i.e., do not have a
+  /// coefficient) identifiers in the range [pos, pos + num).
+  void removeIndependentConstraints(unsigned pos, unsigned num);
+
+  /// Returns true if the set can be trivially detected as being
+  /// hyper-rectangular on the specified contiguous set of identifiers.
+  bool isHyperRectangular(unsigned pos, unsigned num) const;
+
+  /// Removes duplicate constraints, trivially true constraints, and constraints
+  /// that can be detected as redundant as a result of 
diff ering only in their
+  /// constant term part. A constraint of the form <non-negative constant> >= 0
+  /// is considered trivially true. This method is a linear time method on the
+  /// constraints, does a single scan, and updates in place. It also normalizes
+  /// constraints by their GCD and performs GCD tightening on inequalities.
+  void removeTrivialRedundancy();
+
+  /// A more expensive check than `removeTrivialRedundancy` to detect redundant
+  /// inequalities.
+  void removeRedundantInequalities();
+
+  /// Removes redundant constraints using Simplex. Although the algorithm can
+  /// theoretically take exponential time in the worst case (rare), it is known
+  /// to perform much better in the average case. If V is the number of vertices
+  /// in the polytope and C is the number of constraints, the algorithm takes
+  /// O(VC) time.
+  void removeRedundantConstraints();
+
+  /// Converts identifiers in the column range [idStart, idLimit) to local
+  /// variables.
+  void convertDimToLocal(unsigned dimStart, unsigned dimLimit);
+
+  /// Adds additional local ids to the sets such that they both have the union
+  /// of the local ids in each set, without changing the set of points that
+  /// lie in `this` and `other`. The ordering of the local ids in the
+  /// sets may also be changed. After merging, if the `i^th` local variable in
+  /// one set has a known division representation, then the `i^th` local
+  /// variable in the other set either has the same division representation or
+  /// no known division representation.
+  ///
+  /// The number of dimensions and symbol ids in `this` and `other` should
+  /// match.
+  void mergeLocalIds(IntegerPolyhedron &other);
+
   void print(raw_ostream &os) const;
   void dump() const;
 
 protected:
+  /// Checks all rows of equality/inequality constraints for trivial
+  /// contradictions (for example: 1 == 0, 0 >= 1), which may have surfaced
+  /// after elimination. Returns true if an invalid constraint is found;
+  /// false otherwise.
+  bool hasInvalidConstraint() const;
+
+  /// Returns the constant lower bound bound if isLower is true, and the upper
+  /// bound if isLower is false.
+  template <bool isLower>
+  Optional<int64_t> computeConstantLowerOrUpperBound(unsigned pos);
+
+  /// Eliminates a single identifier at `position` from equality and inequality
+  /// constraints. Returns `success` if the identifier was eliminated, and
+  /// `failure` otherwise.
+  inline LogicalResult gaussianEliminateId(unsigned position) {
+    return success(gaussianEliminateIds(position, position + 1) == 1);
+  }
+
+  /// Removes local variables using equalities. Each equality is checked if it
+  /// can be reduced to the form: `e = affine-expr`, where `e` is a local
+  /// variable and `affine-expr` is an affine expression not containing `e`.
+  /// If an equality satisfies this form, the local variable is replaced in
+  /// each constraint and then removed. The equality used to replace this local
+  /// variable is also removed.
+  void removeRedundantLocalVars();
+
+  /// Eliminates identifiers from equality and inequality constraints
+  /// in column range [posStart, posLimit).
+  /// Returns the number of variables eliminated.
+  unsigned gaussianEliminateIds(unsigned posStart, unsigned posLimit);
+
+  /// Eliminates the identifier at the specified position using Fourier-Motzkin
+  /// variable elimination, but uses Gaussian elimination if there is an
+  /// equality involving that identifier. If the result of the elimination is
+  /// integer exact, `*isResultIntegerExact` is set to true. If `darkShadow` is
+  /// set to true, a potential under approximation (subset) of the rational
+  /// shadow / exact integer shadow is computed.
+  // See implementation comments for more details.
+  virtual void fourierMotzkinEliminate(unsigned pos, bool darkShadow = false,
+                                       bool *isResultIntegerExact = nullptr);
+
+  /// Tightens inequalities given that we are dealing with integer spaces. This
+  /// is similar to the GCD test but applied to inequalities. The constant term
+  /// can be reduced to the preceding multiple of the GCD of the coefficients,
+  /// i.e.,
+  ///  64*i - 100 >= 0  =>  64*i - 128 >= 0 (since 'i' is an integer). This is a
+  /// fast method (linear in the number of coefficients).
+  void gcdTightenInequalities();
+
+  /// Normalized each constraints by the GCD of its coefficients.
+  void normalizeConstraintsByGCD();
+
+  /// Searches for a constraint with a non-zero coefficient at `colIdx` in
+  /// equality (isEq=true) or inequality (isEq=false) constraints.
+  /// Returns true and sets row found in search in `rowIdx`, false otherwise.
+  bool findConstraintWithNonZeroAt(unsigned colIdx, bool isEq,
+                                   unsigned *rowIdx) const;
+
+  /// Returns true if the pos^th column is all zero for both inequalities and
+  /// equalities.
+  bool isColZero(unsigned pos) const;
+
   /// Returns false if the fields corresponding to various identifier counts, or
   /// equality/inequality buffer sizes aren't consistent; true otherwise. This
   /// is meant to be used within an assert internally.
@@ -238,6 +483,18 @@ class IntegerPolyhedron {
   /// arrays as needed.
   virtual void removeIdRange(unsigned idStart, unsigned idLimit);
 
+  /// A parameter that controls detection of an unrealistic number of
+  /// constraints. If the number of constraints is this many times the number of
+  /// variables, we consider such a system out of line with the intended use
+  /// case of IntegerPolyhedron.
+  // The rationale for 32 is that in the typical simplest of cases, an
+  // identifier is expected to have one lower bound and one upper bound
+  // constraint. With a level of tiling or a connection to another identifier
+  // through a div or mod, an extra pair of bounds gets added. As a limit, we
+  // don't expect an identifier to have more than 32 lower/upper/equality
+  // constraints. This is conservatively set low and can be raised if needed.
+  constexpr static unsigned kExplosionFactor = 32;
+
   /// Total number of identifiers.
   unsigned numIds;
 

diff  --git a/mlir/lib/Analysis/AffineStructures.cpp b/mlir/lib/Analysis/AffineStructures.cpp
index 6c11bf74a414..0e25c640e1b0 100644
--- a/mlir/lib/Analysis/AffineStructures.cpp
+++ b/mlir/lib/Analysis/AffineStructures.cpp
@@ -31,8 +31,6 @@
 #define DEBUG_TYPE "affine-structures"
 
 using namespace mlir;
-using llvm::SmallDenseMap;
-using llvm::SmallDenseSet;
 
 namespace {
 
@@ -698,627 +696,17 @@ void FlatAffineValueConstraints::addAffineIfOpDomain(AffineIfOp ifOp) {
   append(cst);
 }
 
-// Searches for a constraint with a non-zero coefficient at `colIdx` in
-// equality (isEq=true) or inequality (isEq=false) constraints.
-// Returns true and sets row found in search in `rowIdx`, false otherwise.
-bool FlatAffineConstraints::findConstraintWithNonZeroAt(
-    unsigned colIdx, bool isEq, unsigned *rowIdx) const {
-  assert(colIdx < getNumCols() && "position out of bounds");
-  auto at = [&](unsigned rowIdx) -> int64_t {
-    return isEq ? atEq(rowIdx, colIdx) : atIneq(rowIdx, colIdx);
-  };
-  unsigned e = isEq ? getNumEqualities() : getNumInequalities();
-  for (*rowIdx = 0; *rowIdx < e; ++(*rowIdx)) {
-    if (at(*rowIdx) != 0) {
-      return true;
-    }
-  }
-  return false;
-}
-
-// Normalizes the coefficient values across all columns in `rowIdx` by their
-// GCD in equality or inequality constraints as specified by `isEq`.
-template <bool isEq>
-static void normalizeConstraintByGCD(FlatAffineConstraints *constraints,
-                                     unsigned rowIdx) {
-  auto at = [&](unsigned colIdx) -> int64_t {
-    return isEq ? constraints->atEq(rowIdx, colIdx)
-                : constraints->atIneq(rowIdx, colIdx);
-  };
-  uint64_t gcd = std::abs(at(0));
-  for (unsigned j = 1, e = constraints->getNumCols(); j < e; ++j) {
-    gcd = llvm::GreatestCommonDivisor64(gcd, std::abs(at(j)));
-  }
-  if (gcd > 0 && gcd != 1) {
-    for (unsigned j = 0, e = constraints->getNumCols(); j < e; ++j) {
-      int64_t v = at(j) / static_cast<int64_t>(gcd);
-      isEq ? constraints->atEq(rowIdx, j) = v
-           : constraints->atIneq(rowIdx, j) = v;
-    }
-  }
-}
-
-void FlatAffineConstraints::normalizeConstraintsByGCD() {
-  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
-    normalizeConstraintByGCD</*isEq=*/true>(this, i);
-  }
-  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
-    normalizeConstraintByGCD</*isEq=*/false>(this, i);
-  }
-}
-
 bool FlatAffineValueConstraints::hasConsistentState() const {
   return FlatAffineConstraints::hasConsistentState() &&
          values.size() == getNumIds();
 }
 
-bool FlatAffineConstraints::hasInvalidConstraint() const {
-  assert(hasConsistentState());
-  auto check = [&](bool isEq) -> bool {
-    unsigned numCols = getNumCols();
-    unsigned numRows = isEq ? getNumEqualities() : getNumInequalities();
-    for (unsigned i = 0, e = numRows; i < e; ++i) {
-      unsigned j;
-      for (j = 0; j < numCols - 1; ++j) {
-        int64_t v = isEq ? atEq(i, j) : atIneq(i, j);
-        // Skip rows with non-zero variable coefficients.
-        if (v != 0)
-          break;
-      }
-      if (j < numCols - 1) {
-        continue;
-      }
-      // Check validity of constant term at 'numCols - 1' w.r.t 'isEq'.
-      // Example invalid constraints include: '1 == 0' or '-1 >= 0'
-      int64_t v = isEq ? atEq(i, numCols - 1) : atIneq(i, numCols - 1);
-      if ((isEq && v != 0) || (!isEq && v < 0)) {
-        return true;
-      }
-    }
-    return false;
-  };
-  if (check(/*isEq=*/true))
-    return true;
-  return check(/*isEq=*/false);
-}
-
-/// Eliminate identifier from constraint at `rowIdx` based on coefficient at
-/// pivotRow, pivotCol. Columns in range [elimColStart, pivotCol) will not be
-/// updated as they have already been eliminated.
-static void eliminateFromConstraint(FlatAffineConstraints *constraints,
-                                    unsigned rowIdx, unsigned pivotRow,
-                                    unsigned pivotCol, unsigned elimColStart,
-                                    bool isEq) {
-  // Skip if equality 'rowIdx' if same as 'pivotRow'.
-  if (isEq && rowIdx == pivotRow)
-    return;
-  auto at = [&](unsigned i, unsigned j) -> int64_t {
-    return isEq ? constraints->atEq(i, j) : constraints->atIneq(i, j);
-  };
-  int64_t leadCoeff = at(rowIdx, pivotCol);
-  // Skip if leading coefficient at 'rowIdx' is already zero.
-  if (leadCoeff == 0)
-    return;
-  int64_t pivotCoeff = constraints->atEq(pivotRow, pivotCol);
-  int64_t sign = (leadCoeff * pivotCoeff > 0) ? -1 : 1;
-  int64_t lcm = mlir::lcm(pivotCoeff, leadCoeff);
-  int64_t pivotMultiplier = sign * (lcm / std::abs(pivotCoeff));
-  int64_t rowMultiplier = lcm / std::abs(leadCoeff);
-
-  unsigned numCols = constraints->getNumCols();
-  for (unsigned j = 0; j < numCols; ++j) {
-    // Skip updating column 'j' if it was just eliminated.
-    if (j >= elimColStart && j < pivotCol)
-      continue;
-    int64_t v = pivotMultiplier * constraints->atEq(pivotRow, j) +
-                rowMultiplier * at(rowIdx, j);
-    isEq ? constraints->atEq(rowIdx, j) = v
-         : constraints->atIneq(rowIdx, j) = v;
-  }
-}
-
 void FlatAffineValueConstraints::removeIdRange(unsigned idStart,
                                                unsigned idLimit) {
   FlatAffineConstraints::removeIdRange(idStart, idLimit);
   values.erase(values.begin() + idStart, values.begin() + idLimit);
 }
 
-/// Returns the position of the identifier that has the minimum <number of lower
-/// bounds> times <number of upper bounds> from the specified range of
-/// identifiers [start, end). It is often best to eliminate in the increasing
-/// order of these counts when doing Fourier-Motzkin elimination since FM adds
-/// that many new constraints.
-static unsigned getBestIdToEliminate(const FlatAffineConstraints &cst,
-                                     unsigned start, unsigned end) {
-  assert(start < cst.getNumIds() && end < cst.getNumIds() + 1);
-
-  auto getProductOfNumLowerUpperBounds = [&](unsigned pos) {
-    unsigned numLb = 0;
-    unsigned numUb = 0;
-    for (unsigned r = 0, e = cst.getNumInequalities(); r < e; r++) {
-      if (cst.atIneq(r, pos) > 0) {
-        ++numLb;
-      } else if (cst.atIneq(r, pos) < 0) {
-        ++numUb;
-      }
-    }
-    return numLb * numUb;
-  };
-
-  unsigned minLoc = start;
-  unsigned min = getProductOfNumLowerUpperBounds(start);
-  for (unsigned c = start + 1; c < end; c++) {
-    unsigned numLbUbProduct = getProductOfNumLowerUpperBounds(c);
-    if (numLbUbProduct < min) {
-      min = numLbUbProduct;
-      minLoc = c;
-    }
-  }
-  return minLoc;
-}
-
-// Checks for emptiness of the set by eliminating identifiers successively and
-// using the GCD test (on all equality constraints) and checking for trivially
-// invalid constraints. Returns 'true' if the constraint system is found to be
-// empty; false otherwise.
-bool FlatAffineConstraints::isEmpty() const {
-  if (isEmptyByGCDTest() || hasInvalidConstraint())
-    return true;
-
-  FlatAffineConstraints tmpCst(*this);
-
-  // First, eliminate as many local variables as possible using equalities.
-  tmpCst.removeRedundantLocalVars();
-  if (tmpCst.isEmptyByGCDTest() || tmpCst.hasInvalidConstraint())
-    return true;
-
-  // Eliminate as many identifiers as possible using Gaussian elimination.
-  unsigned currentPos = 0;
-  while (currentPos < tmpCst.getNumIds()) {
-    tmpCst.gaussianEliminateIds(currentPos, tmpCst.getNumIds());
-    ++currentPos;
-    // We check emptiness through trivial checks after eliminating each ID to
-    // detect emptiness early. Since the checks isEmptyByGCDTest() and
-    // hasInvalidConstraint() are linear time and single sweep on the constraint
-    // buffer, this appears reasonable - but can optimize in the future.
-    if (tmpCst.hasInvalidConstraint() || tmpCst.isEmptyByGCDTest())
-      return true;
-  }
-
-  // Eliminate the remaining using FM.
-  for (unsigned i = 0, e = tmpCst.getNumIds(); i < e; i++) {
-    tmpCst.fourierMotzkinEliminate(
-        getBestIdToEliminate(tmpCst, 0, tmpCst.getNumIds()));
-    // Check for a constraint explosion. This rarely happens in practice, but
-    // this check exists as a safeguard against improperly constructed
-    // constraint systems or artificially created arbitrarily complex systems
-    // that aren't the intended use case for FlatAffineConstraints. This is
-    // needed since FM has a worst case exponential complexity in theory.
-    if (tmpCst.getNumConstraints() >= kExplosionFactor * getNumIds()) {
-      LLVM_DEBUG(llvm::dbgs() << "FM constraint explosion detected\n");
-      return false;
-    }
-
-    // FM wouldn't have modified the equalities in any way. So no need to again
-    // run GCD test. Check for trivial invalid constraints.
-    if (tmpCst.hasInvalidConstraint())
-      return true;
-  }
-  return false;
-}
-
-// Runs the GCD test on all equality constraints. Returns 'true' if this test
-// fails on any equality. Returns 'false' otherwise.
-// This test can be used to disprove the existence of a solution. If it returns
-// true, no integer solution to the equality constraints can exist.
-//
-// GCD test definition:
-//
-// The equality constraint:
-//
-//  c_1*x_1 + c_2*x_2 + ... + c_n*x_n = c_0
-//
-// has an integer solution iff:
-//
-//  GCD of c_1, c_2, ..., c_n divides c_0.
-//
-bool FlatAffineConstraints::isEmptyByGCDTest() const {
-  assert(hasConsistentState());
-  unsigned numCols = getNumCols();
-  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
-    uint64_t gcd = std::abs(atEq(i, 0));
-    for (unsigned j = 1; j < numCols - 1; ++j) {
-      gcd = llvm::GreatestCommonDivisor64(gcd, std::abs(atEq(i, j)));
-    }
-    int64_t v = std::abs(atEq(i, numCols - 1));
-    if (gcd > 0 && (v % gcd != 0)) {
-      return true;
-    }
-  }
-  return false;
-}
-
-// Returns a matrix where each row is a vector along which the polytope is
-// bounded. The span of the returned vectors is guaranteed to contain all
-// such vectors. The returned vectors are NOT guaranteed to be linearly
-// independent. This function should not be called on empty sets.
-//
-// It is sufficient to check the perpendiculars of the constraints, as the set
-// of perpendiculars which are bounded must span all bounded directions.
-Matrix FlatAffineConstraints::getBoundedDirections() const {
-  // Note that it is necessary to add the equalities too (which the constructor
-  // does) even though we don't need to check if they are bounded; whether an
-  // inequality is bounded or not depends on what other constraints, including
-  // equalities, are present.
-  Simplex simplex(*this);
-
-  assert(!simplex.isEmpty() && "It is not meaningful to ask whether a "
-                               "direction is bounded in an empty set.");
-
-  SmallVector<unsigned, 8> boundedIneqs;
-  // The constructor adds the inequalities to the simplex first, so this
-  // processes all the inequalities.
-  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
-    if (simplex.isBoundedAlongConstraint(i))
-      boundedIneqs.push_back(i);
-  }
-
-  // The direction vector is given by the coefficients and does not include the
-  // constant term, so the matrix has one fewer column.
-  unsigned dirsNumCols = getNumCols() - 1;
-  Matrix dirs(boundedIneqs.size() + getNumEqualities(), dirsNumCols);
-
-  // Copy the bounded inequalities.
-  unsigned row = 0;
-  for (unsigned i : boundedIneqs) {
-    for (unsigned col = 0; col < dirsNumCols; ++col)
-      dirs(row, col) = atIneq(i, col);
-    ++row;
-  }
-
-  // Copy the equalities. All the equalities' perpendiculars are bounded.
-  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
-    for (unsigned col = 0; col < dirsNumCols; ++col)
-      dirs(row, col) = atEq(i, col);
-    ++row;
-  }
-
-  return dirs;
-}
-
-bool eqInvolvesSuffixDims(const FlatAffineConstraints &fac, unsigned eqIndex,
-                          unsigned numDims) {
-  for (unsigned e = fac.getNumIds(), j = e - numDims; j < e; ++j)
-    if (fac.atEq(eqIndex, j) != 0)
-      return true;
-  return false;
-}
-bool ineqInvolvesSuffixDims(const FlatAffineConstraints &fac,
-                            unsigned ineqIndex, unsigned numDims) {
-  for (unsigned e = fac.getNumIds(), j = e - numDims; j < e; ++j)
-    if (fac.atIneq(ineqIndex, j) != 0)
-      return true;
-  return false;
-}
-
-void removeConstraintsInvolvingSuffixDims(FlatAffineConstraints &fac,
-                                          unsigned unboundedDims) {
-  // We iterate backwards so that whether we remove constraint i - 1 or not, the
-  // next constraint to be tested is always i - 2.
-  for (unsigned i = fac.getNumEqualities(); i > 0; i--)
-    if (eqInvolvesSuffixDims(fac, i - 1, unboundedDims))
-      fac.removeEquality(i - 1);
-  for (unsigned i = fac.getNumInequalities(); i > 0; i--)
-    if (ineqInvolvesSuffixDims(fac, i - 1, unboundedDims))
-      fac.removeInequality(i - 1);
-}
-
-bool FlatAffineConstraints::isIntegerEmpty() const {
-  return !findIntegerSample().hasValue();
-}
-
-/// Let this set be S. If S is bounded then we directly call into the GBR
-/// sampling algorithm. Otherwise, there are some unbounded directions, i.e.,
-/// vectors v such that S extends to infinity along v or -v. In this case we
-/// use an algorithm described in the integer set library (isl) manual and used
-/// by the isl_set_sample function in that library. The algorithm is:
-///
-/// 1) Apply a unimodular transform T to S to obtain S*T, such that all
-/// dimensions in which S*T is bounded lie in the linear span of a prefix of the
-/// dimensions.
-///
-/// 2) Construct a set B by removing all constraints that involve
-/// the unbounded dimensions and then deleting the unbounded dimensions. Note
-/// that B is a Bounded set.
-///
-/// 3) Try to obtain a sample from B using the GBR sampling
-/// algorithm. If no sample is found, return that S is empty.
-///
-/// 4) Otherwise, substitute the obtained sample into S*T to obtain a set
-/// C. C is a full-dimensional Cone and always contains a sample.
-///
-/// 5) Obtain an integer sample from C.
-///
-/// 6) Return T*v, where v is the concatenation of the samples from B and C.
-///
-/// The following is a sketch of a proof that
-/// a) If the algorithm returns empty, then S is empty.
-/// b) If the algorithm returns a sample, it is a valid sample in S.
-///
-/// The algorithm returns empty only if B is empty, in which case S*T is
-/// certainly empty since B was obtained by removing constraints and then
-/// deleting unconstrained dimensions from S*T. Since T is unimodular, a vector
-/// v is in S*T iff T*v is in S. So in this case, since
-/// S*T is empty, S is empty too.
-///
-/// Otherwise, the algorithm substitutes the sample from B into S*T. All the
-/// constraints of S*T that did not involve unbounded dimensions are satisfied
-/// by this substitution. All dimensions in the linear span of the dimensions
-/// outside the prefix are unbounded in S*T (step 1). Substituting values for
-/// the bounded dimensions cannot make these dimensions bounded, and these are
-/// the only remaining dimensions in C, so C is unbounded along every vector (in
-/// the positive or negative direction, or both). C is hence a full-dimensional
-/// cone and therefore always contains an integer point.
-///
-/// Concatenating the samples from B and C gives a sample v in S*T, so the
-/// returned sample T*v is a sample in S.
-Optional<SmallVector<int64_t, 8>>
-FlatAffineConstraints::findIntegerSample() const {
-  // First, try the GCD test heuristic.
-  if (isEmptyByGCDTest())
-    return {};
-
-  Simplex simplex(*this);
-  if (simplex.isEmpty())
-    return {};
-
-  // For a bounded set, we directly call into the GBR sampling algorithm.
-  if (!simplex.isUnbounded())
-    return simplex.findIntegerSample();
-
-  // The set is unbounded. We cannot directly use the GBR algorithm.
-  //
-  // m is a matrix containing, in each row, a vector in which S is
-  // bounded, such that the linear span of all these dimensions contains all
-  // bounded dimensions in S.
-  Matrix m = getBoundedDirections();
-  // In column echelon form, each row of m occupies only the first rank(m)
-  // columns and has zeros on the other columns. The transform T that brings S
-  // to column echelon form is unimodular as well, so this is a suitable
-  // transform to use in step 1 of the algorithm.
-  std::pair<unsigned, LinearTransform> result =
-      LinearTransform::makeTransformToColumnEchelon(std::move(m));
-  const LinearTransform &transform = result.second;
-  // 1) Apply T to S to obtain S*T.
-  IntegerPolyhedron transformedSet = transform.applyTo(*this);
-
-  // 2) Remove the unbounded dimensions and constraints involving them to
-  // obtain a bounded set.
-  FlatAffineConstraints boundedSet(transformedSet);
-  unsigned numBoundedDims = result.first;
-  unsigned numUnboundedDims = getNumIds() - numBoundedDims;
-  removeConstraintsInvolvingSuffixDims(boundedSet, numUnboundedDims);
-  boundedSet.removeIdRange(numBoundedDims, boundedSet.getNumIds());
-
-  // 3) Try to obtain a sample from the bounded set.
-  Optional<SmallVector<int64_t, 8>> boundedSample =
-      Simplex(boundedSet).findIntegerSample();
-  if (!boundedSample)
-    return {};
-  assert(boundedSet.containsPoint(*boundedSample) &&
-         "Simplex returned an invalid sample!");
-
-  // 4) Substitute the values of the bounded dimensions into S*T to obtain a
-  // full-dimensional cone, which necessarily contains an integer sample.
-  transformedSet.setAndEliminate(0, *boundedSample);
-  IntegerPolyhedron &cone = transformedSet;
-
-  // 5) Obtain an integer sample from the cone.
-  //
-  // We shrink the cone such that for any rational point in the shrunken cone,
-  // rounding up each of the point's coordinates produces a point that still
-  // lies in the original cone.
-  //
-  // Rounding up a point x adds a number e_i in [0, 1) to each coordinate x_i.
-  // For each inequality sum_i a_i x_i + c >= 0 in the original cone, the
-  // shrunken cone will have the inequality tightened by some amount s, such
-  // that if x satisfies the shrunken cone's tightened inequality, then x + e
-  // satisfies the original inequality, i.e.,
-  //
-  // sum_i a_i x_i + c + s >= 0 implies sum_i a_i (x_i + e_i) + c >= 0
-  //
-  // for any e_i values in [0, 1). In fact, we will handle the slightly more
-  // general case where e_i can be in [0, 1]. For example, consider the
-  // inequality 2x_1 - 3x_2 - 7x_3 - 6 >= 0, and let x = (3, 0, 0). How low
-  // could the LHS go if we added a number in [0, 1] to each coordinate? The LHS
-  // is minimized when we add 1 to the x_i with negative coefficient a_i and
-  // keep the other x_i the same. In the example, we would get x = (3, 1, 1),
-  // changing the value of the LHS by -3 + -7 = -10.
-  //
-  // In general, the value of the LHS can change by at most the sum of the
-  // negative a_i, so we accomodate this by shifting the inequality by this
-  // amount for the shrunken cone.
-  for (unsigned i = 0, e = cone.getNumInequalities(); i < e; ++i) {
-    for (unsigned j = 0; j < cone.getNumIds(); ++j) {
-      int64_t coeff = cone.atIneq(i, j);
-      if (coeff < 0)
-        cone.atIneq(i, cone.getNumIds()) += coeff;
-    }
-  }
-
-  // Obtain an integer sample in the cone by rounding up a rational point from
-  // the shrunken cone. Shrinking the cone amounts to shifting its apex
-  // "inwards" without changing its "shape"; the shrunken cone is still a
-  // full-dimensional cone and is hence non-empty.
-  Simplex shrunkenConeSimplex(cone);
-  assert(!shrunkenConeSimplex.isEmpty() && "Shrunken cone cannot be empty!");
-  SmallVector<Fraction, 8> shrunkenConeSample =
-      shrunkenConeSimplex.getRationalSample();
-
-  SmallVector<int64_t, 8> coneSample(llvm::map_range(shrunkenConeSample, ceil));
-
-  // 6) Return transform * concat(boundedSample, coneSample).
-  SmallVector<int64_t, 8> &sample = boundedSample.getValue();
-  sample.append(coneSample.begin(), coneSample.end());
-  return transform.preMultiplyColumn(sample);
-}
-
-/// Helper to evaluate an affine expression at a point.
-/// The expression is a list of coefficients for the dimensions followed by the
-/// constant term.
-static int64_t valueAt(ArrayRef<int64_t> expr, ArrayRef<int64_t> point) {
-  assert(expr.size() == 1 + point.size() &&
-         "Dimensionalities of point and expression don't match!");
-  int64_t value = expr.back();
-  for (unsigned i = 0; i < point.size(); ++i)
-    value += expr[i] * point[i];
-  return value;
-}
-
-/// A point satisfies an equality iff the value of the equality at the
-/// expression is zero, and it satisfies an inequality iff the value of the
-/// inequality at that point is non-negative.
-bool FlatAffineConstraints::containsPoint(ArrayRef<int64_t> point) const {
-  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
-    if (valueAt(getEquality(i), point) != 0)
-      return false;
-  }
-  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
-    if (valueAt(getInequality(i), point) < 0)
-      return false;
-  }
-  return true;
-}
-
-void FlatAffineConstraints::getLocalReprs(
-    std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const {
-  std::vector<SmallVector<int64_t, 8>> dividends(getNumLocalIds());
-  SmallVector<unsigned, 4> denominators(getNumLocalIds());
-  getLocalReprs(dividends, denominators, repr);
-}
-
-void FlatAffineConstraints::getLocalReprs(
-    std::vector<SmallVector<int64_t, 8>> &dividends,
-    SmallVector<unsigned, 4> &denominators) const {
-  std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> repr(
-      getNumLocalIds());
-  getLocalReprs(dividends, denominators, repr);
-}
-
-void FlatAffineConstraints::getLocalReprs(
-    std::vector<SmallVector<int64_t, 8>> &dividends,
-    SmallVector<unsigned, 4> &denominators,
-    std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const {
-
-  repr.resize(getNumLocalIds());
-  dividends.resize(getNumLocalIds());
-  denominators.resize(getNumLocalIds());
-
-  SmallVector<bool, 8> foundRepr(getNumIds(), false);
-  for (unsigned i = 0, e = getNumDimAndSymbolIds(); i < e; ++i)
-    foundRepr[i] = true;
-
-  unsigned divOffset = getNumDimAndSymbolIds();
-  bool changed;
-  do {
-    // Each time changed is true, at end of this iteration, one or more local
-    // vars have been detected as floor divs.
-    changed = false;
-    for (unsigned i = 0, e = getNumLocalIds(); i < e; ++i) {
-      if (!foundRepr[i + divOffset]) {
-        if (auto res = presburger_utils::computeSingleVarRepr(
-                *this, foundRepr, divOffset + i, dividends[i],
-                denominators[i])) {
-          foundRepr[i + divOffset] = true;
-          repr[i] = res;
-          changed = true;
-        }
-      }
-    }
-  } while (changed);
-
-  // Set 0 denominator for identifiers for which no division representation
-  // could be found.
-  for (unsigned i = 0, e = repr.size(); i < e; ++i)
-    if (!repr[i].hasValue())
-      denominators[i] = 0;
-}
-
-/// Tightens inequalities given that we are dealing with integer spaces. This is
-/// analogous to the GCD test but applied to inequalities. The constant term can
-/// be reduced to the preceding multiple of the GCD of the coefficients, i.e.,
-///  64*i - 100 >= 0  =>  64*i - 128 >= 0 (since 'i' is an integer). This is a
-/// fast method - linear in the number of coefficients.
-// Example on how this affects practical cases: consider the scenario:
-// 64*i >= 100, j = 64*i; without a tightening, elimination of i would yield
-// j >= 100 instead of the tighter (exact) j >= 128.
-void FlatAffineConstraints::gcdTightenInequalities() {
-  unsigned numCols = getNumCols();
-  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
-    uint64_t gcd = std::abs(atIneq(i, 0));
-    for (unsigned j = 1; j < numCols - 1; ++j) {
-      gcd = llvm::GreatestCommonDivisor64(gcd, std::abs(atIneq(i, j)));
-    }
-    if (gcd > 0 && gcd != 1) {
-      int64_t gcdI = static_cast<int64_t>(gcd);
-      // Tighten the constant term and normalize the constraint by the GCD.
-      atIneq(i, numCols - 1) = mlir::floorDiv(atIneq(i, numCols - 1), gcdI);
-      for (unsigned j = 0, e = numCols - 1; j < e; ++j)
-        atIneq(i, j) /= gcdI;
-    }
-  }
-}
-
-// Eliminates all identifier variables in column range [posStart, posLimit).
-// Returns the number of variables eliminated.
-unsigned FlatAffineConstraints::gaussianEliminateIds(unsigned posStart,
-                                                     unsigned posLimit) {
-  // Return if identifier positions to eliminate are out of range.
-  assert(posLimit <= numIds);
-  assert(hasConsistentState());
-
-  if (posStart >= posLimit)
-    return 0;
-
-  gcdTightenInequalities();
-
-  unsigned pivotCol = 0;
-  for (pivotCol = posStart; pivotCol < posLimit; ++pivotCol) {
-    // Find a row which has a non-zero coefficient in column 'j'.
-    unsigned pivotRow;
-    if (!findConstraintWithNonZeroAt(pivotCol, /*isEq=*/true, &pivotRow)) {
-      // No pivot row in equalities with non-zero at 'pivotCol'.
-      if (!findConstraintWithNonZeroAt(pivotCol, /*isEq=*/false, &pivotRow)) {
-        // If inequalities are also non-zero in 'pivotCol', it can be
-        // eliminated.
-        continue;
-      }
-      break;
-    }
-
-    // Eliminate identifier at 'pivotCol' from each equality row.
-    for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
-      eliminateFromConstraint(this, i, pivotRow, pivotCol, posStart,
-                              /*isEq=*/true);
-      normalizeConstraintByGCD</*isEq=*/true>(this, i);
-    }
-
-    // Eliminate identifier at 'pivotCol' from each inequality row.
-    for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
-      eliminateFromConstraint(this, i, pivotRow, pivotCol, posStart,
-                              /*isEq=*/false);
-      normalizeConstraintByGCD</*isEq=*/false>(this, i);
-    }
-    removeEquality(pivotRow);
-    gcdTightenInequalities();
-  }
-  // Update position limit based on number eliminated.
-  posLimit = pivotCol;
-  // Remove eliminated columns from all constraints.
-  removeIdRange(posStart, posLimit);
-  return posLimit - posStart;
-}
-
 // Determine whether the identifier at 'pos' (say id_r) can be expressed as
 // modulo of another known identifier (say id_n) w.r.t a constant. For example,
 // if the following constraints hold true:
@@ -1481,278 +869,6 @@ static bool detectAsFloorDiv(const FlatAffineConstraints &cst, unsigned pos,
   return true;
 }
 
-// Fills an inequality row with the value 'val'.
-static inline void fillInequality(FlatAffineConstraints *cst, unsigned r,
-                                  int64_t val) {
-  for (unsigned c = 0, f = cst->getNumCols(); c < f; c++) {
-    cst->atIneq(r, c) = val;
-  }
-}
-
-// Negates an inequality.
-static inline void negateInequality(FlatAffineConstraints *cst, unsigned r) {
-  for (unsigned c = 0, f = cst->getNumCols(); c < f; c++) {
-    cst->atIneq(r, c) = -cst->atIneq(r, c);
-  }
-}
-
-// A more complex check to eliminate redundant inequalities. Uses FourierMotzkin
-// to check if a constraint is redundant.
-void FlatAffineConstraints::removeRedundantInequalities() {
-  SmallVector<bool, 32> redun(getNumInequalities(), false);
-  // To check if an inequality is redundant, we replace the inequality by its
-  // complement (for eg., i - 1 >= 0 by i <= 0), and check if the resulting
-  // system is empty. If it is, the inequality is redundant.
-  FlatAffineConstraints tmpCst(*this);
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
-    // Change the inequality to its complement.
-    negateInequality(&tmpCst, r);
-    tmpCst.atIneq(r, tmpCst.getNumCols() - 1)--;
-    if (tmpCst.isEmpty()) {
-      redun[r] = true;
-      // Zero fill the redundant inequality.
-      fillInequality(this, r, /*val=*/0);
-      fillInequality(&tmpCst, r, /*val=*/0);
-    } else {
-      // Reverse the change (to avoid recreating tmpCst each time).
-      tmpCst.atIneq(r, tmpCst.getNumCols() - 1)++;
-      negateInequality(&tmpCst, r);
-    }
-  }
-
-  // Scan to get rid of all rows marked redundant, in-place.
-  auto copyRow = [&](unsigned src, unsigned dest) {
-    if (src == dest)
-      return;
-    for (unsigned c = 0, e = getNumCols(); c < e; c++) {
-      atIneq(dest, c) = atIneq(src, c);
-    }
-  };
-  unsigned pos = 0;
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
-    if (!redun[r])
-      copyRow(r, pos++);
-  }
-  inequalities.resizeVertically(pos);
-}
-
-// A more complex check to eliminate redundant inequalities and equalities. Uses
-// Simplex to check if a constraint is redundant.
-void FlatAffineConstraints::removeRedundantConstraints() {
-  // First, we run gcdTightenInequalities. This allows us to catch some
-  // constraints which are not redundant when considering rational solutions
-  // but are redundant in terms of integer solutions.
-  gcdTightenInequalities();
-  Simplex simplex(*this);
-  simplex.detectRedundant();
-
-  auto copyInequality = [&](unsigned src, unsigned dest) {
-    if (src == dest)
-      return;
-    for (unsigned c = 0, e = getNumCols(); c < e; c++)
-      atIneq(dest, c) = atIneq(src, c);
-  };
-  unsigned pos = 0;
-  unsigned numIneqs = getNumInequalities();
-  // Scan to get rid of all inequalities marked redundant, in-place. In Simplex,
-  // the first constraints added are the inequalities.
-  for (unsigned r = 0; r < numIneqs; r++) {
-    if (!simplex.isMarkedRedundant(r))
-      copyInequality(r, pos++);
-  }
-  inequalities.resizeVertically(pos);
-
-  // Scan to get rid of all equalities marked redundant, in-place. In Simplex,
-  // after the inequalities, a pair of constraints for each equality is added.
-  // An equality is redundant if both the inequalities in its pair are
-  // redundant.
-  auto copyEquality = [&](unsigned src, unsigned dest) {
-    if (src == dest)
-      return;
-    for (unsigned c = 0, e = getNumCols(); c < e; c++)
-      atEq(dest, c) = atEq(src, c);
-  };
-  pos = 0;
-  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
-    if (!(simplex.isMarkedRedundant(numIneqs + 2 * r) &&
-          simplex.isMarkedRedundant(numIneqs + 2 * r + 1)))
-      copyEquality(r, pos++);
-  }
-  equalities.resizeVertically(pos);
-}
-
-/// Eliminate `pos2^th` local identifier, replacing its every instance with
-/// `pos1^th` local identifier. This function is intended to be used to remove
-/// redundancy when local variables at position `pos1` and `pos2` are restricted
-/// to have the same value.
-static void eliminateRedundantLocalId(FlatAffineConstraints &fac, unsigned pos1,
-                                      unsigned pos2) {
-
-  assert(pos1 < fac.getNumLocalIds() && "Invalid local id position");
-  assert(pos2 < fac.getNumLocalIds() && "Invalid local id position");
-
-  unsigned localOffset = fac.getNumDimAndSymbolIds();
-  pos1 += localOffset;
-  pos2 += localOffset;
-  for (unsigned i = 0, e = fac.getNumInequalities(); i < e; ++i)
-    fac.atIneq(i, pos1) += fac.atIneq(i, pos2);
-  for (unsigned i = 0, e = fac.getNumEqualities(); i < e; ++i)
-    fac.atEq(i, pos1) += fac.atEq(i, pos2);
-  fac.removeId(pos2);
-}
-
-/// Adds additional local ids to the sets such that they both have the union
-/// of the local ids in each set, without changing the set of points that
-/// lie in `this` and `other`.
-///
-/// To detect local ids that always take the same in both sets, each local id is
-/// represented as a floordiv with constant denominator in terms of other ids.
-/// After extracting these divisions, local ids with the same division
-/// representation are considered duplicate and are merged. It is possible that
-/// division representation for some local id cannot be obtained, and thus these
-/// local ids are not considered for detecting duplicates.
-void FlatAffineConstraints::mergeLocalIds(FlatAffineConstraints &other) {
-  assert(getNumDimIds() == other.getNumDimIds() &&
-         "Number of dimension ids should match");
-  assert(getNumSymbolIds() == other.getNumSymbolIds() &&
-         "Number of symbol ids should match");
-
-  FlatAffineConstraints &facA = *this;
-  FlatAffineConstraints &facB = other;
-
-  // Merge local ids of facA and facB without using division information,
-  // i.e. append local ids of `facB` to `facA` and insert local ids of `facA`
-  // to `facB` at start of its local ids.
-  unsigned initLocals = facA.getNumLocalIds();
-  insertLocalId(facA.getNumLocalIds(), facB.getNumLocalIds());
-  facB.insertLocalId(0, initLocals);
-
-  // Get division representations from each FAC.
-  std::vector<SmallVector<int64_t, 8>> divsA, divsB;
-  SmallVector<unsigned, 4> denomsA, denomsB;
-  facA.getLocalReprs(divsA, denomsA);
-  facB.getLocalReprs(divsB, denomsB);
-
-  // Copy division information for facB into `divsA` and `denomsA`, so that
-  // these have the combined division information of both FACs. Since newly
-  // added local variables in facA and facB have no constraints, they will not
-  // have any division representation.
-  std::copy(divsB.begin() + initLocals, divsB.end(),
-            divsA.begin() + initLocals);
-  std::copy(denomsB.begin() + initLocals, denomsB.end(),
-            denomsA.begin() + initLocals);
-
-  // Find and merge duplicate divisions.
-  // TODO: Add division normalization to support divisions that 
diff er by
-  // a constant.
-  // TODO: Add division ordering such that a division representation for local
-  // identifier at position `i` only depends on local identifiers at position <
-  // `i`. This would make sure that all divisions depending on other local
-  // variables that can be merged, are merged.
-  unsigned localOffset = getIdKindOffset(IdKind::Local);
-  for (unsigned i = 0; i < divsA.size(); ++i) {
-    // Check if a division representation exists for the `i^th` local id.
-    if (denomsA[i] == 0)
-      continue;
-    // Check if a division exists which is a duplicate of the division at `i`.
-    for (unsigned j = i + 1; j < divsA.size(); ++j) {
-      // Check if a division representation exists for the `j^th` local id.
-      if (denomsA[j] == 0)
-        continue;
-      // Check if the denominators match.
-      if (denomsA[i] != denomsA[j])
-        continue;
-      // Check if the representations are equal.
-      if (divsA[i] != divsA[j])
-        continue;
-
-      // Merge divisions at position `j` into division at position `i`.
-      eliminateRedundantLocalId(facA, i, j);
-      eliminateRedundantLocalId(facB, i, j);
-      for (unsigned k = 0, g = divsA.size(); k < g; ++k) {
-        SmallVector<int64_t, 8> &div = divsA[k];
-        if (denomsA[k] != 0) {
-          div[localOffset + i] += div[localOffset + j];
-          div.erase(div.begin() + localOffset + j);
-        }
-      }
-
-      divsA.erase(divsA.begin() + j);
-      denomsA.erase(denomsA.begin() + j);
-      // Since `j` can never be zero, we do not need to worry about overflows.
-      --j;
-    }
-  }
-}
-
-/// Removes local variables using equalities. Each equality is checked if it
-/// can be reduced to the form: `e = affine-expr`, where `e` is a local
-/// variable and `affine-expr` is an affine expression not containing `e`.
-/// If an equality satisfies this form, the local variable is replaced in
-/// each constraint and then removed. The equality used to replace this local
-/// variable is also removed.
-void FlatAffineConstraints::removeRedundantLocalVars() {
-  // Normalize the equality constraints to reduce coefficients of local
-  // variables to 1 wherever possible.
-  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i)
-    normalizeConstraintByGCD</*isEq=*/true>(this, i);
-
-  while (true) {
-    unsigned i, e, j, f;
-    for (i = 0, e = getNumEqualities(); i < e; ++i) {
-      // Find a local variable to eliminate using ith equality.
-      for (j = getNumDimAndSymbolIds(), f = getNumIds(); j < f; ++j)
-        if (std::abs(atEq(i, j)) == 1)
-          break;
-
-      // Local variable can be eliminated using ith equality.
-      if (j < f)
-        break;
-    }
-
-    // No equality can be used to eliminate a local variable.
-    if (i == e)
-      break;
-
-    // Use the ith equality to simplify other equalities. If any changes
-    // are made to an equality constraint, it is normalized by GCD.
-    for (unsigned k = 0, t = getNumEqualities(); k < t; ++k) {
-      if (atEq(k, j) != 0) {
-        eliminateFromConstraint(this, k, i, j, j, /*isEq=*/true);
-        normalizeConstraintByGCD</*isEq=*/true>(this, k);
-      }
-    }
-
-    // Use the ith equality to simplify inequalities.
-    for (unsigned k = 0, t = getNumInequalities(); k < t; ++k)
-      eliminateFromConstraint(this, k, i, j, j, /*isEq=*/false);
-
-    // Remove the ith equality and the found local variable.
-    removeId(j);
-    removeEquality(i);
-  }
-}
-
-void FlatAffineConstraints::convertDimToLocal(unsigned dimStart,
-                                              unsigned dimLimit) {
-  assert(dimLimit <= getNumDimIds() && "Invalid dim pos range");
-
-  if (dimStart >= dimLimit)
-    return;
-
-  // Append new local variables corresponding to the dimensions to be converted.
-  unsigned convertCount = dimLimit - dimStart;
-  unsigned newLocalIdStart = getNumIds();
-  appendLocalId(convertCount);
-
-  // Swap the new local variables with dimensions.
-  for (unsigned i = 0; i < convertCount; ++i)
-    swapId(i + dimStart, i + newLocalIdStart);
-
-  // Remove dimensions converted to local variables.
-  removeIdRange(dimStart, dimLimit);
-}
-
 std::pair<AffineMap, AffineMap> FlatAffineConstraints::getLowerAndUpperBound(
     unsigned pos, unsigned offset, unsigned num, unsigned symStartPos,
     ArrayRef<AffineExpr> localExprs, MLIRContext *context) const {
@@ -2179,61 +1295,6 @@ LogicalResult FlatAffineValueConstraints::addSliceBounds(
   return success();
 }
 
-void FlatAffineConstraints::addBound(BoundType type, unsigned pos,
-                                     int64_t value) {
-  assert(pos < getNumCols());
-  if (type == BoundType::EQ) {
-    unsigned row = equalities.appendExtraRow();
-    equalities(row, pos) = 1;
-    equalities(row, getNumCols() - 1) = -value;
-  } else {
-    unsigned row = inequalities.appendExtraRow();
-    inequalities(row, pos) = type == BoundType::LB ? 1 : -1;
-    inequalities(row, getNumCols() - 1) =
-        type == BoundType::LB ? -value : value;
-  }
-}
-
-void FlatAffineConstraints::addBound(BoundType type, ArrayRef<int64_t> expr,
-                                     int64_t value) {
-  assert(type != BoundType::EQ && "EQ not implemented");
-  assert(expr.size() == getNumCols());
-  unsigned row = inequalities.appendExtraRow();
-  for (unsigned i = 0, e = expr.size(); i < e; ++i)
-    inequalities(row, i) = type == BoundType::LB ? expr[i] : -expr[i];
-  inequalities(inequalities.getNumRows() - 1, getNumCols() - 1) +=
-      type == BoundType::LB ? -value : value;
-}
-
-/// Adds a new local identifier as the floordiv of an affine function of other
-/// identifiers, the coefficients of which are provided in 'dividend' and with
-/// respect to a positive constant 'divisor'. Two constraints are added to the
-/// system to capture equivalence with the floordiv.
-///      q = expr floordiv c    <=>   c*q <= expr <= c*q + c - 1.
-void FlatAffineConstraints::addLocalFloorDiv(ArrayRef<int64_t> dividend,
-                                             int64_t divisor) {
-  assert(dividend.size() == getNumCols() && "incorrect dividend size");
-  assert(divisor > 0 && "positive divisor expected");
-
-  appendLocalId();
-
-  // Add two constraints for this new identifier 'q'.
-  SmallVector<int64_t, 8> bound(dividend.size() + 1);
-
-  // dividend - q * divisor >= 0
-  std::copy(dividend.begin(), dividend.begin() + dividend.size() - 1,
-            bound.begin());
-  bound.back() = dividend.back();
-  bound[getNumIds() - 1] = -divisor;
-  addInequality(bound);
-
-  // -dividend +qdivisor * q + divisor - 1 >= 0
-  std::transform(bound.begin(), bound.end(), bound.begin(),
-                 std::negate<int64_t>());
-  bound[bound.size() - 1] += divisor - 1;
-  addInequality(bound);
-}
-
 bool FlatAffineValueConstraints::findId(Value val, unsigned *pos) const {
   unsigned i = 0;
   for (const auto &mayBeId : values) {
@@ -2257,13 +1318,6 @@ void FlatAffineValueConstraints::swapId(unsigned posA, unsigned posB) {
   std::swap(values[posA], values[posB]);
 }
 
-void FlatAffineConstraints::setDimSymbolSeparation(unsigned newSymbolCount) {
-  assert(newSymbolCount <= numDims + numSymbols &&
-         "invalid separation position");
-  numDims = numDims + numSymbols - newSymbolCount;
-  numSymbols = newSymbolCount;
-}
-
 void FlatAffineValueConstraints::addBound(BoundType type, Value val,
                                           int64_t value) {
   unsigned pos;
@@ -2273,289 +1327,6 @@ void FlatAffineValueConstraints::addBound(BoundType type, Value val,
   addBound(type, pos, value);
 }
 
-/// Finds an equality that equates the specified identifier to a constant.
-/// Returns the position of the equality row. If 'symbolic' is set to true,
-/// symbols are also treated like a constant, i.e., an affine function of the
-/// symbols is also treated like a constant. Returns -1 if such an equality
-/// could not be found.
-static int findEqualityToConstant(const FlatAffineConstraints &cst,
-                                  unsigned pos, bool symbolic = false) {
-  assert(pos < cst.getNumIds() && "invalid position");
-  for (unsigned r = 0, e = cst.getNumEqualities(); r < e; r++) {
-    int64_t v = cst.atEq(r, pos);
-    if (v * v != 1)
-      continue;
-    unsigned c;
-    unsigned f = symbolic ? cst.getNumDimIds() : cst.getNumIds();
-    // This checks for zeros in all positions other than 'pos' in [0, f)
-    for (c = 0; c < f; c++) {
-      if (c == pos)
-        continue;
-      if (cst.atEq(r, c) != 0) {
-        // Dependent on another identifier.
-        break;
-      }
-    }
-    if (c == f)
-      // Equality is free of other identifiers.
-      return r;
-  }
-  return -1;
-}
-
-LogicalResult FlatAffineConstraints::constantFoldId(unsigned pos) {
-  assert(pos < getNumIds() && "invalid position");
-  int rowIdx;
-  if ((rowIdx = findEqualityToConstant(*this, pos)) == -1)
-    return failure();
-
-  // atEq(rowIdx, pos) is either -1 or 1.
-  assert(atEq(rowIdx, pos) * atEq(rowIdx, pos) == 1);
-  int64_t constVal = -atEq(rowIdx, getNumCols() - 1) / atEq(rowIdx, pos);
-  setAndEliminate(pos, constVal);
-  return success();
-}
-
-void FlatAffineConstraints::constantFoldIdRange(unsigned pos, unsigned num) {
-  for (unsigned s = pos, t = pos, e = pos + num; s < e; s++) {
-    if (failed(constantFoldId(t)))
-      t++;
-  }
-}
-
-/// Returns a non-negative constant bound on the extent (upper bound - lower
-/// bound) of the specified identifier if it is found to be a constant; returns
-/// None if it's not a constant. This methods treats symbolic identifiers
-/// specially, i.e., it looks for constant 
diff erences between affine
-/// expressions involving only the symbolic identifiers. See comments at
-/// function definition for example. 'lb', if provided, is set to the lower
-/// bound associated with the constant 
diff erence. Note that 'lb' is purely
-/// symbolic and thus will contain the coefficients of the symbolic identifiers
-/// and the constant coefficient.
-//  Egs: 0 <= i <= 15, return 16.
-//       s0 + 2 <= i <= s0 + 17, returns 16. (s0 has to be a symbol)
-//       s0 + s1 + 16 <= d0 <= s0 + s1 + 31, returns 16.
-//       s0 - 7 <= 8*j <= s0 returns 1 with lb = s0, lbDivisor = 8 (since lb =
-//       ceil(s0 - 7 / 8) = floor(s0 / 8)).
-Optional<int64_t> FlatAffineConstraints::getConstantBoundOnDimSize(
-    unsigned pos, SmallVectorImpl<int64_t> *lb, int64_t *boundFloorDivisor,
-    SmallVectorImpl<int64_t> *ub, unsigned *minLbPos,
-    unsigned *minUbPos) const {
-  assert(pos < getNumDimIds() && "Invalid identifier position");
-
-  // Find an equality for 'pos'^th identifier that equates it to some function
-  // of the symbolic identifiers (+ constant).
-  int eqPos = findEqualityToConstant(*this, pos, /*symbolic=*/true);
-  if (eqPos != -1) {
-    auto eq = getEquality(eqPos);
-    // If the equality involves a local var, punt for now.
-    // TODO: this can be handled in the future by using the explicit
-    // representation of the local vars.
-    if (!std::all_of(eq.begin() + getNumDimAndSymbolIds(), eq.end() - 1,
-                     [](int64_t coeff) { return coeff == 0; }))
-      return None;
-
-    // This identifier can only take a single value.
-    if (lb) {
-      // Set lb to that symbolic value.
-      lb->resize(getNumSymbolIds() + 1);
-      if (ub)
-        ub->resize(getNumSymbolIds() + 1);
-      for (unsigned c = 0, f = getNumSymbolIds() + 1; c < f; c++) {
-        int64_t v = atEq(eqPos, pos);
-        // atEq(eqRow, pos) is either -1 or 1.
-        assert(v * v == 1);
-        (*lb)[c] = v < 0 ? atEq(eqPos, getNumDimIds() + c) / -v
-                         : -atEq(eqPos, getNumDimIds() + c) / v;
-        // Since this is an equality, ub = lb.
-        if (ub)
-          (*ub)[c] = (*lb)[c];
-      }
-      assert(boundFloorDivisor &&
-             "both lb and divisor or none should be provided");
-      *boundFloorDivisor = 1;
-    }
-    if (minLbPos)
-      *minLbPos = eqPos;
-    if (minUbPos)
-      *minUbPos = eqPos;
-    return 1;
-  }
-
-  // Check if the identifier appears at all in any of the inequalities.
-  unsigned r, e;
-  for (r = 0, e = getNumInequalities(); r < e; r++) {
-    if (atIneq(r, pos) != 0)
-      break;
-  }
-  if (r == e)
-    // If it doesn't, there isn't a bound on it.
-    return None;
-
-  // Positions of constraints that are lower/upper bounds on the variable.
-  SmallVector<unsigned, 4> lbIndices, ubIndices;
-
-  // Gather all symbolic lower bounds and upper bounds of the variable, i.e.,
-  // the bounds can only involve symbolic (and local) identifiers. Since the
-  // canonical form c_1*x_1 + c_2*x_2 + ... + c_0 >= 0, a constraint is a lower
-  // bound for x_i if c_i >= 1, and an upper bound if c_i <= -1.
-  getLowerAndUpperBoundIndices(pos, &lbIndices, &ubIndices,
-                               /*eqIndices=*/nullptr, /*offset=*/0,
-                               /*num=*/getNumDimIds());
-
-  Optional<int64_t> minDiff = None;
-  unsigned minLbPosition = 0, minUbPosition = 0;
-  for (auto ubPos : ubIndices) {
-    for (auto lbPos : lbIndices) {
-      // Look for a lower bound and an upper bound that only 
diff er by a
-      // constant, i.e., pairs of the form  0 <= c_pos - f(c_i's) <= 
diff Const.
-      // For example, if ii is the pos^th variable, we are looking for
-      // constraints like ii >= i, ii <= ii + 50, 50 being the 
diff erence. The
-      // minimum among all such constant 
diff erences is kept since that's the
-      // constant bounding the extent of the pos^th variable.
-      unsigned j, e;
-      for (j = 0, e = getNumCols() - 1; j < e; j++)
-        if (atIneq(ubPos, j) != -atIneq(lbPos, j)) {
-          break;
-        }
-      if (j < getNumCols() - 1)
-        continue;
-      int64_t 
diff  = ceilDiv(atIneq(ubPos, getNumCols() - 1) +
-                                 atIneq(lbPos, getNumCols() - 1) + 1,
-                             atIneq(lbPos, pos));
-      // This bound is non-negative by definition.
-      
diff  = std::max<int64_t>(
diff , 0);
-      if (minDiff == None || 
diff  < minDiff) {
-        minDiff = 
diff ;
-        minLbPosition = lbPos;
-        minUbPosition = ubPos;
-      }
-    }
-  }
-  if (lb && minDiff.hasValue()) {
-    // Set lb to the symbolic lower bound.
-    lb->resize(getNumSymbolIds() + 1);
-    if (ub)
-      ub->resize(getNumSymbolIds() + 1);
-    // The lower bound is the ceildiv of the lb constraint over the coefficient
-    // of the variable at 'pos'. We express the ceildiv equivalently as a floor
-    // for uniformity. For eg., if the lower bound constraint was: 32*d0 - N +
-    // 31 >= 0, the lower bound for d0 is ceil(N - 31, 32), i.e., floor(N, 32).
-    *boundFloorDivisor = atIneq(minLbPosition, pos);
-    assert(*boundFloorDivisor == -atIneq(minUbPosition, pos));
-    for (unsigned c = 0, e = getNumSymbolIds() + 1; c < e; c++) {
-      (*lb)[c] = -atIneq(minLbPosition, getNumDimIds() + c);
-    }
-    if (ub) {
-      for (unsigned c = 0, e = getNumSymbolIds() + 1; c < e; c++)
-        (*ub)[c] = atIneq(minUbPosition, getNumDimIds() + c);
-    }
-    // The lower bound leads to a ceildiv while the upper bound is a floordiv
-    // whenever the coefficient at pos != 1. ceildiv (val / d) = floordiv (val +
-    // d - 1 / d); hence, the addition of 'atIneq(minLbPosition, pos) - 1' to
-    // the constant term for the lower bound.
-    (*lb)[getNumSymbolIds()] += atIneq(minLbPosition, pos) - 1;
-  }
-  if (minLbPos)
-    *minLbPos = minLbPosition;
-  if (minUbPos)
-    *minUbPos = minUbPosition;
-  return minDiff;
-}
-
-template <bool isLower>
-Optional<int64_t>
-FlatAffineConstraints::computeConstantLowerOrUpperBound(unsigned pos) {
-  assert(pos < getNumIds() && "invalid position");
-  // Project to 'pos'.
-  projectOut(0, pos);
-  projectOut(1, getNumIds() - 1);
-  // Check if there's an equality equating the '0'^th identifier to a constant.
-  int eqRowIdx = findEqualityToConstant(*this, 0, /*symbolic=*/false);
-  if (eqRowIdx != -1)
-    // atEq(rowIdx, 0) is either -1 or 1.
-    return -atEq(eqRowIdx, getNumCols() - 1) / atEq(eqRowIdx, 0);
-
-  // Check if the identifier appears at all in any of the inequalities.
-  unsigned r, e;
-  for (r = 0, e = getNumInequalities(); r < e; r++) {
-    if (atIneq(r, 0) != 0)
-      break;
-  }
-  if (r == e)
-    // If it doesn't, there isn't a bound on it.
-    return None;
-
-  Optional<int64_t> minOrMaxConst = None;
-
-  // Take the max across all const lower bounds (or min across all constant
-  // upper bounds).
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
-    if (isLower) {
-      if (atIneq(r, 0) <= 0)
-        // Not a lower bound.
-        continue;
-    } else if (atIneq(r, 0) >= 0) {
-      // Not an upper bound.
-      continue;
-    }
-    unsigned c, f;
-    for (c = 0, f = getNumCols() - 1; c < f; c++)
-      if (c != 0 && atIneq(r, c) != 0)
-        break;
-    if (c < getNumCols() - 1)
-      // Not a constant bound.
-      continue;
-
-    int64_t boundConst =
-        isLower ? mlir::ceilDiv(-atIneq(r, getNumCols() - 1), atIneq(r, 0))
-                : mlir::floorDiv(atIneq(r, getNumCols() - 1), -atIneq(r, 0));
-    if (isLower) {
-      if (minOrMaxConst == None || boundConst > minOrMaxConst)
-        minOrMaxConst = boundConst;
-    } else {
-      if (minOrMaxConst == None || boundConst < minOrMaxConst)
-        minOrMaxConst = boundConst;
-    }
-  }
-  return minOrMaxConst;
-}
-
-Optional<int64_t> FlatAffineConstraints::getConstantBound(BoundType type,
-                                                          unsigned pos) const {
-  assert(type != BoundType::EQ && "EQ not implemented");
-  FlatAffineConstraints tmpCst(*this);
-  if (type == BoundType::LB)
-    return tmpCst.computeConstantLowerOrUpperBound</*isLower=*/true>(pos);
-  return tmpCst.computeConstantLowerOrUpperBound</*isLower=*/false>(pos);
-}
-
-// A simple (naive and conservative) check for hyper-rectangularity.
-bool FlatAffineConstraints::isHyperRectangular(unsigned pos,
-                                               unsigned num) const {
-  assert(pos < getNumCols() - 1);
-  // Check for two non-zero coefficients in the range [pos, pos + sum).
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
-    unsigned sum = 0;
-    for (unsigned c = pos; c < pos + num; c++) {
-      if (atIneq(r, c) != 0)
-        sum++;
-    }
-    if (sum > 1)
-      return false;
-  }
-  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
-    unsigned sum = 0;
-    for (unsigned c = pos; c < pos + num; c++) {
-      if (atEq(r, c) != 0)
-        sum++;
-    }
-    if (sum > 1)
-      return false;
-  }
-  return true;
-}
-
 void FlatAffineConstraints::printSpace(raw_ostream &os) const {
   IntegerPolyhedron::printSpace(os);
   os << "(";
@@ -2572,77 +1343,6 @@ void FlatAffineConstraints::printSpace(raw_ostream &os) const {
   os << " const)\n";
 }
 
-/// Removes duplicate constraints, trivially true constraints, and constraints
-/// that can be detected as redundant as a result of 
diff ering only in their
-/// constant term part. A constraint of the form <non-negative constant> >= 0 is
-/// considered trivially true.
-//  Uses a DenseSet to hash and detect duplicates followed by a linear scan to
-//  remove duplicates in place.
-void FlatAffineConstraints::removeTrivialRedundancy() {
-  gcdTightenInequalities();
-  normalizeConstraintsByGCD();
-
-  // A map used to detect redundancy stemming from constraints that only 
diff er
-  // in their constant term. The value stored is <row position, const term>
-  // for a given row.
-  SmallDenseMap<ArrayRef<int64_t>, std::pair<unsigned, int64_t>>
-      rowsWithoutConstTerm;
-  // To unique rows.
-  SmallDenseSet<ArrayRef<int64_t>, 8> rowSet;
-
-  // Check if constraint is of the form <non-negative-constant> >= 0.
-  auto isTriviallyValid = [&](unsigned r) -> bool {
-    for (unsigned c = 0, e = getNumCols() - 1; c < e; c++) {
-      if (atIneq(r, c) != 0)
-        return false;
-    }
-    return atIneq(r, getNumCols() - 1) >= 0;
-  };
-
-  // Detect and mark redundant constraints.
-  SmallVector<bool, 256> redunIneq(getNumInequalities(), false);
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
-    int64_t *rowStart = &inequalities(r, 0);
-    auto row = ArrayRef<int64_t>(rowStart, getNumCols());
-    if (isTriviallyValid(r) || !rowSet.insert(row).second) {
-      redunIneq[r] = true;
-      continue;
-    }
-
-    // Among constraints that only 
diff er in the constant term part, mark
-    // everything other than the one with the smallest constant term redundant.
-    // (eg: among i - 16j - 5 >= 0, i - 16j - 1 >=0, i - 16j - 7 >= 0, the
-    // former two are redundant).
-    int64_t constTerm = atIneq(r, getNumCols() - 1);
-    auto rowWithoutConstTerm = ArrayRef<int64_t>(rowStart, getNumCols() - 1);
-    const auto &ret =
-        rowsWithoutConstTerm.insert({rowWithoutConstTerm, {r, constTerm}});
-    if (!ret.second) {
-      // Check if the other constraint has a higher constant term.
-      auto &val = ret.first->second;
-      if (val.second > constTerm) {
-        // The stored row is redundant. Mark it so, and update with this one.
-        redunIneq[val.first] = true;
-        val = {r, constTerm};
-      } else {
-        // The one stored makes this one redundant.
-        redunIneq[r] = true;
-      }
-    }
-  }
-
-  // Scan to get rid of all rows marked redundant, in-place.
-  unsigned pos = 0;
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++)
-    if (!redunIneq[r])
-      inequalities.copyRow(r, pos++);
-
-  inequalities.resizeVertically(pos);
-
-  // TODO: consider doing this for equalities as well, but probably not worth
-  // the savings.
-}
-
 void FlatAffineConstraints::clearAndCopyFrom(const IntegerPolyhedron &other) {
   if (auto *otherValueSet = dyn_cast<const FlatAffineValueConstraints>(&other))
     assert(!otherValueSet->hasValues() &&
@@ -2674,232 +1374,6 @@ void FlatAffineValueConstraints::clearAndCopyFrom(
   values.resize(numIds, None);
 }
 
-static std::pair<unsigned, unsigned>
-getNewNumDimsSymbols(unsigned pos, const FlatAffineConstraints &cst) {
-  unsigned numDims = cst.getNumDimIds();
-  unsigned numSymbols = cst.getNumSymbolIds();
-  unsigned newNumDims, newNumSymbols;
-  if (pos < numDims) {
-    newNumDims = numDims - 1;
-    newNumSymbols = numSymbols;
-  } else if (pos < numDims + numSymbols) {
-    assert(numSymbols >= 1);
-    newNumDims = numDims;
-    newNumSymbols = numSymbols - 1;
-  } else {
-    newNumDims = numDims;
-    newNumSymbols = numSymbols;
-  }
-  return {newNumDims, newNumSymbols};
-}
-
-#undef DEBUG_TYPE
-#define DEBUG_TYPE "fm"
-
-/// Eliminates identifier at the specified position using Fourier-Motzkin
-/// variable elimination. This technique is exact for rational spaces but
-/// conservative (in "rare" cases) for integer spaces. The operation corresponds
-/// to a projection operation yielding the (convex) set of integer points
-/// contained in the rational shadow of the set. An emptiness test that relies
-/// on this method will guarantee emptiness, i.e., it disproves the existence of
-/// a solution if it says it's empty.
-/// If a non-null isResultIntegerExact is passed, it is set to true if the
-/// result is also integer exact. If it's set to false, the obtained solution
-/// *may* not be exact, i.e., it may contain integer points that do not have an
-/// integer pre-image in the original set.
-///
-/// Eg:
-/// j >= 0, j <= i + 1
-/// i >= 0, i <= N + 1
-/// Eliminating i yields,
-///   j >= 0, 0 <= N + 1, j - 1 <= N + 1
-///
-/// If darkShadow = true, this method computes the dark shadow on elimination;
-/// the dark shadow is a convex integer subset of the exact integer shadow. A
-/// non-empty dark shadow proves the existence of an integer solution. The
-/// elimination in such a case could however be an under-approximation, and thus
-/// should not be used for scanning sets or used by itself for dependence
-/// checking.
-///
-/// Eg: 2-d set, * represents grid points, 'o' represents a point in the set.
-///            ^
-///            |
-///            | * * * * o o
-///         i  | * * o o o o
-///            | o * * * * *
-///            --------------->
-///                 j ->
-///
-/// Eliminating i from this system (projecting on the j dimension):
-/// rational shadow / integer light shadow:  1 <= j <= 6
-/// dark shadow:                             3 <= j <= 6
-/// exact integer shadow:                    j = 1 \union  3 <= j <= 6
-/// holes/splinters:                         j = 2
-///
-/// darkShadow = false, isResultIntegerExact = nullptr are default values.
-// TODO: a slight modification to yield dark shadow version of FM (tightened),
-// which can prove the existence of a solution if there is one.
-void FlatAffineConstraints::fourierMotzkinEliminate(
-    unsigned pos, bool darkShadow, bool *isResultIntegerExact) {
-  LLVM_DEBUG(llvm::dbgs() << "FM input (eliminate pos " << pos << "):\n");
-  LLVM_DEBUG(dump());
-  assert(pos < getNumIds() && "invalid position");
-  assert(hasConsistentState());
-
-  // Check if this identifier can be eliminated through a substitution.
-  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
-    if (atEq(r, pos) != 0) {
-      // Use Gaussian elimination here (since we have an equality).
-      LogicalResult ret = gaussianEliminateId(pos);
-      (void)ret;
-      assert(succeeded(ret) && "Gaussian elimination guaranteed to succeed");
-      LLVM_DEBUG(llvm::dbgs() << "FM output (through Gaussian elimination):\n");
-      LLVM_DEBUG(dump());
-      return;
-    }
-  }
-
-  // A fast linear time tightening.
-  gcdTightenInequalities();
-
-  // Check if the identifier appears at all in any of the inequalities.
-  unsigned r, e;
-  for (r = 0, e = getNumInequalities(); r < e; r++) {
-    if (atIneq(r, pos) != 0)
-      break;
-  }
-  if (r == getNumInequalities()) {
-    // If it doesn't appear, just remove the column and return.
-    // TODO: refactor removeColumns to use it from here.
-    removeId(pos);
-    LLVM_DEBUG(llvm::dbgs() << "FM output:\n");
-    LLVM_DEBUG(dump());
-    return;
-  }
-
-  // Positions of constraints that are lower bounds on the variable.
-  SmallVector<unsigned, 4> lbIndices;
-  // Positions of constraints that are lower bounds on the variable.
-  SmallVector<unsigned, 4> ubIndices;
-  // Positions of constraints that do not involve the variable.
-  std::vector<unsigned> nbIndices;
-  nbIndices.reserve(getNumInequalities());
-
-  // Gather all lower bounds and upper bounds of the variable. Since the
-  // canonical form c_1*x_1 + c_2*x_2 + ... + c_0 >= 0, a constraint is a lower
-  // bound for x_i if c_i >= 1, and an upper bound if c_i <= -1.
-  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
-    if (atIneq(r, pos) == 0) {
-      // Id does not appear in bound.
-      nbIndices.push_back(r);
-    } else if (atIneq(r, pos) >= 1) {
-      // Lower bound.
-      lbIndices.push_back(r);
-    } else {
-      // Upper bound.
-      ubIndices.push_back(r);
-    }
-  }
-
-  // Set the number of dimensions, symbols in the resulting system.
-  const auto &dimsSymbols = getNewNumDimsSymbols(pos, *this);
-  unsigned newNumDims = dimsSymbols.first;
-  unsigned newNumSymbols = dimsSymbols.second;
-
-  /// Create the new system which has one identifier less.
-  FlatAffineConstraints newFac(
-      lbIndices.size() * ubIndices.size() + nbIndices.size(),
-      getNumEqualities(), getNumCols() - 1, newNumDims, newNumSymbols,
-      /*numLocals=*/getNumIds() - 1 - newNumDims - newNumSymbols);
-
-  // This will be used to check if the elimination was integer exact.
-  unsigned lcmProducts = 1;
-
-  // Let x be the variable we are eliminating.
-  // For each lower bound, lb <= c_l*x, and each upper bound c_u*x <= ub, (note
-  // that c_l, c_u >= 1) we have:
-  // lb*lcm(c_l, c_u)/c_l <= lcm(c_l, c_u)*x <= ub*lcm(c_l, c_u)/c_u
-  // We thus generate a constraint:
-  // lcm(c_l, c_u)/c_l*lb <= lcm(c_l, c_u)/c_u*ub.
-  // Note if c_l = c_u = 1, all integer points captured by the resulting
-  // constraint correspond to integer points in the original system (i.e., they
-  // have integer pre-images). Hence, if the lcm's are all 1, the elimination is
-  // integer exact.
-  for (auto ubPos : ubIndices) {
-    for (auto lbPos : lbIndices) {
-      SmallVector<int64_t, 4> ineq;
-      ineq.reserve(newFac.getNumCols());
-      int64_t lbCoeff = atIneq(lbPos, pos);
-      // Note that in the comments above, ubCoeff is the negation of the
-      // coefficient in the canonical form as the view taken here is that of the
-      // term being moved to the other size of '>='.
-      int64_t ubCoeff = -atIneq(ubPos, pos);
-      // TODO: refactor this loop to avoid all branches inside.
-      for (unsigned l = 0, e = getNumCols(); l < e; l++) {
-        if (l == pos)
-          continue;
-        assert(lbCoeff >= 1 && ubCoeff >= 1 && "bounds wrongly identified");
-        int64_t lcm = mlir::lcm(lbCoeff, ubCoeff);
-        ineq.push_back(atIneq(ubPos, l) * (lcm / ubCoeff) +
-                       atIneq(lbPos, l) * (lcm / lbCoeff));
-        lcmProducts *= lcm;
-      }
-      if (darkShadow) {
-        // The dark shadow is a convex subset of the exact integer shadow. If
-        // there is a point here, it proves the existence of a solution.
-        ineq[ineq.size() - 1] += lbCoeff * ubCoeff - lbCoeff - ubCoeff + 1;
-      }
-      // TODO: we need to have a way to add inequalities in-place in
-      // FlatAffineConstraints instead of creating and copying over.
-      newFac.addInequality(ineq);
-    }
-  }
-
-  LLVM_DEBUG(llvm::dbgs() << "FM isResultIntegerExact: " << (lcmProducts == 1)
-                          << "\n");
-  if (lcmProducts == 1 && isResultIntegerExact)
-    *isResultIntegerExact = true;
-
-  // Copy over the constraints not involving this variable.
-  for (auto nbPos : nbIndices) {
-    SmallVector<int64_t, 4> ineq;
-    ineq.reserve(getNumCols() - 1);
-    for (unsigned l = 0, e = getNumCols(); l < e; l++) {
-      if (l == pos)
-        continue;
-      ineq.push_back(atIneq(nbPos, l));
-    }
-    newFac.addInequality(ineq);
-  }
-
-  assert(newFac.getNumConstraints() ==
-         lbIndices.size() * ubIndices.size() + nbIndices.size());
-
-  // Copy over the equalities.
-  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
-    SmallVector<int64_t, 4> eq;
-    eq.reserve(newFac.getNumCols());
-    for (unsigned l = 0, e = getNumCols(); l < e; l++) {
-      if (l == pos)
-        continue;
-      eq.push_back(atEq(r, l));
-    }
-    newFac.addEquality(eq);
-  }
-
-  // GCD tightening and normalization allows detection of more trivially
-  // redundant constraints.
-  newFac.gcdTightenInequalities();
-  newFac.normalizeConstraintsByGCD();
-  newFac.removeTrivialRedundancy();
-  clearAndCopyFrom(newFac);
-  LLVM_DEBUG(llvm::dbgs() << "FM output:\n");
-  LLVM_DEBUG(dump());
-}
-
-#undef DEBUG_TYPE
-#define DEBUG_TYPE "affine-structures"
-
 void FlatAffineValueConstraints::fourierMotzkinEliminate(
     unsigned pos, bool darkShadow, bool *isResultIntegerExact) {
   SmallVector<Optional<Value>, 8> newVals;
@@ -2913,41 +1387,6 @@ void FlatAffineValueConstraints::fourierMotzkinEliminate(
   assert(values.size() == getNumIds());
 }
 
-void FlatAffineConstraints::projectOut(unsigned pos, unsigned num) {
-  if (num == 0)
-    return;
-
-  // 'pos' can be at most getNumCols() - 2 if num > 0.
-  assert((getNumCols() < 2 || pos <= getNumCols() - 2) && "invalid position");
-  assert(pos + num < getNumCols() && "invalid range");
-
-  // Eliminate as many identifiers as possible using Gaussian elimination.
-  unsigned currentPos = pos;
-  unsigned numToEliminate = num;
-  unsigned numGaussianEliminated = 0;
-
-  while (currentPos < getNumIds()) {
-    unsigned curNumEliminated =
-        gaussianEliminateIds(currentPos, currentPos + numToEliminate);
-    ++currentPos;
-    numToEliminate -= curNumEliminated + 1;
-    numGaussianEliminated += curNumEliminated;
-  }
-
-  // Eliminate the remaining using Fourier-Motzkin.
-  for (unsigned i = 0; i < num - numGaussianEliminated; i++) {
-    unsigned numToEliminate = num - numGaussianEliminated - i;
-    fourierMotzkinEliminate(
-        getBestIdToEliminate(*this, pos, pos + numToEliminate));
-  }
-
-  // Fast/trivial simplifications.
-  gcdTightenInequalities();
-  // Normalize constraints after tightening since the latter impacts this, but
-  // not the other way round.
-  normalizeConstraintsByGCD();
-}
-
 void FlatAffineValueConstraints::projectOut(Value val) {
   unsigned pos;
   bool ret = findId(val, &pos);
@@ -2956,167 +1395,6 @@ void FlatAffineValueConstraints::projectOut(Value val) {
   fourierMotzkinEliminate(pos);
 }
 
-namespace {
-
-enum BoundCmpResult { Greater, Less, Equal, Unknown };
-
-/// Compares two affine bounds whose coefficients are provided in 'first' and
-/// 'second'. The last coefficient is the constant term.
-static BoundCmpResult compareBounds(ArrayRef<int64_t> a, ArrayRef<int64_t> b) {
-  assert(a.size() == b.size());
-
-  // For the bounds to be comparable, their corresponding identifier
-  // coefficients should be equal; the constant terms are then compared to
-  // determine less/greater/equal.
-
-  if (!std::equal(a.begin(), a.end() - 1, b.begin()))
-    return Unknown;
-
-  if (a.back() == b.back())
-    return Equal;
-
-  return a.back() < b.back() ? Less : Greater;
-}
-} // namespace
-
-// Returns constraints that are common to both A & B.
-static void getCommonConstraints(const FlatAffineConstraints &a,
-                                 const FlatAffineConstraints &b,
-                                 FlatAffineConstraints &c) {
-  c.reset(a.getNumDimIds(), a.getNumSymbolIds(), a.getNumLocalIds());
-  // a naive O(n^2) check should be enough here given the input sizes.
-  for (unsigned r = 0, e = a.getNumInequalities(); r < e; ++r) {
-    for (unsigned s = 0, f = b.getNumInequalities(); s < f; ++s) {
-      if (a.getInequality(r) == b.getInequality(s)) {
-        c.addInequality(a.getInequality(r));
-        break;
-      }
-    }
-  }
-  for (unsigned r = 0, e = a.getNumEqualities(); r < e; ++r) {
-    for (unsigned s = 0, f = b.getNumEqualities(); s < f; ++s) {
-      if (a.getEquality(r) == b.getEquality(s)) {
-        c.addEquality(a.getEquality(r));
-        break;
-      }
-    }
-  }
-}
-
-// Computes the bounding box with respect to 'other' by finding the min of the
-// lower bounds and the max of the upper bounds along each of the dimensions.
-LogicalResult
-FlatAffineConstraints::unionBoundingBox(const FlatAffineConstraints &otherCst) {
-  assert(otherCst.getNumDimIds() == numDims && "dims mismatch");
-  assert(otherCst.getNumLocalIds() == 0 && "local ids not supported here");
-  assert(getNumLocalIds() == 0 && "local ids not supported yet here");
-
-  // Get the constraints common to both systems; these will be added as is to
-  // the union.
-  FlatAffineConstraints commonCst;
-  getCommonConstraints(*this, otherCst, commonCst);
-
-  std::vector<SmallVector<int64_t, 8>> boundingLbs;
-  std::vector<SmallVector<int64_t, 8>> boundingUbs;
-  boundingLbs.reserve(2 * getNumDimIds());
-  boundingUbs.reserve(2 * getNumDimIds());
-
-  // To hold lower and upper bounds for each dimension.
-  SmallVector<int64_t, 4> lb, otherLb, ub, otherUb;
-  // To compute min of lower bounds and max of upper bounds for each dimension.
-  SmallVector<int64_t, 4> minLb(getNumSymbolIds() + 1);
-  SmallVector<int64_t, 4> maxUb(getNumSymbolIds() + 1);
-  // To compute final new lower and upper bounds for the union.
-  SmallVector<int64_t, 8> newLb(getNumCols()), newUb(getNumCols());
-
-  int64_t lbFloorDivisor, otherLbFloorDivisor;
-  for (unsigned d = 0, e = getNumDimIds(); d < e; ++d) {
-    auto extent = getConstantBoundOnDimSize(d, &lb, &lbFloorDivisor, &ub);
-    if (!extent.hasValue())
-      // TODO: symbolic extents when necessary.
-      // TODO: handle union if a dimension is unbounded.
-      return failure();
-
-    auto otherExtent = otherCst.getConstantBoundOnDimSize(
-        d, &otherLb, &otherLbFloorDivisor, &otherUb);
-    if (!otherExtent.hasValue() || lbFloorDivisor != otherLbFloorDivisor)
-      // TODO: symbolic extents when necessary.
-      return failure();
-
-    assert(lbFloorDivisor > 0 && "divisor always expected to be positive");
-
-    auto res = compareBounds(lb, otherLb);
-    // Identify min.
-    if (res == BoundCmpResult::Less || res == BoundCmpResult::Equal) {
-      minLb = lb;
-      // Since the divisor is for a floordiv, we need to convert to ceildiv,
-      // i.e., i >= expr floordiv div <=> i >= (expr - div + 1) ceildiv div <=>
-      // div * i >= expr - div + 1.
-      minLb.back() -= lbFloorDivisor - 1;
-    } else if (res == BoundCmpResult::Greater) {
-      minLb = otherLb;
-      minLb.back() -= otherLbFloorDivisor - 1;
-    } else {
-      // Uncomparable - check for constant lower/upper bounds.
-      auto constLb = getConstantBound(BoundType::LB, d);
-      auto constOtherLb = otherCst.getConstantBound(BoundType::LB, d);
-      if (!constLb.hasValue() || !constOtherLb.hasValue())
-        return failure();
-      std::fill(minLb.begin(), minLb.end(), 0);
-      minLb.back() = std::min(constLb.getValue(), constOtherLb.getValue());
-    }
-
-    // Do the same for ub's but max of upper bounds. Identify max.
-    auto uRes = compareBounds(ub, otherUb);
-    if (uRes == BoundCmpResult::Greater || uRes == BoundCmpResult::Equal) {
-      maxUb = ub;
-    } else if (uRes == BoundCmpResult::Less) {
-      maxUb = otherUb;
-    } else {
-      // Uncomparable - check for constant lower/upper bounds.
-      auto constUb = getConstantBound(BoundType::UB, d);
-      auto constOtherUb = otherCst.getConstantBound(BoundType::UB, d);
-      if (!constUb.hasValue() || !constOtherUb.hasValue())
-        return failure();
-      std::fill(maxUb.begin(), maxUb.end(), 0);
-      maxUb.back() = std::max(constUb.getValue(), constOtherUb.getValue());
-    }
-
-    std::fill(newLb.begin(), newLb.end(), 0);
-    std::fill(newUb.begin(), newUb.end(), 0);
-
-    // The divisor for lb, ub, otherLb, otherUb at this point is lbDivisor,
-    // and so it's the divisor for newLb and newUb as well.
-    newLb[d] = lbFloorDivisor;
-    newUb[d] = -lbFloorDivisor;
-    // Copy over the symbolic part + constant term.
-    std::copy(minLb.begin(), minLb.end(), newLb.begin() + getNumDimIds());
-    std::transform(newLb.begin() + getNumDimIds(), newLb.end(),
-                   newLb.begin() + getNumDimIds(), std::negate<int64_t>());
-    std::copy(maxUb.begin(), maxUb.end(), newUb.begin() + getNumDimIds());
-
-    boundingLbs.push_back(newLb);
-    boundingUbs.push_back(newUb);
-  }
-
-  // Clear all constraints and add the lower/upper bounds for the bounding box.
-  clearConstraints();
-  for (unsigned d = 0, e = getNumDimIds(); d < e; ++d) {
-    addInequality(boundingLbs[d]);
-    addInequality(boundingUbs[d]);
-  }
-
-  // Add the constraints that were common to both systems.
-  append(commonCst);
-  removeTrivialRedundancy();
-
-  // TODO: copy over pure symbolic constraints from this and 'other' over to the
-  // union (since the above are just the union along dimensions); we shouldn't
-  // be discarding any other constraints on the symbols.
-
-  return success();
-}
-
 LogicalResult FlatAffineValueConstraints::unionBoundingBox(
     const FlatAffineValueConstraints &otherCst) {
   assert(otherCst.getNumDimIds() == numDims && "dims mismatch");
@@ -3215,12 +1493,6 @@ void FlatAffineValueConstraints::getIneqAsAffineValueMap(
   vmap.reset(AffineMap::get(numDims - 1, numSyms, boundExpr), operands);
 }
 
-bool FlatAffineConstraints::isColZero(unsigned pos) const {
-  unsigned rowPos;
-  return !findConstraintWithNonZeroAt(pos, /*isEq=*/false, &rowPos) &&
-         !findConstraintWithNonZeroAt(pos, /*isEq=*/true, &rowPos);
-}
-
 IntegerSet FlatAffineConstraints::getAsIntegerSet(MLIRContext *context) const {
   if (getNumConstraints() == 0)
     // Return universal set (always true): 0 == 0.
@@ -3274,55 +1546,6 @@ IntegerSet FlatAffineConstraints::getAsIntegerSet(MLIRContext *context) const {
   return IntegerSet::get(numDims, numSyms, exprs, eqFlags);
 }
 
-/// Find positions of inequalities and equalities that do not have a coefficient
-/// for [pos, pos + num) identifiers.
-static void getIndependentConstraints(const FlatAffineConstraints &cst,
-                                      unsigned pos, unsigned num,
-                                      SmallVectorImpl<unsigned> &nbIneqIndices,
-                                      SmallVectorImpl<unsigned> &nbEqIndices) {
-  assert(pos < cst.getNumIds() && "invalid start position");
-  assert(pos + num <= cst.getNumIds() && "invalid limit");
-
-  for (unsigned r = 0, e = cst.getNumInequalities(); r < e; r++) {
-    // The bounds are to be independent of [offset, offset + num) columns.
-    unsigned c;
-    for (c = pos; c < pos + num; ++c) {
-      if (cst.atIneq(r, c) != 0)
-        break;
-    }
-    if (c == pos + num)
-      nbIneqIndices.push_back(r);
-  }
-
-  for (unsigned r = 0, e = cst.getNumEqualities(); r < e; r++) {
-    // The bounds are to be independent of [offset, offset + num) columns.
-    unsigned c;
-    for (c = pos; c < pos + num; ++c) {
-      if (cst.atEq(r, c) != 0)
-        break;
-    }
-    if (c == pos + num)
-      nbEqIndices.push_back(r);
-  }
-}
-
-void FlatAffineConstraints::removeIndependentConstraints(unsigned pos,
-                                                         unsigned num) {
-  assert(pos + num <= getNumIds() && "invalid range");
-
-  // Remove constraints that are independent of these identifiers.
-  SmallVector<unsigned, 4> nbIneqIndices, nbEqIndices;
-  getIndependentConstraints(*this, /*pos=*/0, num, nbIneqIndices, nbEqIndices);
-
-  // Iterate in reverse so that indices don't have to be updated.
-  // TODO: This method can be made more efficient (because removal of each
-  // inequality leads to much shifting/copying in the underlying buffer).
-  for (auto nbIndex : llvm::reverse(nbIneqIndices))
-    removeInequality(nbIndex);
-  for (auto nbIndex : llvm::reverse(nbEqIndices))
-    removeEquality(nbIndex);
-}
-
 AffineMap mlir::alignAffineMapWithValues(AffineMap map, ValueRange operands,
                                          ValueRange dims, ValueRange syms,
                                          SmallVector<Value> *newSyms) {

diff  --git a/mlir/lib/Analysis/Presburger/IntegerPolyhedron.cpp b/mlir/lib/Analysis/Presburger/IntegerPolyhedron.cpp
index 06b39cbc0f8b..520d17a31bab 100644
--- a/mlir/lib/Analysis/Presburger/IntegerPolyhedron.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerPolyhedron.cpp
@@ -11,8 +11,18 @@
 //===----------------------------------------------------------------------===//
 
 #include "mlir/Analysis/Presburger/IntegerPolyhedron.h"
+#include "mlir/Analysis/Presburger/LinearTransform.h"
+#include "mlir/Analysis/Presburger/Simplex.h"
+#include "mlir/Analysis/Presburger/Utils.h"
+#include "llvm/ADT/DenseMap.h"
+#include "llvm/ADT/DenseSet.h"
+#include "llvm/Support/Debug.h"
+
+#define DEBUG_TYPE "presburger"
 
 using namespace mlir;
+using llvm::SmallDenseMap;
+using llvm::SmallDenseSet;
 
 std::unique_ptr<IntegerPolyhedron> IntegerPolyhedron::clone() const {
   return std::make_unique<IntegerPolyhedron>(*this);
@@ -307,6 +317,1778 @@ void IntegerPolyhedron::clearAndCopyFrom(const IntegerPolyhedron &other) {
   *this = other;
 }
 
+// Searches for a constraint with a non-zero coefficient at `colIdx` in
+// equality (isEq=true) or inequality (isEq=false) constraints.
+// Returns true and sets row found in search in `rowIdx`, false otherwise.
+bool IntegerPolyhedron::findConstraintWithNonZeroAt(unsigned colIdx, bool isEq,
+                                                    unsigned *rowIdx) const {
+  assert(colIdx < getNumCols() && "position out of bounds");
+  auto at = [&](unsigned rowIdx) -> int64_t {
+    return isEq ? atEq(rowIdx, colIdx) : atIneq(rowIdx, colIdx);
+  };
+  unsigned e = isEq ? getNumEqualities() : getNumInequalities();
+  for (*rowIdx = 0; *rowIdx < e; ++(*rowIdx)) {
+    if (at(*rowIdx) != 0) {
+      return true;
+    }
+  }
+  return false;
+}
+
+// Normalizes the coefficient values across all columns in `rowIdx` by their
+// GCD in equality or inequality constraints as specified by `isEq`.
+template <bool isEq>
+static void normalizeConstraintByGCD(IntegerPolyhedron *constraints,
+                                     unsigned rowIdx) {
+  auto at = [&](unsigned colIdx) -> int64_t {
+    return isEq ? constraints->atEq(rowIdx, colIdx)
+                : constraints->atIneq(rowIdx, colIdx);
+  };
+  uint64_t gcd = std::abs(at(0));
+  for (unsigned j = 1, e = constraints->getNumCols(); j < e; ++j) {
+    gcd = llvm::GreatestCommonDivisor64(gcd, std::abs(at(j)));
+  }
+  if (gcd > 0 && gcd != 1) {
+    for (unsigned j = 0, e = constraints->getNumCols(); j < e; ++j) {
+      int64_t v = at(j) / static_cast<int64_t>(gcd);
+      isEq ? constraints->atEq(rowIdx, j) = v
+           : constraints->atIneq(rowIdx, j) = v;
+    }
+  }
+}
+
+void IntegerPolyhedron::normalizeConstraintsByGCD() {
+  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
+    normalizeConstraintByGCD</*isEq=*/true>(this, i);
+  }
+  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
+    normalizeConstraintByGCD</*isEq=*/false>(this, i);
+  }
+}
+
+bool IntegerPolyhedron::hasInvalidConstraint() const {
+  assert(hasConsistentState());
+  auto check = [&](bool isEq) -> bool {
+    unsigned numCols = getNumCols();
+    unsigned numRows = isEq ? getNumEqualities() : getNumInequalities();
+    for (unsigned i = 0, e = numRows; i < e; ++i) {
+      unsigned j;
+      for (j = 0; j < numCols - 1; ++j) {
+        int64_t v = isEq ? atEq(i, j) : atIneq(i, j);
+        // Skip rows with non-zero variable coefficients.
+        if (v != 0)
+          break;
+      }
+      if (j < numCols - 1) {
+        continue;
+      }
+      // Check validity of constant term at 'numCols - 1' w.r.t 'isEq'.
+      // Example invalid constraints include: '1 == 0' or '-1 >= 0'
+      int64_t v = isEq ? atEq(i, numCols - 1) : atIneq(i, numCols - 1);
+      if ((isEq && v != 0) || (!isEq && v < 0)) {
+        return true;
+      }
+    }
+    return false;
+  };
+  if (check(/*isEq=*/true))
+    return true;
+  return check(/*isEq=*/false);
+}
+
+/// Eliminate identifier from constraint at `rowIdx` based on coefficient at
+/// pivotRow, pivotCol. Columns in range [elimColStart, pivotCol) will not be
+/// updated as they have already been eliminated.
+static void eliminateFromConstraint(IntegerPolyhedron *constraints,
+                                    unsigned rowIdx, unsigned pivotRow,
+                                    unsigned pivotCol, unsigned elimColStart,
+                                    bool isEq) {
+  // Skip if equality 'rowIdx' if same as 'pivotRow'.
+  if (isEq && rowIdx == pivotRow)
+    return;
+  auto at = [&](unsigned i, unsigned j) -> int64_t {
+    return isEq ? constraints->atEq(i, j) : constraints->atIneq(i, j);
+  };
+  int64_t leadCoeff = at(rowIdx, pivotCol);
+  // Skip if leading coefficient at 'rowIdx' is already zero.
+  if (leadCoeff == 0)
+    return;
+  int64_t pivotCoeff = constraints->atEq(pivotRow, pivotCol);
+  int64_t sign = (leadCoeff * pivotCoeff > 0) ? -1 : 1;
+  int64_t lcm = mlir::lcm(pivotCoeff, leadCoeff);
+  int64_t pivotMultiplier = sign * (lcm / std::abs(pivotCoeff));
+  int64_t rowMultiplier = lcm / std::abs(leadCoeff);
+
+  unsigned numCols = constraints->getNumCols();
+  for (unsigned j = 0; j < numCols; ++j) {
+    // Skip updating column 'j' if it was just eliminated.
+    if (j >= elimColStart && j < pivotCol)
+      continue;
+    int64_t v = pivotMultiplier * constraints->atEq(pivotRow, j) +
+                rowMultiplier * at(rowIdx, j);
+    isEq ? constraints->atEq(rowIdx, j) = v
+         : constraints->atIneq(rowIdx, j) = v;
+  }
+}
+
+/// Returns the position of the identifier that has the minimum <number of lower
+/// bounds> times <number of upper bounds> from the specified range of
+/// identifiers [start, end). It is often best to eliminate in the increasing
+/// order of these counts when doing Fourier-Motzkin elimination since FM adds
+/// that many new constraints.
+static unsigned getBestIdToEliminate(const IntegerPolyhedron &cst,
+                                     unsigned start, unsigned end) {
+  assert(start < cst.getNumIds() && end < cst.getNumIds() + 1);
+
+  auto getProductOfNumLowerUpperBounds = [&](unsigned pos) {
+    unsigned numLb = 0;
+    unsigned numUb = 0;
+    for (unsigned r = 0, e = cst.getNumInequalities(); r < e; r++) {
+      if (cst.atIneq(r, pos) > 0) {
+        ++numLb;
+      } else if (cst.atIneq(r, pos) < 0) {
+        ++numUb;
+      }
+    }
+    return numLb * numUb;
+  };
+
+  unsigned minLoc = start;
+  unsigned min = getProductOfNumLowerUpperBounds(start);
+  for (unsigned c = start + 1; c < end; c++) {
+    unsigned numLbUbProduct = getProductOfNumLowerUpperBounds(c);
+    if (numLbUbProduct < min) {
+      min = numLbUbProduct;
+      minLoc = c;
+    }
+  }
+  return minLoc;
+}
+
+// Checks for emptiness of the set by eliminating identifiers successively and
+// using the GCD test (on all equality constraints) and checking for trivially
+// invalid constraints. Returns 'true' if the constraint system is found to be
+// empty; false otherwise.
+bool IntegerPolyhedron::isEmpty() const {
+  if (isEmptyByGCDTest() || hasInvalidConstraint())
+    return true;
+
+  IntegerPolyhedron tmpCst(*this);
+
+  // First, eliminate as many local variables as possible using equalities.
+  tmpCst.removeRedundantLocalVars();
+  if (tmpCst.isEmptyByGCDTest() || tmpCst.hasInvalidConstraint())
+    return true;
+
+  // Eliminate as many identifiers as possible using Gaussian elimination.
+  unsigned currentPos = 0;
+  while (currentPos < tmpCst.getNumIds()) {
+    tmpCst.gaussianEliminateIds(currentPos, tmpCst.getNumIds());
+    ++currentPos;
+    // We check emptiness through trivial checks after eliminating each ID to
+    // detect emptiness early. Since the checks isEmptyByGCDTest() and
+    // hasInvalidConstraint() are linear time and single sweep on the constraint
+    // buffer, this appears reasonable - but can optimize in the future.
+    if (tmpCst.hasInvalidConstraint() || tmpCst.isEmptyByGCDTest())
+      return true;
+  }
+
+  // Eliminate the remaining using FM.
+  for (unsigned i = 0, e = tmpCst.getNumIds(); i < e; i++) {
+    tmpCst.fourierMotzkinEliminate(
+        getBestIdToEliminate(tmpCst, 0, tmpCst.getNumIds()));
+    // Check for a constraint explosion. This rarely happens in practice, but
+    // this check exists as a safeguard against improperly constructed
+    // constraint systems or artificially created arbitrarily complex systems
+    // that aren't the intended use case for IntegerPolyhedron. This is
+    // needed since FM has a worst case exponential complexity in theory.
+    if (tmpCst.getNumConstraints() >= kExplosionFactor * getNumIds()) {
+      LLVM_DEBUG(llvm::dbgs() << "FM constraint explosion detected\n");
+      return false;
+    }
+
+    // FM wouldn't have modified the equalities in any way. So no need to again
+    // run GCD test. Check for trivial invalid constraints.
+    if (tmpCst.hasInvalidConstraint())
+      return true;
+  }
+  return false;
+}
+
+// Runs the GCD test on all equality constraints. Returns 'true' if this test
+// fails on any equality. Returns 'false' otherwise.
+// This test can be used to disprove the existence of a solution. If it returns
+// true, no integer solution to the equality constraints can exist.
+//
+// GCD test definition:
+//
+// The equality constraint:
+//
+//  c_1*x_1 + c_2*x_2 + ... + c_n*x_n = c_0
+//
+// has an integer solution iff:
+//
+//  GCD of c_1, c_2, ..., c_n divides c_0.
+//
+bool IntegerPolyhedron::isEmptyByGCDTest() const {
+  assert(hasConsistentState());
+  unsigned numCols = getNumCols();
+  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
+    uint64_t gcd = std::abs(atEq(i, 0));
+    for (unsigned j = 1; j < numCols - 1; ++j) {
+      gcd = llvm::GreatestCommonDivisor64(gcd, std::abs(atEq(i, j)));
+    }
+    int64_t v = std::abs(atEq(i, numCols - 1));
+    if (gcd > 0 && (v % gcd != 0)) {
+      return true;
+    }
+  }
+  return false;
+}
+
+// Returns a matrix where each row is a vector along which the polytope is
+// bounded. The span of the returned vectors is guaranteed to contain all
+// such vectors. The returned vectors are NOT guaranteed to be linearly
+// independent. This function should not be called on empty sets.
+//
+// It is sufficient to check the perpendiculars of the constraints, as the set
+// of perpendiculars which are bounded must span all bounded directions.
+Matrix IntegerPolyhedron::getBoundedDirections() const {
+  // Note that it is necessary to add the equalities too (which the constructor
+  // does) even though we don't need to check if they are bounded; whether an
+  // inequality is bounded or not depends on what other constraints, including
+  // equalities, are present.
+  Simplex simplex(*this);
+
+  assert(!simplex.isEmpty() && "It is not meaningful to ask whether a "
+                               "direction is bounded in an empty set.");
+
+  SmallVector<unsigned, 8> boundedIneqs;
+  // The constructor adds the inequalities to the simplex first, so this
+  // processes all the inequalities.
+  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
+    if (simplex.isBoundedAlongConstraint(i))
+      boundedIneqs.push_back(i);
+  }
+
+  // The direction vector is given by the coefficients and does not include the
+  // constant term, so the matrix has one fewer column.
+  unsigned dirsNumCols = getNumCols() - 1;
+  Matrix dirs(boundedIneqs.size() + getNumEqualities(), dirsNumCols);
+
+  // Copy the bounded inequalities.
+  unsigned row = 0;
+  for (unsigned i : boundedIneqs) {
+    for (unsigned col = 0; col < dirsNumCols; ++col)
+      dirs(row, col) = atIneq(i, col);
+    ++row;
+  }
+
+  // Copy the equalities. All the equalities' perpendiculars are bounded.
+  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
+    for (unsigned col = 0; col < dirsNumCols; ++col)
+      dirs(row, col) = atEq(i, col);
+    ++row;
+  }
+
+  return dirs;
+}
+
+bool eqInvolvesSuffixDims(const IntegerPolyhedron &poly, unsigned eqIndex,
+                          unsigned numDims) {
+  for (unsigned e = poly.getNumIds(), j = e - numDims; j < e; ++j)
+    if (poly.atEq(eqIndex, j) != 0)
+      return true;
+  return false;
+}
+bool ineqInvolvesSuffixDims(const IntegerPolyhedron &poly, unsigned ineqIndex,
+                            unsigned numDims) {
+  for (unsigned e = poly.getNumIds(), j = e - numDims; j < e; ++j)
+    if (poly.atIneq(ineqIndex, j) != 0)
+      return true;
+  return false;
+}
+
+void removeConstraintsInvolvingSuffixDims(IntegerPolyhedron &poly,
+                                          unsigned unboundedDims) {
+  // We iterate backwards so that whether we remove constraint i - 1 or not, the
+  // next constraint to be tested is always i - 2.
+  for (unsigned i = poly.getNumEqualities(); i > 0; i--)
+    if (eqInvolvesSuffixDims(poly, i - 1, unboundedDims))
+      poly.removeEquality(i - 1);
+  for (unsigned i = poly.getNumInequalities(); i > 0; i--)
+    if (ineqInvolvesSuffixDims(poly, i - 1, unboundedDims))
+      poly.removeInequality(i - 1);
+}
+
+bool IntegerPolyhedron::isIntegerEmpty() const {
+  return !findIntegerSample().hasValue();
+}
+
+/// Let this set be S. If S is bounded then we directly call into the GBR
+/// sampling algorithm. Otherwise, there are some unbounded directions, i.e.,
+/// vectors v such that S extends to infinity along v or -v. In this case we
+/// use an algorithm described in the integer set library (isl) manual and used
+/// by the isl_set_sample function in that library. The algorithm is:
+///
+/// 1) Apply a unimodular transform T to S to obtain S*T, such that all
+/// dimensions in which S*T is bounded lie in the linear span of a prefix of the
+/// dimensions.
+///
+/// 2) Construct a set B by removing all constraints that involve
+/// the unbounded dimensions and then deleting the unbounded dimensions. Note
+/// that B is a Bounded set.
+///
+/// 3) Try to obtain a sample from B using the GBR sampling
+/// algorithm. If no sample is found, return that S is empty.
+///
+/// 4) Otherwise, substitute the obtained sample into S*T to obtain a set
+/// C. C is a full-dimensional Cone and always contains a sample.
+///
+/// 5) Obtain an integer sample from C.
+///
+/// 6) Return T*v, where v is the concatenation of the samples from B and C.
+///
+/// The following is a sketch of a proof that
+/// a) If the algorithm returns empty, then S is empty.
+/// b) If the algorithm returns a sample, it is a valid sample in S.
+///
+/// The algorithm returns empty only if B is empty, in which case S*T is
+/// certainly empty since B was obtained by removing constraints and then
+/// deleting unconstrained dimensions from S*T. Since T is unimodular, a vector
+/// v is in S*T iff T*v is in S. So in this case, since
+/// S*T is empty, S is empty too.
+///
+/// Otherwise, the algorithm substitutes the sample from B into S*T. All the
+/// constraints of S*T that did not involve unbounded dimensions are satisfied
+/// by this substitution. All dimensions in the linear span of the dimensions
+/// outside the prefix are unbounded in S*T (step 1). Substituting values for
+/// the bounded dimensions cannot make these dimensions bounded, and these are
+/// the only remaining dimensions in C, so C is unbounded along every vector (in
+/// the positive or negative direction, or both). C is hence a full-dimensional
+/// cone and therefore always contains an integer point.
+///
+/// Concatenating the samples from B and C gives a sample v in S*T, so the
+/// returned sample T*v is a sample in S.
+Optional<SmallVector<int64_t, 8>> IntegerPolyhedron::findIntegerSample() const {
+  // First, try the GCD test heuristic.
+  if (isEmptyByGCDTest())
+    return {};
+
+  Simplex simplex(*this);
+  if (simplex.isEmpty())
+    return {};
+
+  // For a bounded set, we directly call into the GBR sampling algorithm.
+  if (!simplex.isUnbounded())
+    return simplex.findIntegerSample();
+
+  // The set is unbounded. We cannot directly use the GBR algorithm.
+  //
+  // m is a matrix containing, in each row, a vector in which S is
+  // bounded, such that the linear span of all these dimensions contains all
+  // bounded dimensions in S.
+  Matrix m = getBoundedDirections();
+  // In column echelon form, each row of m occupies only the first rank(m)
+  // columns and has zeros on the other columns. The transform T that brings S
+  // to column echelon form is unimodular as well, so this is a suitable
+  // transform to use in step 1 of the algorithm.
+  std::pair<unsigned, LinearTransform> result =
+      LinearTransform::makeTransformToColumnEchelon(std::move(m));
+  const LinearTransform &transform = result.second;
+  // 1) Apply T to S to obtain S*T.
+  IntegerPolyhedron transformedSet = transform.applyTo(*this);
+
+  // 2) Remove the unbounded dimensions and constraints involving them to
+  // obtain a bounded set.
+  IntegerPolyhedron boundedSet(transformedSet);
+  unsigned numBoundedDims = result.first;
+  unsigned numUnboundedDims = getNumIds() - numBoundedDims;
+  removeConstraintsInvolvingSuffixDims(boundedSet, numUnboundedDims);
+  boundedSet.removeIdRange(numBoundedDims, boundedSet.getNumIds());
+
+  // 3) Try to obtain a sample from the bounded set.
+  Optional<SmallVector<int64_t, 8>> boundedSample =
+      Simplex(boundedSet).findIntegerSample();
+  if (!boundedSample)
+    return {};
+  assert(boundedSet.containsPoint(*boundedSample) &&
+         "Simplex returned an invalid sample!");
+
+  // 4) Substitute the values of the bounded dimensions into S*T to obtain a
+  // full-dimensional cone, which necessarily contains an integer sample.
+  transformedSet.setAndEliminate(0, *boundedSample);
+  IntegerPolyhedron &cone = transformedSet;
+
+  // 5) Obtain an integer sample from the cone.
+  //
+  // We shrink the cone such that for any rational point in the shrunken cone,
+  // rounding up each of the point's coordinates produces a point that still
+  // lies in the original cone.
+  //
+  // Rounding up a point x adds a number e_i in [0, 1) to each coordinate x_i.
+  // For each inequality sum_i a_i x_i + c >= 0 in the original cone, the
+  // shrunken cone will have the inequality tightened by some amount s, such
+  // that if x satisfies the shrunken cone's tightened inequality, then x + e
+  // satisfies the original inequality, i.e.,
+  //
+  // sum_i a_i x_i + c + s >= 0 implies sum_i a_i (x_i + e_i) + c >= 0
+  //
+  // for any e_i values in [0, 1). In fact, we will handle the slightly more
+  // general case where e_i can be in [0, 1]. For example, consider the
+  // inequality 2x_1 - 3x_2 - 7x_3 - 6 >= 0, and let x = (3, 0, 0). How low
+  // could the LHS go if we added a number in [0, 1] to each coordinate? The LHS
+  // is minimized when we add 1 to the x_i with negative coefficient a_i and
+  // keep the other x_i the same. In the example, we would get x = (3, 1, 1),
+  // changing the value of the LHS by -3 + -7 = -10.
+  //
+  // In general, the value of the LHS can change by at most the sum of the
+  // negative a_i, so we accomodate this by shifting the inequality by this
+  // amount for the shrunken cone.
+  for (unsigned i = 0, e = cone.getNumInequalities(); i < e; ++i) {
+    for (unsigned j = 0; j < cone.getNumIds(); ++j) {
+      int64_t coeff = cone.atIneq(i, j);
+      if (coeff < 0)
+        cone.atIneq(i, cone.getNumIds()) += coeff;
+    }
+  }
+
+  // Obtain an integer sample in the cone by rounding up a rational point from
+  // the shrunken cone. Shrinking the cone amounts to shifting its apex
+  // "inwards" without changing its "shape"; the shrunken cone is still a
+  // full-dimensional cone and is hence non-empty.
+  Simplex shrunkenConeSimplex(cone);
+  assert(!shrunkenConeSimplex.isEmpty() && "Shrunken cone cannot be empty!");
+  SmallVector<Fraction, 8> shrunkenConeSample =
+      shrunkenConeSimplex.getRationalSample();
+
+  SmallVector<int64_t, 8> coneSample(llvm::map_range(shrunkenConeSample, ceil));
+
+  // 6) Return transform * concat(boundedSample, coneSample).
+  SmallVector<int64_t, 8> &sample = boundedSample.getValue();
+  sample.append(coneSample.begin(), coneSample.end());
+  return transform.preMultiplyColumn(sample);
+}
+
+/// Helper to evaluate an affine expression at a point.
+/// The expression is a list of coefficients for the dimensions followed by the
+/// constant term.
+static int64_t valueAt(ArrayRef<int64_t> expr, ArrayRef<int64_t> point) {
+  assert(expr.size() == 1 + point.size() &&
+         "Dimensionalities of point and expression don't match!");
+  int64_t value = expr.back();
+  for (unsigned i = 0; i < point.size(); ++i)
+    value += expr[i] * point[i];
+  return value;
+}
+
+/// A point satisfies an equality iff the value of the equality at the
+/// expression is zero, and it satisfies an inequality iff the value of the
+/// inequality at that point is non-negative.
+bool IntegerPolyhedron::containsPoint(ArrayRef<int64_t> point) const {
+  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
+    if (valueAt(getEquality(i), point) != 0)
+      return false;
+  }
+  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
+    if (valueAt(getInequality(i), point) < 0)
+      return false;
+  }
+  return true;
+}
+
+void IntegerPolyhedron::getLocalReprs(
+    std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const {
+  std::vector<SmallVector<int64_t, 8>> dividends(getNumLocalIds());
+  SmallVector<unsigned, 4> denominators(getNumLocalIds());
+  getLocalReprs(dividends, denominators, repr);
+}
+
+void IntegerPolyhedron::getLocalReprs(
+    std::vector<SmallVector<int64_t, 8>> &dividends,
+    SmallVector<unsigned, 4> &denominators) const {
+  std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> repr(
+      getNumLocalIds());
+  getLocalReprs(dividends, denominators, repr);
+}
+
+void IntegerPolyhedron::getLocalReprs(
+    std::vector<SmallVector<int64_t, 8>> &dividends,
+    SmallVector<unsigned, 4> &denominators,
+    std::vector<llvm::Optional<std::pair<unsigned, unsigned>>> &repr) const {
+
+  repr.resize(getNumLocalIds());
+  dividends.resize(getNumLocalIds());
+  denominators.resize(getNumLocalIds());
+
+  SmallVector<bool, 8> foundRepr(getNumIds(), false);
+  for (unsigned i = 0, e = getNumDimAndSymbolIds(); i < e; ++i)
+    foundRepr[i] = true;
+
+  unsigned divOffset = getNumDimAndSymbolIds();
+  bool changed;
+  do {
+    // Each time changed is true, at end of this iteration, one or more local
+    // vars have been detected as floor divs.
+    changed = false;
+    for (unsigned i = 0, e = getNumLocalIds(); i < e; ++i) {
+      if (!foundRepr[i + divOffset]) {
+        if (auto res = presburger_utils::computeSingleVarRepr(
+                *this, foundRepr, divOffset + i, dividends[i],
+                denominators[i])) {
+          foundRepr[i + divOffset] = true;
+          repr[i] = res;
+          changed = true;
+        }
+      }
+    }
+  } while (changed);
+
+  // Set 0 denominator for identifiers for which no division representation
+  // could be found.
+  for (unsigned i = 0, e = repr.size(); i < e; ++i)
+    if (!repr[i].hasValue())
+      denominators[i] = 0;
+}
+
+/// Tightens inequalities given that we are dealing with integer spaces. This is
+/// analogous to the GCD test but applied to inequalities. The constant term can
+/// be reduced to the preceding multiple of the GCD of the coefficients, i.e.,
+///  64*i - 100 >= 0  =>  64*i - 128 >= 0 (since 'i' is an integer). This is a
+/// fast method - linear in the number of coefficients.
+// Example on how this affects practical cases: consider the scenario:
+// 64*i >= 100, j = 64*i; without a tightening, elimination of i would yield
+// j >= 100 instead of the tighter (exact) j >= 128.
+void IntegerPolyhedron::gcdTightenInequalities() {
+  unsigned numCols = getNumCols();
+  for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
+    uint64_t gcd = std::abs(atIneq(i, 0));
+    for (unsigned j = 1; j < numCols - 1; ++j) {
+      gcd = llvm::GreatestCommonDivisor64(gcd, std::abs(atIneq(i, j)));
+    }
+    if (gcd > 0 && gcd != 1) {
+      int64_t gcdI = static_cast<int64_t>(gcd);
+      // Tighten the constant term and normalize the constraint by the GCD.
+      atIneq(i, numCols - 1) = mlir::floorDiv(atIneq(i, numCols - 1), gcdI);
+      for (unsigned j = 0, e = numCols - 1; j < e; ++j)
+        atIneq(i, j) /= gcdI;
+    }
+  }
+}
+
+// Eliminates all identifier variables in column range [posStart, posLimit).
+// Returns the number of variables eliminated.
+unsigned IntegerPolyhedron::gaussianEliminateIds(unsigned posStart,
+                                                 unsigned posLimit) {
+  // Return if identifier positions to eliminate are out of range.
+  assert(posLimit <= numIds);
+  assert(hasConsistentState());
+
+  if (posStart >= posLimit)
+    return 0;
+
+  gcdTightenInequalities();
+
+  unsigned pivotCol = 0;
+  for (pivotCol = posStart; pivotCol < posLimit; ++pivotCol) {
+    // Find a row which has a non-zero coefficient in column 'j'.
+    unsigned pivotRow;
+    if (!findConstraintWithNonZeroAt(pivotCol, /*isEq=*/true, &pivotRow)) {
+      // No pivot row in equalities with non-zero at 'pivotCol'.
+      if (!findConstraintWithNonZeroAt(pivotCol, /*isEq=*/false, &pivotRow)) {
+        // If inequalities are also non-zero in 'pivotCol', it can be
+        // eliminated.
+        continue;
+      }
+      break;
+    }
+
+    // Eliminate identifier at 'pivotCol' from each equality row.
+    for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
+      eliminateFromConstraint(this, i, pivotRow, pivotCol, posStart,
+                              /*isEq=*/true);
+      normalizeConstraintByGCD</*isEq=*/true>(this, i);
+    }
+
+    // Eliminate identifier at 'pivotCol' from each inequality row.
+    for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
+      eliminateFromConstraint(this, i, pivotRow, pivotCol, posStart,
+                              /*isEq=*/false);
+      normalizeConstraintByGCD</*isEq=*/false>(this, i);
+    }
+    removeEquality(pivotRow);
+    gcdTightenInequalities();
+  }
+  // Update position limit based on number eliminated.
+  posLimit = pivotCol;
+  // Remove eliminated columns from all constraints.
+  removeIdRange(posStart, posLimit);
+  return posLimit - posStart;
+}
+
+// Fills an inequality row with the value 'val'.
+static inline void fillInequality(IntegerPolyhedron *cst, unsigned r,
+                                  int64_t val) {
+  for (unsigned c = 0, f = cst->getNumCols(); c < f; c++) {
+    cst->atIneq(r, c) = val;
+  }
+}
+
+// Negates an inequality.
+static inline void negateInequality(IntegerPolyhedron *cst, unsigned r) {
+  for (unsigned c = 0, f = cst->getNumCols(); c < f; c++) {
+    cst->atIneq(r, c) = -cst->atIneq(r, c);
+  }
+}
+
+// A more complex check to eliminate redundant inequalities. Uses FourierMotzkin
+// to check if a constraint is redundant.
+void IntegerPolyhedron::removeRedundantInequalities() {
+  SmallVector<bool, 32> redun(getNumInequalities(), false);
+  // To check if an inequality is redundant, we replace the inequality by its
+  // complement (for eg., i - 1 >= 0 by i <= 0), and check if the resulting
+  // system is empty. If it is, the inequality is redundant.
+  IntegerPolyhedron tmpCst(*this);
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
+    // Change the inequality to its complement.
+    negateInequality(&tmpCst, r);
+    tmpCst.atIneq(r, tmpCst.getNumCols() - 1)--;
+    if (tmpCst.isEmpty()) {
+      redun[r] = true;
+      // Zero fill the redundant inequality.
+      fillInequality(this, r, /*val=*/0);
+      fillInequality(&tmpCst, r, /*val=*/0);
+    } else {
+      // Reverse the change (to avoid recreating tmpCst each time).
+      tmpCst.atIneq(r, tmpCst.getNumCols() - 1)++;
+      negateInequality(&tmpCst, r);
+    }
+  }
+
+  // Scan to get rid of all rows marked redundant, in-place.
+  auto copyRow = [&](unsigned src, unsigned dest) {
+    if (src == dest)
+      return;
+    for (unsigned c = 0, e = getNumCols(); c < e; c++) {
+      atIneq(dest, c) = atIneq(src, c);
+    }
+  };
+  unsigned pos = 0;
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
+    if (!redun[r])
+      copyRow(r, pos++);
+  }
+  inequalities.resizeVertically(pos);
+}
+
+// A more complex check to eliminate redundant inequalities and equalities. Uses
+// Simplex to check if a constraint is redundant.
+void IntegerPolyhedron::removeRedundantConstraints() {
+  // First, we run gcdTightenInequalities. This allows us to catch some
+  // constraints which are not redundant when considering rational solutions
+  // but are redundant in terms of integer solutions.
+  gcdTightenInequalities();
+  Simplex simplex(*this);
+  simplex.detectRedundant();
+
+  auto copyInequality = [&](unsigned src, unsigned dest) {
+    if (src == dest)
+      return;
+    for (unsigned c = 0, e = getNumCols(); c < e; c++)
+      atIneq(dest, c) = atIneq(src, c);
+  };
+  unsigned pos = 0;
+  unsigned numIneqs = getNumInequalities();
+  // Scan to get rid of all inequalities marked redundant, in-place. In Simplex,
+  // the first constraints added are the inequalities.
+  for (unsigned r = 0; r < numIneqs; r++) {
+    if (!simplex.isMarkedRedundant(r))
+      copyInequality(r, pos++);
+  }
+  inequalities.resizeVertically(pos);
+
+  // Scan to get rid of all equalities marked redundant, in-place. In Simplex,
+  // after the inequalities, a pair of constraints for each equality is added.
+  // An equality is redundant if both the inequalities in its pair are
+  // redundant.
+  auto copyEquality = [&](unsigned src, unsigned dest) {
+    if (src == dest)
+      return;
+    for (unsigned c = 0, e = getNumCols(); c < e; c++)
+      atEq(dest, c) = atEq(src, c);
+  };
+  pos = 0;
+  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
+    if (!(simplex.isMarkedRedundant(numIneqs + 2 * r) &&
+          simplex.isMarkedRedundant(numIneqs + 2 * r + 1)))
+      copyEquality(r, pos++);
+  }
+  equalities.resizeVertically(pos);
+}
+
+/// Eliminate `pos2^th` local identifier, replacing its every instance with
+/// `pos1^th` local identifier. This function is intended to be used to remove
+/// redundancy when local variables at position `pos1` and `pos2` are restricted
+/// to have the same value.
+static void eliminateRedundantLocalId(IntegerPolyhedron &poly, unsigned pos1,
+                                      unsigned pos2) {
+
+  assert(pos1 < poly.getNumLocalIds() && "Invalid local id position");
+  assert(pos2 < poly.getNumLocalIds() && "Invalid local id position");
+
+  unsigned localOffset = poly.getNumDimAndSymbolIds();
+  pos1 += localOffset;
+  pos2 += localOffset;
+  for (unsigned i = 0, e = poly.getNumInequalities(); i < e; ++i)
+    poly.atIneq(i, pos1) += poly.atIneq(i, pos2);
+  for (unsigned i = 0, e = poly.getNumEqualities(); i < e; ++i)
+    poly.atEq(i, pos1) += poly.atEq(i, pos2);
+  poly.removeId(pos2);
+}
+
+/// Adds additional local ids to the sets such that they both have the union
+/// of the local ids in each set, without changing the set of points that
+/// lie in `this` and `other`.
+///
+/// To detect local ids that always take the same in both sets, each local id is
+/// represented as a floordiv with constant denominator in terms of other ids.
+/// After extracting these divisions, local ids with the same division
+/// representation are considered duplicate and are merged. It is possible that
+/// division representation for some local id cannot be obtained, and thus these
+/// local ids are not considered for detecting duplicates.
+void IntegerPolyhedron::mergeLocalIds(IntegerPolyhedron &other) {
+  assert(getNumDimIds() == other.getNumDimIds() &&
+         "Number of dimension ids should match");
+  assert(getNumSymbolIds() == other.getNumSymbolIds() &&
+         "Number of symbol ids should match");
+
+  IntegerPolyhedron &polyA = *this;
+  IntegerPolyhedron &polyB = other;
+
+  // Merge local ids of polyA and polyB without using division information,
+  // i.e. append local ids of `polyB` to `polyA` and insert local ids of `polyA`
+  // to `polyB` at start of its local ids.
+  unsigned initLocals = polyA.getNumLocalIds();
+  insertLocalId(polyA.getNumLocalIds(), polyB.getNumLocalIds());
+  polyB.insertLocalId(0, initLocals);
+
+  // Get division representations from each poly.
+  std::vector<SmallVector<int64_t, 8>> divsA, divsB;
+  SmallVector<unsigned, 4> denomsA, denomsB;
+  polyA.getLocalReprs(divsA, denomsA);
+  polyB.getLocalReprs(divsB, denomsB);
+
+  // Copy division information for polyB into `divsA` and `denomsA`, so that
+  // these have the combined division information of both polys. Since newly
+  // added local variables in polyA and polyB have no constraints, they will not
+  // have any division representation.
+  std::copy(divsB.begin() + initLocals, divsB.end(),
+            divsA.begin() + initLocals);
+  std::copy(denomsB.begin() + initLocals, denomsB.end(),
+            denomsA.begin() + initLocals);
+
+  // Find and merge duplicate divisions.
+  // TODO: Add division normalization to support divisions that 
diff er by
+  // a constant.
+  // TODO: Add division ordering such that a division representation for local
+  // identifier at position `i` only depends on local identifiers at position <
+  // `i`. This would make sure that all divisions depending on other local
+  // variables that can be merged, are merged.
+  unsigned localOffset = getIdKindOffset(IdKind::Local);
+  for (unsigned i = 0; i < divsA.size(); ++i) {
+    // Check if a division representation exists for the `i^th` local id.
+    if (denomsA[i] == 0)
+      continue;
+    // Check if a division exists which is a duplicate of the division at `i`.
+    for (unsigned j = i + 1; j < divsA.size(); ++j) {
+      // Check if a division representation exists for the `j^th` local id.
+      if (denomsA[j] == 0)
+        continue;
+      // Check if the denominators match.
+      if (denomsA[i] != denomsA[j])
+        continue;
+      // Check if the representations are equal.
+      if (divsA[i] != divsA[j])
+        continue;
+
+      // Merge divisions at position `j` into division at position `i`.
+      eliminateRedundantLocalId(polyA, i, j);
+      eliminateRedundantLocalId(polyB, i, j);
+      for (unsigned k = 0, g = divsA.size(); k < g; ++k) {
+        SmallVector<int64_t, 8> &div = divsA[k];
+        if (denomsA[k] != 0) {
+          div[localOffset + i] += div[localOffset + j];
+          div.erase(div.begin() + localOffset + j);
+        }
+      }
+
+      divsA.erase(divsA.begin() + j);
+      denomsA.erase(denomsA.begin() + j);
+      // Since `j` can never be zero, we do not need to worry about overflows.
+      --j;
+    }
+  }
+}
+
+/// Removes local variables using equalities. Each equality is checked if it
+/// can be reduced to the form: `e = affine-expr`, where `e` is a local
+/// variable and `affine-expr` is an affine expression not containing `e`.
+/// If an equality satisfies this form, the local variable is replaced in
+/// each constraint and then removed. The equality used to replace this local
+/// variable is also removed.
+void IntegerPolyhedron::removeRedundantLocalVars() {
+  // Normalize the equality constraints to reduce coefficients of local
+  // variables to 1 wherever possible.
+  for (unsigned i = 0, e = getNumEqualities(); i < e; ++i)
+    normalizeConstraintByGCD</*isEq=*/true>(this, i);
+
+  while (true) {
+    unsigned i, e, j, f;
+    for (i = 0, e = getNumEqualities(); i < e; ++i) {
+      // Find a local variable to eliminate using ith equality.
+      for (j = getNumDimAndSymbolIds(), f = getNumIds(); j < f; ++j)
+        if (std::abs(atEq(i, j)) == 1)
+          break;
+
+      // Local variable can be eliminated using ith equality.
+      if (j < f)
+        break;
+    }
+
+    // No equality can be used to eliminate a local variable.
+    if (i == e)
+      break;
+
+    // Use the ith equality to simplify other equalities. If any changes
+    // are made to an equality constraint, it is normalized by GCD.
+    for (unsigned k = 0, t = getNumEqualities(); k < t; ++k) {
+      if (atEq(k, j) != 0) {
+        eliminateFromConstraint(this, k, i, j, j, /*isEq=*/true);
+        normalizeConstraintByGCD</*isEq=*/true>(this, k);
+      }
+    }
+
+    // Use the ith equality to simplify inequalities.
+    for (unsigned k = 0, t = getNumInequalities(); k < t; ++k)
+      eliminateFromConstraint(this, k, i, j, j, /*isEq=*/false);
+
+    // Remove the ith equality and the found local variable.
+    removeId(j);
+    removeEquality(i);
+  }
+}
+
+void IntegerPolyhedron::convertDimToLocal(unsigned dimStart,
+                                          unsigned dimLimit) {
+  assert(dimLimit <= getNumDimIds() && "Invalid dim pos range");
+
+  if (dimStart >= dimLimit)
+    return;
+
+  // Append new local variables corresponding to the dimensions to be converted.
+  unsigned convertCount = dimLimit - dimStart;
+  unsigned newLocalIdStart = getNumIds();
+  appendLocalId(convertCount);
+
+  // Swap the new local variables with dimensions.
+  for (unsigned i = 0; i < convertCount; ++i)
+    swapId(i + dimStart, i + newLocalIdStart);
+
+  // Remove dimensions converted to local variables.
+  removeIdRange(dimStart, dimLimit);
+}
+
+void IntegerPolyhedron::addBound(BoundType type, unsigned pos, int64_t value) {
+  assert(pos < getNumCols());
+  if (type == BoundType::EQ) {
+    unsigned row = equalities.appendExtraRow();
+    equalities(row, pos) = 1;
+    equalities(row, getNumCols() - 1) = -value;
+  } else {
+    unsigned row = inequalities.appendExtraRow();
+    inequalities(row, pos) = type == BoundType::LB ? 1 : -1;
+    inequalities(row, getNumCols() - 1) =
+        type == BoundType::LB ? -value : value;
+  }
+}
+
+void IntegerPolyhedron::addBound(BoundType type, ArrayRef<int64_t> expr,
+                                 int64_t value) {
+  assert(type != BoundType::EQ && "EQ not implemented");
+  assert(expr.size() == getNumCols());
+  unsigned row = inequalities.appendExtraRow();
+  for (unsigned i = 0, e = expr.size(); i < e; ++i)
+    inequalities(row, i) = type == BoundType::LB ? expr[i] : -expr[i];
+  inequalities(inequalities.getNumRows() - 1, getNumCols() - 1) +=
+      type == BoundType::LB ? -value : value;
+}
+
+/// Adds a new local identifier as the floordiv of an affine function of other
+/// identifiers, the coefficients of which are provided in 'dividend' and with
+/// respect to a positive constant 'divisor'. Two constraints are added to the
+/// system to capture equivalence with the floordiv.
+///      q = expr floordiv c    <=>   c*q <= expr <= c*q + c - 1.
+void IntegerPolyhedron::addLocalFloorDiv(ArrayRef<int64_t> dividend,
+                                         int64_t divisor) {
+  assert(dividend.size() == getNumCols() && "incorrect dividend size");
+  assert(divisor > 0 && "positive divisor expected");
+
+  appendLocalId();
+
+  // Add two constraints for this new identifier 'q'.
+  SmallVector<int64_t, 8> bound(dividend.size() + 1);
+
+  // dividend - q * divisor >= 0
+  std::copy(dividend.begin(), dividend.begin() + dividend.size() - 1,
+            bound.begin());
+  bound.back() = dividend.back();
+  bound[getNumIds() - 1] = -divisor;
+  addInequality(bound);
+
+  // -dividend +qdivisor * q + divisor - 1 >= 0
+  std::transform(bound.begin(), bound.end(), bound.begin(),
+                 std::negate<int64_t>());
+  bound[bound.size() - 1] += divisor - 1;
+  addInequality(bound);
+}
+
+void IntegerPolyhedron::setDimSymbolSeparation(unsigned newSymbolCount) {
+  assert(newSymbolCount <= numDims + numSymbols &&
+         "invalid separation position");
+  numDims = numDims + numSymbols - newSymbolCount;
+  numSymbols = newSymbolCount;
+}
+
+/// Finds an equality that equates the specified identifier to a constant.
+/// Returns the position of the equality row. If 'symbolic' is set to true,
+/// symbols are also treated like a constant, i.e., an affine function of the
+/// symbols is also treated like a constant. Returns -1 if such an equality
+/// could not be found.
+static int findEqualityToConstant(const IntegerPolyhedron &cst, unsigned pos,
+                                  bool symbolic = false) {
+  assert(pos < cst.getNumIds() && "invalid position");
+  for (unsigned r = 0, e = cst.getNumEqualities(); r < e; r++) {
+    int64_t v = cst.atEq(r, pos);
+    if (v * v != 1)
+      continue;
+    unsigned c;
+    unsigned f = symbolic ? cst.getNumDimIds() : cst.getNumIds();
+    // This checks for zeros in all positions other than 'pos' in [0, f)
+    for (c = 0; c < f; c++) {
+      if (c == pos)
+        continue;
+      if (cst.atEq(r, c) != 0) {
+        // Dependent on another identifier.
+        break;
+      }
+    }
+    if (c == f)
+      // Equality is free of other identifiers.
+      return r;
+  }
+  return -1;
+}
+
+LogicalResult IntegerPolyhedron::constantFoldId(unsigned pos) {
+  assert(pos < getNumIds() && "invalid position");
+  int rowIdx;
+  if ((rowIdx = findEqualityToConstant(*this, pos)) == -1)
+    return failure();
+
+  // atEq(rowIdx, pos) is either -1 or 1.
+  assert(atEq(rowIdx, pos) * atEq(rowIdx, pos) == 1);
+  int64_t constVal = -atEq(rowIdx, getNumCols() - 1) / atEq(rowIdx, pos);
+  setAndEliminate(pos, constVal);
+  return success();
+}
+
+void IntegerPolyhedron::constantFoldIdRange(unsigned pos, unsigned num) {
+  for (unsigned s = pos, t = pos, e = pos + num; s < e; s++) {
+    if (failed(constantFoldId(t)))
+      t++;
+  }
+}
+
+/// Returns a non-negative constant bound on the extent (upper bound - lower
+/// bound) of the specified identifier if it is found to be a constant; returns
+/// None if it's not a constant. This methods treats symbolic identifiers
+/// specially, i.e., it looks for constant 
diff erences between affine
+/// expressions involving only the symbolic identifiers. See comments at
+/// function definition for example. 'lb', if provided, is set to the lower
+/// bound associated with the constant 
diff erence. Note that 'lb' is purely
+/// symbolic and thus will contain the coefficients of the symbolic identifiers
+/// and the constant coefficient.
+//  Egs: 0 <= i <= 15, return 16.
+//       s0 + 2 <= i <= s0 + 17, returns 16. (s0 has to be a symbol)
+//       s0 + s1 + 16 <= d0 <= s0 + s1 + 31, returns 16.
+//       s0 - 7 <= 8*j <= s0 returns 1 with lb = s0, lbDivisor = 8 (since lb =
+//       ceil(s0 - 7 / 8) = floor(s0 / 8)).
+Optional<int64_t> IntegerPolyhedron::getConstantBoundOnDimSize(
+    unsigned pos, SmallVectorImpl<int64_t> *lb, int64_t *boundFloorDivisor,
+    SmallVectorImpl<int64_t> *ub, unsigned *minLbPos,
+    unsigned *minUbPos) const {
+  assert(pos < getNumDimIds() && "Invalid identifier position");
+
+  // Find an equality for 'pos'^th identifier that equates it to some function
+  // of the symbolic identifiers (+ constant).
+  int eqPos = findEqualityToConstant(*this, pos, /*symbolic=*/true);
+  if (eqPos != -1) {
+    auto eq = getEquality(eqPos);
+    // If the equality involves a local var, punt for now.
+    // TODO: this can be handled in the future by using the explicit
+    // representation of the local vars.
+    if (!std::all_of(eq.begin() + getNumDimAndSymbolIds(), eq.end() - 1,
+                     [](int64_t coeff) { return coeff == 0; }))
+      return None;
+
+    // This identifier can only take a single value.
+    if (lb) {
+      // Set lb to that symbolic value.
+      lb->resize(getNumSymbolIds() + 1);
+      if (ub)
+        ub->resize(getNumSymbolIds() + 1);
+      for (unsigned c = 0, f = getNumSymbolIds() + 1; c < f; c++) {
+        int64_t v = atEq(eqPos, pos);
+        // atEq(eqRow, pos) is either -1 or 1.
+        assert(v * v == 1);
+        (*lb)[c] = v < 0 ? atEq(eqPos, getNumDimIds() + c) / -v
+                         : -atEq(eqPos, getNumDimIds() + c) / v;
+        // Since this is an equality, ub = lb.
+        if (ub)
+          (*ub)[c] = (*lb)[c];
+      }
+      assert(boundFloorDivisor &&
+             "both lb and divisor or none should be provided");
+      *boundFloorDivisor = 1;
+    }
+    if (minLbPos)
+      *minLbPos = eqPos;
+    if (minUbPos)
+      *minUbPos = eqPos;
+    return 1;
+  }
+
+  // Check if the identifier appears at all in any of the inequalities.
+  unsigned r, e;
+  for (r = 0, e = getNumInequalities(); r < e; r++) {
+    if (atIneq(r, pos) != 0)
+      break;
+  }
+  if (r == e)
+    // If it doesn't, there isn't a bound on it.
+    return None;
+
+  // Positions of constraints that are lower/upper bounds on the variable.
+  SmallVector<unsigned, 4> lbIndices, ubIndices;
+
+  // Gather all symbolic lower bounds and upper bounds of the variable, i.e.,
+  // the bounds can only involve symbolic (and local) identifiers. Since the
+  // canonical form c_1*x_1 + c_2*x_2 + ... + c_0 >= 0, a constraint is a lower
+  // bound for x_i if c_i >= 1, and an upper bound if c_i <= -1.
+  getLowerAndUpperBoundIndices(pos, &lbIndices, &ubIndices,
+                               /*eqIndices=*/nullptr, /*offset=*/0,
+                               /*num=*/getNumDimIds());
+
+  Optional<int64_t> minDiff = None;
+  unsigned minLbPosition = 0, minUbPosition = 0;
+  for (auto ubPos : ubIndices) {
+    for (auto lbPos : lbIndices) {
+      // Look for a lower bound and an upper bound that only 
diff er by a
+      // constant, i.e., pairs of the form  0 <= c_pos - f(c_i's) <= 
diff Const.
+      // For example, if ii is the pos^th variable, we are looking for
+      // constraints like ii >= i, ii <= ii + 50, 50 being the 
diff erence. The
+      // minimum among all such constant 
diff erences is kept since that's the
+      // constant bounding the extent of the pos^th variable.
+      unsigned j, e;
+      for (j = 0, e = getNumCols() - 1; j < e; j++)
+        if (atIneq(ubPos, j) != -atIneq(lbPos, j)) {
+          break;
+        }
+      if (j < getNumCols() - 1)
+        continue;
+      int64_t 
diff  = ceilDiv(atIneq(ubPos, getNumCols() - 1) +
+                                 atIneq(lbPos, getNumCols() - 1) + 1,
+                             atIneq(lbPos, pos));
+      // This bound is non-negative by definition.
+      
diff  = std::max<int64_t>(
diff , 0);
+      if (minDiff == None || 
diff  < minDiff) {
+        minDiff = 
diff ;
+        minLbPosition = lbPos;
+        minUbPosition = ubPos;
+      }
+    }
+  }
+  if (lb && minDiff.hasValue()) {
+    // Set lb to the symbolic lower bound.
+    lb->resize(getNumSymbolIds() + 1);
+    if (ub)
+      ub->resize(getNumSymbolIds() + 1);
+    // The lower bound is the ceildiv of the lb constraint over the coefficient
+    // of the variable at 'pos'. We express the ceildiv equivalently as a floor
+    // for uniformity. For eg., if the lower bound constraint was: 32*d0 - N +
+    // 31 >= 0, the lower bound for d0 is ceil(N - 31, 32), i.e., floor(N, 32).
+    *boundFloorDivisor = atIneq(minLbPosition, pos);
+    assert(*boundFloorDivisor == -atIneq(minUbPosition, pos));
+    for (unsigned c = 0, e = getNumSymbolIds() + 1; c < e; c++) {
+      (*lb)[c] = -atIneq(minLbPosition, getNumDimIds() + c);
+    }
+    if (ub) {
+      for (unsigned c = 0, e = getNumSymbolIds() + 1; c < e; c++)
+        (*ub)[c] = atIneq(minUbPosition, getNumDimIds() + c);
+    }
+    // The lower bound leads to a ceildiv while the upper bound is a floordiv
+    // whenever the coefficient at pos != 1. ceildiv (val / d) = floordiv (val +
+    // d - 1 / d); hence, the addition of 'atIneq(minLbPosition, pos) - 1' to
+    // the constant term for the lower bound.
+    (*lb)[getNumSymbolIds()] += atIneq(minLbPosition, pos) - 1;
+  }
+  if (minLbPos)
+    *minLbPos = minLbPosition;
+  if (minUbPos)
+    *minUbPos = minUbPosition;
+  return minDiff;
+}
+
+template <bool isLower>
+Optional<int64_t>
+IntegerPolyhedron::computeConstantLowerOrUpperBound(unsigned pos) {
+  assert(pos < getNumIds() && "invalid position");
+  // Project to 'pos'.
+  projectOut(0, pos);
+  projectOut(1, getNumIds() - 1);
+  // Check if there's an equality equating the '0'^th identifier to a constant.
+  int eqRowIdx = findEqualityToConstant(*this, 0, /*symbolic=*/false);
+  if (eqRowIdx != -1)
+    // atEq(rowIdx, 0) is either -1 or 1.
+    return -atEq(eqRowIdx, getNumCols() - 1) / atEq(eqRowIdx, 0);
+
+  // Check if the identifier appears at all in any of the inequalities.
+  unsigned r, e;
+  for (r = 0, e = getNumInequalities(); r < e; r++) {
+    if (atIneq(r, 0) != 0)
+      break;
+  }
+  if (r == e)
+    // If it doesn't, there isn't a bound on it.
+    return None;
+
+  Optional<int64_t> minOrMaxConst = None;
+
+  // Take the max across all const lower bounds (or min across all constant
+  // upper bounds).
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
+    if (isLower) {
+      if (atIneq(r, 0) <= 0)
+        // Not a lower bound.
+        continue;
+    } else if (atIneq(r, 0) >= 0) {
+      // Not an upper bound.
+      continue;
+    }
+    unsigned c, f;
+    for (c = 0, f = getNumCols() - 1; c < f; c++)
+      if (c != 0 && atIneq(r, c) != 0)
+        break;
+    if (c < getNumCols() - 1)
+      // Not a constant bound.
+      continue;
+
+    int64_t boundConst =
+        isLower ? mlir::ceilDiv(-atIneq(r, getNumCols() - 1), atIneq(r, 0))
+                : mlir::floorDiv(atIneq(r, getNumCols() - 1), -atIneq(r, 0));
+    if (isLower) {
+      if (minOrMaxConst == None || boundConst > minOrMaxConst)
+        minOrMaxConst = boundConst;
+    } else {
+      if (minOrMaxConst == None || boundConst < minOrMaxConst)
+        minOrMaxConst = boundConst;
+    }
+  }
+  return minOrMaxConst;
+}
+
+Optional<int64_t> IntegerPolyhedron::getConstantBound(BoundType type,
+                                                      unsigned pos) const {
+  assert(type != BoundType::EQ && "EQ not implemented");
+  IntegerPolyhedron tmpCst(*this);
+  if (type == BoundType::LB)
+    return tmpCst.computeConstantLowerOrUpperBound</*isLower=*/true>(pos);
+  return tmpCst.computeConstantLowerOrUpperBound</*isLower=*/false>(pos);
+}
+
+// A simple (naive and conservative) check for hyper-rectangularity.
+bool IntegerPolyhedron::isHyperRectangular(unsigned pos, unsigned num) const {
+  assert(pos < getNumCols() - 1);
+  // Check for two non-zero coefficients in the range [pos, pos + sum).
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
+    unsigned sum = 0;
+    for (unsigned c = pos; c < pos + num; c++) {
+      if (atIneq(r, c) != 0)
+        sum++;
+    }
+    if (sum > 1)
+      return false;
+  }
+  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
+    unsigned sum = 0;
+    for (unsigned c = pos; c < pos + num; c++) {
+      if (atEq(r, c) != 0)
+        sum++;
+    }
+    if (sum > 1)
+      return false;
+  }
+  return true;
+}
+
+/// Removes duplicate constraints, trivially true constraints, and constraints
+/// that can be detected as redundant as a result of 
diff ering only in their
+/// constant term part. A constraint of the form <non-negative constant> >= 0 is
+/// considered trivially true.
+//  Uses a DenseSet to hash and detect duplicates followed by a linear scan to
+//  remove duplicates in place.
+void IntegerPolyhedron::removeTrivialRedundancy() {
+  gcdTightenInequalities();
+  normalizeConstraintsByGCD();
+
+  // A map used to detect redundancy stemming from constraints that only 
diff er
+  // in their constant term. The value stored is <row position, const term>
+  // for a given row.
+  SmallDenseMap<ArrayRef<int64_t>, std::pair<unsigned, int64_t>>
+      rowsWithoutConstTerm;
+  // To unique rows.
+  SmallDenseSet<ArrayRef<int64_t>, 8> rowSet;
+
+  // Check if constraint is of the form <non-negative-constant> >= 0.
+  auto isTriviallyValid = [&](unsigned r) -> bool {
+    for (unsigned c = 0, e = getNumCols() - 1; c < e; c++) {
+      if (atIneq(r, c) != 0)
+        return false;
+    }
+    return atIneq(r, getNumCols() - 1) >= 0;
+  };
+
+  // Detect and mark redundant constraints.
+  SmallVector<bool, 256> redunIneq(getNumInequalities(), false);
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
+    int64_t *rowStart = &inequalities(r, 0);
+    auto row = ArrayRef<int64_t>(rowStart, getNumCols());
+    if (isTriviallyValid(r) || !rowSet.insert(row).second) {
+      redunIneq[r] = true;
+      continue;
+    }
+
+    // Among constraints that only 
diff er in the constant term part, mark
+    // everything other than the one with the smallest constant term redundant.
+    // (eg: among i - 16j - 5 >= 0, i - 16j - 1 >=0, i - 16j - 7 >= 0, the
+    // former two are redundant).
+    int64_t constTerm = atIneq(r, getNumCols() - 1);
+    auto rowWithoutConstTerm = ArrayRef<int64_t>(rowStart, getNumCols() - 1);
+    const auto &ret =
+        rowsWithoutConstTerm.insert({rowWithoutConstTerm, {r, constTerm}});
+    if (!ret.second) {
+      // Check if the other constraint has a higher constant term.
+      auto &val = ret.first->second;
+      if (val.second > constTerm) {
+        // The stored row is redundant. Mark it so, and update with this one.
+        redunIneq[val.first] = true;
+        val = {r, constTerm};
+      } else {
+        // The one stored makes this one redundant.
+        redunIneq[r] = true;
+      }
+    }
+  }
+
+  // Scan to get rid of all rows marked redundant, in-place.
+  unsigned pos = 0;
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++)
+    if (!redunIneq[r])
+      inequalities.copyRow(r, pos++);
+
+  inequalities.resizeVertically(pos);
+
+  // TODO: consider doing this for equalities as well, but probably not worth
+  // the savings.
+}
+
+static std::pair<unsigned, unsigned>
+getNewNumDimsSymbols(unsigned pos, const IntegerPolyhedron &cst) {
+  unsigned numDims = cst.getNumDimIds();
+  unsigned numSymbols = cst.getNumSymbolIds();
+  unsigned newNumDims, newNumSymbols;
+  if (pos < numDims) {
+    newNumDims = numDims - 1;
+    newNumSymbols = numSymbols;
+  } else if (pos < numDims + numSymbols) {
+    assert(numSymbols >= 1);
+    newNumDims = numDims;
+    newNumSymbols = numSymbols - 1;
+  } else {
+    newNumDims = numDims;
+    newNumSymbols = numSymbols;
+  }
+  return {newNumDims, newNumSymbols};
+}
+
+#undef DEBUG_TYPE
+#define DEBUG_TYPE "fm"
+
+/// Eliminates identifier at the specified position using Fourier-Motzkin
+/// variable elimination. This technique is exact for rational spaces but
+/// conservative (in "rare" cases) for integer spaces. The operation corresponds
+/// to a projection operation yielding the (convex) set of integer points
+/// contained in the rational shadow of the set. An emptiness test that relies
+/// on this method will guarantee emptiness, i.e., it disproves the existence of
+/// a solution if it says it's empty.
+/// If a non-null isResultIntegerExact is passed, it is set to true if the
+/// result is also integer exact. If it's set to false, the obtained solution
+/// *may* not be exact, i.e., it may contain integer points that do not have an
+/// integer pre-image in the original set.
+///
+/// Eg:
+/// j >= 0, j <= i + 1
+/// i >= 0, i <= N + 1
+/// Eliminating i yields,
+///   j >= 0, 0 <= N + 1, j - 1 <= N + 1
+///
+/// If darkShadow = true, this method computes the dark shadow on elimination;
+/// the dark shadow is a convex integer subset of the exact integer shadow. A
+/// non-empty dark shadow proves the existence of an integer solution. The
+/// elimination in such a case could however be an under-approximation, and thus
+/// should not be used for scanning sets or used by itself for dependence
+/// checking.
+///
+/// Eg: 2-d set, * represents grid points, 'o' represents a point in the set.
+///            ^
+///            |
+///            | * * * * o o
+///         i  | * * o o o o
+///            | o * * * * *
+///            --------------->
+///                 j ->
+///
+/// Eliminating i from this system (projecting on the j dimension):
+/// rational shadow / integer light shadow:  1 <= j <= 6
+/// dark shadow:                             3 <= j <= 6
+/// exact integer shadow:                    j = 1 \union  3 <= j <= 6
+/// holes/splinters:                         j = 2
+///
+/// darkShadow = false, isResultIntegerExact = nullptr are default values.
+// TODO: a slight modification to yield dark shadow version of FM (tightened),
+// which can prove the existence of a solution if there is one.
+void IntegerPolyhedron::fourierMotzkinEliminate(unsigned pos, bool darkShadow,
+                                                bool *isResultIntegerExact) {
+  LLVM_DEBUG(llvm::dbgs() << "FM input (eliminate pos " << pos << "):\n");
+  LLVM_DEBUG(dump());
+  assert(pos < getNumIds() && "invalid position");
+  assert(hasConsistentState());
+
+  // Check if this identifier can be eliminated through a substitution.
+  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
+    if (atEq(r, pos) != 0) {
+      // Use Gaussian elimination here (since we have an equality).
+      LogicalResult ret = gaussianEliminateId(pos);
+      (void)ret;
+      assert(succeeded(ret) && "Gaussian elimination guaranteed to succeed");
+      LLVM_DEBUG(llvm::dbgs() << "FM output (through Gaussian elimination):\n");
+      LLVM_DEBUG(dump());
+      return;
+    }
+  }
+
+  // A fast linear time tightening.
+  gcdTightenInequalities();
+
+  // Check if the identifier appears at all in any of the inequalities.
+  unsigned r, e;
+  for (r = 0, e = getNumInequalities(); r < e; r++) {
+    if (atIneq(r, pos) != 0)
+      break;
+  }
+  if (r == getNumInequalities()) {
+    // If it doesn't appear, just remove the column and return.
+    // TODO: refactor removeColumns to use it from here.
+    removeId(pos);
+    LLVM_DEBUG(llvm::dbgs() << "FM output:\n");
+    LLVM_DEBUG(dump());
+    return;
+  }
+
+  // Positions of constraints that are lower bounds on the variable.
+  SmallVector<unsigned, 4> lbIndices;
+  // Positions of constraints that are lower bounds on the variable.
+  SmallVector<unsigned, 4> ubIndices;
+  // Positions of constraints that do not involve the variable.
+  std::vector<unsigned> nbIndices;
+  nbIndices.reserve(getNumInequalities());
+
+  // Gather all lower bounds and upper bounds of the variable. Since the
+  // canonical form c_1*x_1 + c_2*x_2 + ... + c_0 >= 0, a constraint is a lower
+  // bound for x_i if c_i >= 1, and an upper bound if c_i <= -1.
+  for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
+    if (atIneq(r, pos) == 0) {
+      // Id does not appear in bound.
+      nbIndices.push_back(r);
+    } else if (atIneq(r, pos) >= 1) {
+      // Lower bound.
+      lbIndices.push_back(r);
+    } else {
+      // Upper bound.
+      ubIndices.push_back(r);
+    }
+  }
+
+  // Set the number of dimensions, symbols in the resulting system.
+  const auto &dimsSymbols = getNewNumDimsSymbols(pos, *this);
+  unsigned newNumDims = dimsSymbols.first;
+  unsigned newNumSymbols = dimsSymbols.second;
+
+  /// Create the new system which has one identifier less.
+  IntegerPolyhedron newPoly(
+      lbIndices.size() * ubIndices.size() + nbIndices.size(),
+      getNumEqualities(), getNumCols() - 1, newNumDims, newNumSymbols,
+      /*numLocals=*/getNumIds() - 1 - newNumDims - newNumSymbols);
+
+  // This will be used to check if the elimination was integer exact.
+  unsigned lcmProducts = 1;
+
+  // Let x be the variable we are eliminating.
+  // For each lower bound, lb <= c_l*x, and each upper bound c_u*x <= ub, (note
+  // that c_l, c_u >= 1) we have:
+  // lb*lcm(c_l, c_u)/c_l <= lcm(c_l, c_u)*x <= ub*lcm(c_l, c_u)/c_u
+  // We thus generate a constraint:
+  // lcm(c_l, c_u)/c_l*lb <= lcm(c_l, c_u)/c_u*ub.
+  // Note if c_l = c_u = 1, all integer points captured by the resulting
+  // constraint correspond to integer points in the original system (i.e., they
+  // have integer pre-images). Hence, if the lcm's are all 1, the elimination is
+  // integer exact.
+  for (auto ubPos : ubIndices) {
+    for (auto lbPos : lbIndices) {
+      SmallVector<int64_t, 4> ineq;
+      ineq.reserve(newPoly.getNumCols());
+      int64_t lbCoeff = atIneq(lbPos, pos);
+      // Note that in the comments above, ubCoeff is the negation of the
+      // coefficient in the canonical form as the view taken here is that of the
+      // term being moved to the other size of '>='.
+      int64_t ubCoeff = -atIneq(ubPos, pos);
+      // TODO: refactor this loop to avoid all branches inside.
+      for (unsigned l = 0, e = getNumCols(); l < e; l++) {
+        if (l == pos)
+          continue;
+        assert(lbCoeff >= 1 && ubCoeff >= 1 && "bounds wrongly identified");
+        int64_t lcm = mlir::lcm(lbCoeff, ubCoeff);
+        ineq.push_back(atIneq(ubPos, l) * (lcm / ubCoeff) +
+                       atIneq(lbPos, l) * (lcm / lbCoeff));
+        lcmProducts *= lcm;
+      }
+      if (darkShadow) {
+        // The dark shadow is a convex subset of the exact integer shadow. If
+        // there is a point here, it proves the existence of a solution.
+        ineq[ineq.size() - 1] += lbCoeff * ubCoeff - lbCoeff - ubCoeff + 1;
+      }
+      // TODO: we need to have a way to add inequalities in-place in
+      // IntegerPolyhedron instead of creating and copying over.
+      newPoly.addInequality(ineq);
+    }
+  }
+
+  LLVM_DEBUG(llvm::dbgs() << "FM isResultIntegerExact: " << (lcmProducts == 1)
+                          << "\n");
+  if (lcmProducts == 1 && isResultIntegerExact)
+    *isResultIntegerExact = true;
+
+  // Copy over the constraints not involving this variable.
+  for (auto nbPos : nbIndices) {
+    SmallVector<int64_t, 4> ineq;
+    ineq.reserve(getNumCols() - 1);
+    for (unsigned l = 0, e = getNumCols(); l < e; l++) {
+      if (l == pos)
+        continue;
+      ineq.push_back(atIneq(nbPos, l));
+    }
+    newPoly.addInequality(ineq);
+  }
+
+  assert(newPoly.getNumConstraints() ==
+         lbIndices.size() * ubIndices.size() + nbIndices.size());
+
+  // Copy over the equalities.
+  for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
+    SmallVector<int64_t, 4> eq;
+    eq.reserve(newPoly.getNumCols());
+    for (unsigned l = 0, e = getNumCols(); l < e; l++) {
+      if (l == pos)
+        continue;
+      eq.push_back(atEq(r, l));
+    }
+    newPoly.addEquality(eq);
+  }
+
+  // GCD tightening and normalization allows detection of more trivially
+  // redundant constraints.
+  newPoly.gcdTightenInequalities();
+  newPoly.normalizeConstraintsByGCD();
+  newPoly.removeTrivialRedundancy();
+  clearAndCopyFrom(newPoly);
+  LLVM_DEBUG(llvm::dbgs() << "FM output:\n");
+  LLVM_DEBUG(dump());
+}
+
+#undef DEBUG_TYPE
+#define DEBUG_TYPE "presburger"
+
+void IntegerPolyhedron::projectOut(unsigned pos, unsigned num) {
+  if (num == 0)
+    return;
+
+  // 'pos' can be at most getNumCols() - 2 if num > 0.
+  assert((getNumCols() < 2 || pos <= getNumCols() - 2) && "invalid position");
+  assert(pos + num < getNumCols() && "invalid range");
+
+  // Eliminate as many identifiers as possible using Gaussian elimination.
+  unsigned currentPos = pos;
+  unsigned numToEliminate = num;
+  unsigned numGaussianEliminated = 0;
+
+  while (currentPos < getNumIds()) {
+    unsigned curNumEliminated =
+        gaussianEliminateIds(currentPos, currentPos + numToEliminate);
+    ++currentPos;
+    numToEliminate -= curNumEliminated + 1;
+    numGaussianEliminated += curNumEliminated;
+  }
+
+  // Eliminate the remaining using Fourier-Motzkin.
+  for (unsigned i = 0; i < num - numGaussianEliminated; i++) {
+    unsigned numToEliminate = num - numGaussianEliminated - i;
+    fourierMotzkinEliminate(
+        getBestIdToEliminate(*this, pos, pos + numToEliminate));
+  }
+
+  // Fast/trivial simplifications.
+  gcdTightenInequalities();
+  // Normalize constraints after tightening since the latter impacts this, but
+  // not the other way round.
+  normalizeConstraintsByGCD();
+}
+
+namespace {
+
+enum BoundCmpResult { Greater, Less, Equal, Unknown };
+
+/// Compares two affine bounds whose coefficients are provided in 'first' and
+/// 'second'. The last coefficient is the constant term.
+static BoundCmpResult compareBounds(ArrayRef<int64_t> a, ArrayRef<int64_t> b) {
+  assert(a.size() == b.size());
+
+  // For the bounds to be comparable, their corresponding identifier
+  // coefficients should be equal; the constant terms are then compared to
+  // determine less/greater/equal.
+
+  if (!std::equal(a.begin(), a.end() - 1, b.begin()))
+    return Unknown;
+
+  if (a.back() == b.back())
+    return Equal;
+
+  return a.back() < b.back() ? Less : Greater;
+}
+} // namespace
+
+// Returns constraints that are common to both A & B.
+static void getCommonConstraints(const IntegerPolyhedron &a,
+                                 const IntegerPolyhedron &b,
+                                 IntegerPolyhedron &c) {
+  c.reset(a.getNumDimIds(), a.getNumSymbolIds(), a.getNumLocalIds());
+  // a naive O(n^2) check should be enough here given the input sizes.
+  for (unsigned r = 0, e = a.getNumInequalities(); r < e; ++r) {
+    for (unsigned s = 0, f = b.getNumInequalities(); s < f; ++s) {
+      if (a.getInequality(r) == b.getInequality(s)) {
+        c.addInequality(a.getInequality(r));
+        break;
+      }
+    }
+  }
+  for (unsigned r = 0, e = a.getNumEqualities(); r < e; ++r) {
+    for (unsigned s = 0, f = b.getNumEqualities(); s < f; ++s) {
+      if (a.getEquality(r) == b.getEquality(s)) {
+        c.addEquality(a.getEquality(r));
+        break;
+      }
+    }
+  }
+}
+
+// Computes the bounding box with respect to 'other' by finding the min of the
+// lower bounds and the max of the upper bounds along each of the dimensions.
+LogicalResult
+IntegerPolyhedron::unionBoundingBox(const IntegerPolyhedron &otherCst) {
+  assert(otherCst.getNumDimIds() == numDims && "dims mismatch");
+  assert(otherCst.getNumLocalIds() == 0 && "local ids not supported here");
+  assert(getNumLocalIds() == 0 && "local ids not supported yet here");
+
+  // Get the constraints common to both systems; these will be added as is to
+  // the union.
+  IntegerPolyhedron commonCst;
+  getCommonConstraints(*this, otherCst, commonCst);
+
+  std::vector<SmallVector<int64_t, 8>> boundingLbs;
+  std::vector<SmallVector<int64_t, 8>> boundingUbs;
+  boundingLbs.reserve(2 * getNumDimIds());
+  boundingUbs.reserve(2 * getNumDimIds());
+
+  // To hold lower and upper bounds for each dimension.
+  SmallVector<int64_t, 4> lb, otherLb, ub, otherUb;
+  // To compute min of lower bounds and max of upper bounds for each dimension.
+  SmallVector<int64_t, 4> minLb(getNumSymbolIds() + 1);
+  SmallVector<int64_t, 4> maxUb(getNumSymbolIds() + 1);
+  // To compute final new lower and upper bounds for the union.
+  SmallVector<int64_t, 8> newLb(getNumCols()), newUb(getNumCols());
+
+  int64_t lbFloorDivisor, otherLbFloorDivisor;
+  for (unsigned d = 0, e = getNumDimIds(); d < e; ++d) {
+    auto extent = getConstantBoundOnDimSize(d, &lb, &lbFloorDivisor, &ub);
+    if (!extent.hasValue())
+      // TODO: symbolic extents when necessary.
+      // TODO: handle union if a dimension is unbounded.
+      return failure();
+
+    auto otherExtent = otherCst.getConstantBoundOnDimSize(
+        d, &otherLb, &otherLbFloorDivisor, &otherUb);
+    if (!otherExtent.hasValue() || lbFloorDivisor != otherLbFloorDivisor)
+      // TODO: symbolic extents when necessary.
+      return failure();
+
+    assert(lbFloorDivisor > 0 && "divisor always expected to be positive");
+
+    auto res = compareBounds(lb, otherLb);
+    // Identify min.
+    if (res == BoundCmpResult::Less || res == BoundCmpResult::Equal) {
+      minLb = lb;
+      // Since the divisor is for a floordiv, we need to convert to ceildiv,
+      // i.e., i >= expr floordiv div <=> i >= (expr - div + 1) ceildiv div <=>
+      // div * i >= expr - div + 1.
+      minLb.back() -= lbFloorDivisor - 1;
+    } else if (res == BoundCmpResult::Greater) {
+      minLb = otherLb;
+      minLb.back() -= otherLbFloorDivisor - 1;
+    } else {
+      // Uncomparable - check for constant lower/upper bounds.
+      auto constLb = getConstantBound(BoundType::LB, d);
+      auto constOtherLb = otherCst.getConstantBound(BoundType::LB, d);
+      if (!constLb.hasValue() || !constOtherLb.hasValue())
+        return failure();
+      std::fill(minLb.begin(), minLb.end(), 0);
+      minLb.back() = std::min(constLb.getValue(), constOtherLb.getValue());
+    }
+
+    // Do the same for ub's but max of upper bounds. Identify max.
+    auto uRes = compareBounds(ub, otherUb);
+    if (uRes == BoundCmpResult::Greater || uRes == BoundCmpResult::Equal) {
+      maxUb = ub;
+    } else if (uRes == BoundCmpResult::Less) {
+      maxUb = otherUb;
+    } else {
+      // Uncomparable - check for constant lower/upper bounds.
+      auto constUb = getConstantBound(BoundType::UB, d);
+      auto constOtherUb = otherCst.getConstantBound(BoundType::UB, d);
+      if (!constUb.hasValue() || !constOtherUb.hasValue())
+        return failure();
+      std::fill(maxUb.begin(), maxUb.end(), 0);
+      maxUb.back() = std::max(constUb.getValue(), constOtherUb.getValue());
+    }
+
+    std::fill(newLb.begin(), newLb.end(), 0);
+    std::fill(newUb.begin(), newUb.end(), 0);
+
+    // The divisor for lb, ub, otherLb, otherUb at this point is lbDivisor,
+    // and so it's the divisor for newLb and newUb as well.
+    newLb[d] = lbFloorDivisor;
+    newUb[d] = -lbFloorDivisor;
+    // Copy over the symbolic part + constant term.
+    std::copy(minLb.begin(), minLb.end(), newLb.begin() + getNumDimIds());
+    std::transform(newLb.begin() + getNumDimIds(), newLb.end(),
+                   newLb.begin() + getNumDimIds(), std::negate<int64_t>());
+    std::copy(maxUb.begin(), maxUb.end(), newUb.begin() + getNumDimIds());
+
+    boundingLbs.push_back(newLb);
+    boundingUbs.push_back(newUb);
+  }
+
+  // Clear all constraints and add the lower/upper bounds for the bounding box.
+  clearConstraints();
+  for (unsigned d = 0, e = getNumDimIds(); d < e; ++d) {
+    addInequality(boundingLbs[d]);
+    addInequality(boundingUbs[d]);
+  }
+
+  // Add the constraints that were common to both systems.
+  append(commonCst);
+  removeTrivialRedundancy();
+
+  // TODO: copy over pure symbolic constraints from this and 'other' over to the
+  // union (since the above are just the union along dimensions); we shouldn't
+  // be discarding any other constraints on the symbols.
+
+  return success();
+}
+
+bool IntegerPolyhedron::isColZero(unsigned pos) const {
+  unsigned rowPos;
+  return !findConstraintWithNonZeroAt(pos, /*isEq=*/false, &rowPos) &&
+         !findConstraintWithNonZeroAt(pos, /*isEq=*/true, &rowPos);
+}
+
+/// Find positions of inequalities and equalities that do not have a coefficient
+/// for [pos, pos + num) identifiers.
+static void getIndependentConstraints(const IntegerPolyhedron &cst,
+                                      unsigned pos, unsigned num,
+                                      SmallVectorImpl<unsigned> &nbIneqIndices,
+                                      SmallVectorImpl<unsigned> &nbEqIndices) {
+  assert(pos < cst.getNumIds() && "invalid start position");
+  assert(pos + num <= cst.getNumIds() && "invalid limit");
+
+  for (unsigned r = 0, e = cst.getNumInequalities(); r < e; r++) {
+    // The bounds are to be independent of [offset, offset + num) columns.
+    unsigned c;
+    for (c = pos; c < pos + num; ++c) {
+      if (cst.atIneq(r, c) != 0)
+        break;
+    }
+    if (c == pos + num)
+      nbIneqIndices.push_back(r);
+  }
+
+  for (unsigned r = 0, e = cst.getNumEqualities(); r < e; r++) {
+    // The bounds are to be independent of [offset, offset + num) columns.
+    unsigned c;
+    for (c = pos; c < pos + num; ++c) {
+      if (cst.atEq(r, c) != 0)
+        break;
+    }
+    if (c == pos + num)
+      nbEqIndices.push_back(r);
+  }
+}
+
+void IntegerPolyhedron::removeIndependentConstraints(unsigned pos,
+                                                     unsigned num) {
+  assert(pos + num <= getNumIds() && "invalid range");
+
+  // Remove constraints that are independent of these identifiers.
+  SmallVector<unsigned, 4> nbIneqIndices, nbEqIndices;
+  getIndependentConstraints(*this, /*pos=*/0, num, nbIneqIndices, nbEqIndices);
+
+  // Iterate in reverse so that indices don't have to be updated.
+  // TODO: This method can be made more efficient (because removal of each
+  // inequality leads to much shifting/copying in the underlying buffer).
+  for (auto nbIndex : llvm::reverse(nbIneqIndices))
+    removeInequality(nbIndex);
+  for (auto nbIndex : llvm::reverse(nbEqIndices))
+    removeEquality(nbIndex);
+}
+
 void IntegerPolyhedron::printSpace(raw_ostream &os) const {
   os << "\nConstraints (" << getNumDimIds() << " dims, " << getNumSymbolIds()
      << " symbols, " << getNumLocalIds() << " locals), (" << getNumConstraints()

diff  --git a/mlir/unittests/Analysis/AffineStructuresTest.cpp b/mlir/unittests/Analysis/AffineStructuresTest.cpp
deleted file mode 100644
index 0d1d642cbf82..000000000000
--- a/mlir/unittests/Analysis/AffineStructuresTest.cpp
+++ /dev/null
@@ -1,902 +0,0 @@
-//===- AffineStructuresTest.cpp - Tests for AffineStructures ----*- C++ -*-===//
-//
-// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
-// See https://llvm.org/LICENSE.txt for license information.
-// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
-//
-//===----------------------------------------------------------------------===//
-
-#include "mlir/Analysis/AffineStructures.h"
-#include "./AffineStructuresParser.h"
-#include "mlir/IR/IntegerSet.h"
-#include "mlir/IR/MLIRContext.h"
-
-#include <gmock/gmock.h>
-#include <gtest/gtest.h>
-
-#include <numeric>
-
-namespace mlir {
-
-using testing::ElementsAre;
-
-enum class TestFunction { Sample, Empty };
-
-/// If fn is TestFunction::Sample (default):
-/// If hasSample is true, check that findIntegerSample returns a valid sample
-/// for the FlatAffineConstraints fac.
-/// If hasSample is false, check that findIntegerSample returns None.
-///
-/// If fn is TestFunction::Empty, check that isIntegerEmpty returns the
-/// opposite of hasSample.
-static void checkSample(bool hasSample, const FlatAffineConstraints &fac,
-                        TestFunction fn = TestFunction::Sample) {
-  Optional<SmallVector<int64_t, 8>> maybeSample;
-  switch (fn) {
-  case TestFunction::Sample:
-    maybeSample = fac.findIntegerSample();
-    if (!hasSample) {
-      EXPECT_FALSE(maybeSample.hasValue());
-      if (maybeSample.hasValue()) {
-        for (auto x : *maybeSample)
-          llvm::errs() << x << ' ';
-        llvm::errs() << '\n';
-      }
-    } else {
-      ASSERT_TRUE(maybeSample.hasValue());
-      EXPECT_TRUE(fac.containsPoint(*maybeSample));
-    }
-    break;
-  case TestFunction::Empty:
-    EXPECT_EQ(!hasSample, fac.isIntegerEmpty());
-    break;
-  }
-}
-
-/// Construct a FlatAffineConstraints from a set of inequality and
-/// equality constraints.
-static FlatAffineConstraints
-makeFACFromConstraints(unsigned ids, ArrayRef<SmallVector<int64_t, 4>> ineqs,
-                       ArrayRef<SmallVector<int64_t, 4>> eqs,
-                       unsigned syms = 0) {
-  FlatAffineConstraints fac(ineqs.size(), eqs.size(), ids + 1, ids - syms, syms,
-                            /*numLocals=*/0);
-  for (const auto &eq : eqs)
-    fac.addEquality(eq);
-  for (const auto &ineq : ineqs)
-    fac.addInequality(ineq);
-  return fac;
-}
-
-/// Check sampling for all the permutations of the dimensions for the given
-/// constraint set. Since the GBR algorithm progresses dimension-wise, 
diff erent
-/// orderings may cause the algorithm to proceed 
diff erently. At least some of
-///.these permutations should make it past the heuristics and test the
-/// implementation of the GBR algorithm itself.
-/// Use TestFunction fn to test.
-static void checkPermutationsSample(bool hasSample, unsigned nDim,
-                                    ArrayRef<SmallVector<int64_t, 4>> ineqs,
-                                    ArrayRef<SmallVector<int64_t, 4>> eqs,
-                                    TestFunction fn = TestFunction::Sample) {
-  SmallVector<unsigned, 4> perm(nDim);
-  std::iota(perm.begin(), perm.end(), 0);
-  auto permute = [&perm](ArrayRef<int64_t> coeffs) {
-    SmallVector<int64_t, 4> permuted;
-    for (unsigned id : perm)
-      permuted.push_back(coeffs[id]);
-    permuted.push_back(coeffs.back());
-    return permuted;
-  };
-  do {
-    SmallVector<SmallVector<int64_t, 4>, 4> permutedIneqs, permutedEqs;
-    for (const auto &ineq : ineqs)
-      permutedIneqs.push_back(permute(ineq));
-    for (const auto &eq : eqs)
-      permutedEqs.push_back(permute(eq));
-
-    checkSample(hasSample,
-                makeFACFromConstraints(nDim, permutedIneqs, permutedEqs), fn);
-  } while (std::next_permutation(perm.begin(), perm.end()));
-}
-
-/// Parses a FlatAffineConstraints from a StringRef. It is expected that the
-/// string represents a valid IntegerSet, otherwise it will violate a gtest
-/// assertion.
-static FlatAffineConstraints parseFAC(StringRef str, MLIRContext *context) {
-  FailureOr<FlatAffineConstraints> fac = parseIntegerSetToFAC(str, context);
-
-  EXPECT_TRUE(succeeded(fac));
-
-  return *fac;
-}
-
-TEST(FlatAffineConstraintsTest, FindSampleTest) {
-  // Bounded sets with only inequalities.
-
-  MLIRContext context;
-
-  // 0 <= 7x <= 5
-  checkSample(true, parseFAC("(x) : (7 * x >= 0, -7 * x + 5 >= 0)", &context));
-
-  // 1 <= 5x and 5x <= 4 (no solution).
-  checkSample(false,
-              parseFAC("(x) : (5 * x - 1 >= 0, -5 * x + 4 >= 0)", &context));
-
-  // 1 <= 5x and 5x <= 9 (solution: x = 1).
-  checkSample(true,
-              parseFAC("(x) : (5 * x - 1 >= 0, -5 * x + 9 >= 0)", &context));
-
-  // Bounded sets with equalities.
-  // x >= 8 and 40 >= y and x = y.
-  checkSample(true, parseFAC("(x,y) : (x - 8 >= 0, -y + 40 >= 0, x - y == 0)",
-                             &context));
-
-  // x <= 10 and y <= 10 and 10 <= z and x + 2y = 3z.
-  // solution: x = y = z = 10.
-  checkSample(true, parseFAC("(x,y,z) : (-x + 10 >= 0, -y + 10 >= 0, "
-                             "z - 10 >= 0, x + 2 * y - 3 * z == 0)",
-                             &context));
-
-  // x <= 10 and y <= 10 and 11 <= z and x + 2y = 3z.
-  // This implies x + 2y >= 33 and x + 2y <= 30, which has no solution.
-  checkSample(false, parseFAC("(x,y,z) : (-x + 10 >= 0, -y + 10 >= 0, "
-                              "z - 11 >= 0, x + 2 * y - 3 * z == 0)",
-                              &context));
-
-  // 0 <= r and r <= 3 and 4q + r = 7.
-  // Solution: q = 1, r = 3.
-  checkSample(
-      true,
-      parseFAC("(q,r) : (r >= 0, -r + 3 >= 0, 4 * q + r - 7 == 0)", &context));
-
-  // 4q + r = 7 and r = 0.
-  // Solution: q = 1, r = 3.
-  checkSample(false,
-              parseFAC("(q,r) : (4 * q + r - 7 == 0, r == 0)", &context));
-
-  // The next two sets are large sets that should take a long time to sample
-  // with a naive branch and bound algorithm but can be sampled efficiently with
-  // the GBR algorithm.
-  //
-  // This is a triangle with vertices at (1/3, 0), (2/3, 0) and (10000, 10000).
-  checkSample(true, parseFAC("(x,y) : (y >= 0, "
-                             "300000 * x - 299999 * y - 100000 >= 0, "
-                             "-300000 * x + 299998 * y + 200000 >= 0)",
-                             &context));
-
-  // This is a tetrahedron with vertices at
-  // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 10000), and (10000, 10000, 10000).
-  // The first three points form a triangular base on the xz plane with the
-  // apex at the fourth point, which is the only integer point.
-  checkPermutationsSample(
-      true, 3,
-      {
-          {0, 1, 0, 0},  // y >= 0
-          {0, -1, 1, 0}, // z >= y
-          {300000, -299998, -1,
-           -100000},                    // -300000x + 299998y + 100000 + z <= 0.
-          {-150000, 149999, 0, 100000}, // -150000x + 149999y + 100000 >= 0.
-      },
-      {});
-
-  // Same thing with some spurious extra dimensions equated to constants.
-  checkSample(
-      true,
-      parseFAC("(a,b,c,d,e) : (b + d - e >= 0, -b + c - d + e >= 0, "
-               "300000 * a - 299998 * b - c - 9 * d + 21 * e - 112000 >= 0, "
-               "-150000 * a + 149999 * b - 15 * d + 47 * e + 68000 >= 0, "
-               "d - e == 0, d + e - 2000 == 0)",
-               &context));
-
-  // This is a tetrahedron with vertices at
-  // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 100), (100, 100 - 1/3, 100).
-  checkPermutationsSample(false, 3,
-                          {
-                              {0, 1, 0, 0},
-                              {0, -300, 299, 0},
-                              {300 * 299, -89400, -299, -100 * 299},
-                              {-897, 894, 0, 598},
-                          },
-                          {});
-
-  // Two tests involving equalities that are integer empty but not rational
-  // empty.
-
-  // This is a line segment from (0, 1/3) to (100, 100 + 1/3).
-  checkSample(
-      false, parseFAC("(x,y) : (x >= 0, -x + 100 >= 0, 3 * x - 3 * y + 1 == 0)",
-                      &context));
-
-  // A thin parallelogram. 0 <= x <= 100 and x + 1/3 <= y <= x + 2/3.
-  checkSample(false,
-              parseFAC("(x,y) : (x >= 0, -x + 100 >= 0, "
-                       "3 * x - 3 * y + 2 >= 0, -3 * x + 3 * y - 1 >= 0)",
-                       &context));
-
-  checkSample(true, parseFAC("(x,y) : (2 * x >= 0, -2 * x + 99 >= 0, "
-                             "2 * y >= 0, -2 * y + 99 >= 0)",
-                             &context));
-
-  // 2D cone with apex at (10000, 10000) and
-  // edges passing through (1/3, 0) and (2/3, 0).
-  checkSample(true, parseFAC("(x,y) : (300000 * x - 299999 * y - 100000 >= 0, "
-                             "-300000 * x + 299998 * y + 200000 >= 0)",
-                             &context));
-
-  // Cartesian product of a tetrahedron and a 2D cone.
-  // The tetrahedron has vertices at
-  // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 10000), and (10000, 10000, 10000).
-  // The first three points form a triangular base on the xz plane with the
-  // apex at the fourth point, which is the only integer point.
-  // The cone has apex at (10000, 10000) and
-  // edges passing through (1/3, 0) and (2/3, 0).
-  checkPermutationsSample(
-      true /* not empty */, 5,
-      {
-          // Tetrahedron contraints:
-          {0, 1, 0, 0, 0, 0},  // y >= 0
-          {0, -1, 1, 0, 0, 0}, // z >= y
-                               // -300000x + 299998y + 100000 + z <= 0.
-          {300000, -299998, -1, 0, 0, -100000},
-          // -150000x + 149999y + 100000 >= 0.
-          {-150000, 149999, 0, 0, 0, 100000},
-
-          // Triangle constraints:
-          // 300000p - 299999q >= 100000
-          {0, 0, 0, 300000, -299999, -100000},
-          // -300000p + 299998q + 200000 >= 0
-          {0, 0, 0, -300000, 299998, 200000},
-      },
-      {});
-
-  // Cartesian product of same tetrahedron as above and {(p, q) : 1/3 <= p <=
-  // 2/3}. Since the second set is empty, the whole set is too.
-  checkPermutationsSample(
-      false /* empty */, 5,
-      {
-          // Tetrahedron contraints:
-          {0, 1, 0, 0, 0, 0},  // y >= 0
-          {0, -1, 1, 0, 0, 0}, // z >= y
-                               // -300000x + 299998y + 100000 + z <= 0.
-          {300000, -299998, -1, 0, 0, -100000},
-          // -150000x + 149999y + 100000 >= 0.
-          {-150000, 149999, 0, 0, 0, 100000},
-
-          // Second set constraints:
-          // 3p >= 1
-          {0, 0, 0, 3, 0, -1},
-          // 3p <= 2
-          {0, 0, 0, -3, 0, 2},
-      },
-      {});
-
-  // Cartesian product of same tetrahedron as above and
-  // {(p, q, r) : 1 <= p <= 2 and p = 3q + 3r}.
-  // Since the second set is empty, the whole set is too.
-  checkPermutationsSample(
-      false /* empty */, 5,
-      {
-          // Tetrahedron contraints:
-          {0, 1, 0, 0, 0, 0, 0},  // y >= 0
-          {0, -1, 1, 0, 0, 0, 0}, // z >= y
-                                  // -300000x + 299998y + 100000 + z <= 0.
-          {300000, -299998, -1, 0, 0, 0, -100000},
-          // -150000x + 149999y + 100000 >= 0.
-          {-150000, 149999, 0, 0, 0, 0, 100000},
-
-          // Second set constraints:
-          // p >= 1
-          {0, 0, 0, 1, 0, 0, -1},
-          // p <= 2
-          {0, 0, 0, -1, 0, 0, 2},
-      },
-      {
-          {0, 0, 0, 1, -3, -3, 0}, // p = 3q + 3r
-      });
-
-  // Cartesian product of a tetrahedron and a 2D cone.
-  // The tetrahedron is empty and has vertices at
-  // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 100), and (100, 100 - 1/3, 100).
-  // The cone has apex at (10000, 10000) and
-  // edges passing through (1/3, 0) and (2/3, 0).
-  // Since the tetrahedron is empty, the Cartesian product is too.
-  checkPermutationsSample(false /* empty */, 5,
-                          {
-                              // Tetrahedron contraints:
-                              {0, 1, 0, 0, 0, 0},
-                              {0, -300, 299, 0, 0, 0},
-                              {300 * 299, -89400, -299, 0, 0, -100 * 299},
-                              {-897, 894, 0, 0, 0, 598},
-
-                              // Triangle constraints:
-                              // 300000p - 299999q >= 100000
-                              {0, 0, 0, 300000, -299999, -100000},
-                              // -300000p + 299998q + 200000 >= 0
-                              {0, 0, 0, -300000, 299998, 200000},
-                          },
-                          {});
-
-  // Cartesian product of same tetrahedron as above and
-  // {(p, q) : 1/3 <= p <= 2/3}.
-  checkPermutationsSample(false /* empty */, 5,
-                          {
-                              // Tetrahedron contraints:
-                              {0, 1, 0, 0, 0, 0},
-                              {0, -300, 299, 0, 0, 0},
-                              {300 * 299, -89400, -299, 0, 0, -100 * 299},
-                              {-897, 894, 0, 0, 0, 598},
-
-                              // Second set constraints:
-                              // 3p >= 1
-                              {0, 0, 0, 3, 0, -1},
-                              // 3p <= 2
-                              {0, 0, 0, -3, 0, 2},
-                          },
-                          {});
-
-  checkSample(true, parseFAC("(x, y, z) : (2 * x - 1 >= 0, x - y - 1 == 0, "
-                             "y - z == 0)",
-                             &context));
-
-  // Regression tests for the computation of dual coefficients.
-  checkSample(false, parseFAC("(x, y, z) : ("
-                              "6*x - 4*y + 9*z + 2 >= 0,"
-                              "x + 5*y + z + 5 >= 0,"
-                              "-4*x + y + 2*z - 1 >= 0,"
-                              "-3*x - 2*y - 7*z - 1 >= 0,"
-                              "-7*x - 5*y - 9*z - 1 >= 0)",
-                              &context));
-  checkSample(true, parseFAC("(x, y, z) : ("
-                             "3*x + 3*y + 3 >= 0,"
-                             "-4*x - 8*y - z + 4 >= 0,"
-                             "-7*x - 4*y + z + 1 >= 0,"
-                             "2*x - 7*y - 8*z - 7 >= 0,"
-                             "9*x + 8*y - 9*z - 7 >= 0)",
-                             &context));
-}
-
-TEST(FlatAffineConstraintsTest, IsIntegerEmptyTest) {
-
-  MLIRContext context;
-
-  // 1 <= 5x and 5x <= 4 (no solution).
-  EXPECT_TRUE(parseFAC("(x) : (5 * x - 1 >= 0, -5 * x + 4 >= 0)", &context)
-                  .isIntegerEmpty());
-  // 1 <= 5x and 5x <= 9 (solution: x = 1).
-  EXPECT_FALSE(parseFAC("(x) : (5 * x - 1 >= 0, -5 * x + 9 >= 0)", &context)
-                   .isIntegerEmpty());
-
-  // Unbounded sets.
-  EXPECT_TRUE(parseFAC("(x,y,z) : (2 * y - 1 >= 0, -2 * y + 1 >= 0, "
-                       "2 * z - 1 >= 0, 2 * x - 1 == 0)",
-                       &context)
-                  .isIntegerEmpty());
-
-  EXPECT_FALSE(parseFAC("(x,y,z) : (2 * x - 1 >= 0, -3 * x + 3 >= 0, "
-                        "5 * z - 6 >= 0, -7 * z + 17 >= 0, 3 * y - 2 >= 0)",
-                        &context)
-                   .isIntegerEmpty());
-
-  EXPECT_FALSE(
-      parseFAC("(x,y,z) : (2 * x - 1 >= 0, x - y - 1 == 0, y - z == 0)",
-               &context)
-          .isIntegerEmpty());
-
-  // FlatAffineConstraints::isEmpty() does not detect the following sets to be
-  // empty.
-
-  // 3x + 7y = 1 and 0 <= x, y <= 10.
-  // Since x and y are non-negative, 3x + 7y can never be 1.
-  EXPECT_TRUE(parseFAC("(x,y) : (x >= 0, -x + 10 >= 0, y >= 0, -y + 10 >= 0, "
-                       "3 * x + 7 * y - 1 == 0)",
-                       &context)
-                  .isIntegerEmpty());
-
-  // 2x = 3y and y = x - 1 and x + y = 6z + 2 and 0 <= x, y <= 100.
-  // Substituting y = x - 1 in 3y = 2x, we obtain x = 3 and hence y = 2.
-  // Since x + y = 5 cannot be equal to 6z + 2 for any z, the set is empty.
-  EXPECT_TRUE(
-      parseFAC("(x,y,z) : (x >= 0, -x + 100 >= 0, y >= 0, -y + 100 >= 0, "
-               "2 * x - 3 * y == 0, x - y - 1 == 0, x + y - 6 * z - 2 == 0)",
-               &context)
-          .isIntegerEmpty());
-
-  // 2x = 3y and y = x - 1 + 6z and x + y = 6q + 2 and 0 <= x, y <= 100.
-  // 2x = 3y implies x is a multiple of 3 and y is even.
-  // Now y = x - 1 + 6z implies y = 2 mod 3. In fact, since y is even, we have
-  // y = 2 mod 6. Then since x = y + 1 + 6z, we have x = 3 mod 6, implying
-  // x + y = 5 mod 6, which contradicts x + y = 6q + 2, so the set is empty.
-  EXPECT_TRUE(
-      parseFAC(
-          "(x,y,z,q) : (x >= 0, -x + 100 >= 0, y >= 0, -y + 100 >= 0, "
-          "2 * x - 3 * y == 0, x - y + 6 * z - 1 == 0, x + y - 6 * q - 2 == 0)",
-          &context)
-          .isIntegerEmpty());
-
-  // Set with symbols.
-  EXPECT_FALSE(
-      parseFAC("(x)[s] : (x + s >= 0, x - s == 0)", &context).isIntegerEmpty());
-}
-
-TEST(FlatAffineConstraintsTest, removeRedundantConstraintsTest) {
-  MLIRContext context;
-
-  FlatAffineConstraints fac =
-      parseFAC("(x) : (x - 2 >= 0, -x + 2 >= 0, x - 2 == 0)", &context);
-  fac.removeRedundantConstraints();
-
-  // Both inequalities are redundant given the equality. Both have been removed.
-  EXPECT_EQ(fac.getNumInequalities(), 0u);
-  EXPECT_EQ(fac.getNumEqualities(), 1u);
-
-  FlatAffineConstraints fac2 =
-      parseFAC("(x,y) : (x - 3 >= 0, y - 2 >= 0, x - y == 0)", &context);
-  fac2.removeRedundantConstraints();
-
-  // The second inequality is redundant and should have been removed. The
-  // remaining inequality should be the first one.
-  EXPECT_EQ(fac2.getNumInequalities(), 1u);
-  EXPECT_THAT(fac2.getInequality(0), ElementsAre(1, 0, -3));
-  EXPECT_EQ(fac2.getNumEqualities(), 1u);
-
-  FlatAffineConstraints fac3 =
-      parseFAC("(x,y,z) : (x - y == 0, x - z == 0, y - z == 0)", &context);
-  fac3.removeRedundantConstraints();
-
-  // One of the three equalities can be removed.
-  EXPECT_EQ(fac3.getNumInequalities(), 0u);
-  EXPECT_EQ(fac3.getNumEqualities(), 2u);
-
-  FlatAffineConstraints fac4 =
-      parseFAC("(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q) : ("
-               "b - 1 >= 0,"
-               "-b + 500 >= 0,"
-               "-16 * d + f >= 0,"
-               "f - 1 >= 0,"
-               "-f + 998 >= 0,"
-               "16 * d - f + 15 >= 0,"
-               "-16 * e + g >= 0,"
-               "g - 1 >= 0,"
-               "-g + 998 >= 0,"
-               "16 * e - g + 15 >= 0,"
-               "h >= 0,"
-               "-h + 1 >= 0,"
-               "j - 1 >= 0,"
-               "-j + 500 >= 0,"
-               "-f + 16 * l + 15 >= 0,"
-               "f - 16 * l >= 0,"
-               "-16 * m + o >= 0,"
-               "o - 1 >= 0,"
-               "-o + 998 >= 0,"
-               "16 * m - o + 15 >= 0,"
-               "p >= 0,"
-               "-p + 1 >= 0,"
-               "-g - h + 8 * q + 8 >= 0,"
-               "-o - p + 8 * q + 8 >= 0,"
-               "o + p - 8 * q - 1 >= 0,"
-               "g + h - 8 * q - 1 >= 0,"
-               "-f + n >= 0,"
-               "f - n >= 0,"
-               "k - 10 >= 0,"
-               "-k + 10 >= 0,"
-               "i - 13 >= 0,"
-               "-i + 13 >= 0,"
-               "c - 10 >= 0,"
-               "-c + 10 >= 0,"
-               "a - 13 >= 0,"
-               "-a + 13 >= 0"
-               ")",
-               &context);
-
-  // The above is a large set of constraints without any redundant constraints,
-  // as verified by the Fourier-Motzkin based removeRedundantInequalities.
-  unsigned nIneq = fac4.getNumInequalities();
-  unsigned nEq = fac4.getNumEqualities();
-  fac4.removeRedundantInequalities();
-  ASSERT_EQ(fac4.getNumInequalities(), nIneq);
-  ASSERT_EQ(fac4.getNumEqualities(), nEq);
-  // Now we test that removeRedundantConstraints does not find any constraints
-  // to be redundant either.
-  fac4.removeRedundantConstraints();
-  EXPECT_EQ(fac4.getNumInequalities(), nIneq);
-  EXPECT_EQ(fac4.getNumEqualities(), nEq);
-
-  FlatAffineConstraints fac5 = parseFAC(
-      "(x,y) : (128 * x + 127 >= 0, -x + 7 >= 0, -128 * x + y >= 0, y >= 0)",
-      &context);
-  // 128x + 127 >= 0  implies that 128x >= 0, since x has to be an integer.
-  // (This should be caught by GCDTightenInqualities().)
-  // So -128x + y >= 0 and 128x + 127 >= 0 imply y >= 0 since we have
-  // y >= 128x >= 0.
-  fac5.removeRedundantConstraints();
-  EXPECT_EQ(fac5.getNumInequalities(), 3u);
-  SmallVector<int64_t, 8> redundantConstraint = {0, 1, 0};
-  for (unsigned i = 0; i < 3; ++i) {
-    // Ensure that the removed constraint was the redundant constraint [3].
-    EXPECT_NE(fac5.getInequality(i), ArrayRef<int64_t>(redundantConstraint));
-  }
-}
-
-TEST(FlatAffineConstraintsTest, addConstantUpperBound) {
-  FlatAffineConstraints fac(2);
-  fac.addBound(FlatAffineConstraints::UB, 0, 1);
-  EXPECT_EQ(fac.atIneq(0, 0), -1);
-  EXPECT_EQ(fac.atIneq(0, 1), 0);
-  EXPECT_EQ(fac.atIneq(0, 2), 1);
-
-  fac.addBound(FlatAffineConstraints::UB, {1, 2, 3}, 1);
-  EXPECT_EQ(fac.atIneq(1, 0), -1);
-  EXPECT_EQ(fac.atIneq(1, 1), -2);
-  EXPECT_EQ(fac.atIneq(1, 2), -2);
-}
-
-TEST(FlatAffineConstraintsTest, addConstantLowerBound) {
-  FlatAffineConstraints fac(2);
-  fac.addBound(FlatAffineConstraints::LB, 0, 1);
-  EXPECT_EQ(fac.atIneq(0, 0), 1);
-  EXPECT_EQ(fac.atIneq(0, 1), 0);
-  EXPECT_EQ(fac.atIneq(0, 2), -1);
-
-  fac.addBound(FlatAffineConstraints::LB, {1, 2, 3}, 1);
-  EXPECT_EQ(fac.atIneq(1, 0), 1);
-  EXPECT_EQ(fac.atIneq(1, 1), 2);
-  EXPECT_EQ(fac.atIneq(1, 2), 2);
-}
-
-/// Check if the expected division representation of local variables matches the
-/// computed representation. The expected division representation is given as
-/// a vector of expressions set in `expectedDividends` and the corressponding
-/// denominator in `expectedDenominators`. The `denominators` and `dividends`
-/// obtained through `getLocalRepr` function is verified against the
-/// `expectedDenominators` and `expectedDividends` respectively.
-static void checkDivisionRepresentation(
-    FlatAffineConstraints &fac,
-    const std::vector<SmallVector<int64_t, 8>> &expectedDividends,
-    const SmallVectorImpl<unsigned> &expectedDenominators) {
-
-  std::vector<SmallVector<int64_t, 8>> dividends;
-  SmallVector<unsigned, 4> denominators;
-
-  fac.getLocalReprs(dividends, denominators);
-
-  // Check that the `denominators` and `expectedDenominators` match.
-  EXPECT_TRUE(expectedDenominators == denominators);
-
-  // Check that the `dividends` and `expectedDividends` match. If the
-  // denominator for a division is zero, we ignore its dividend.
-  EXPECT_TRUE(dividends.size() == expectedDividends.size());
-  for (unsigned i = 0, e = dividends.size(); i < e; ++i)
-    if (denominators[i] != 0)
-      EXPECT_TRUE(expectedDividends[i] == dividends[i]);
-}
-
-TEST(FlatAffineConstraintsTest, computeLocalReprSimple) {
-  FlatAffineConstraints fac(1);
-
-  fac.addLocalFloorDiv({1, 4}, 10);
-  fac.addLocalFloorDiv({1, 0, 100}, 10);
-
-  std::vector<SmallVector<int64_t, 8>> divisions = {{1, 0, 0, 4},
-                                                    {1, 0, 0, 100}};
-  SmallVector<unsigned, 8> denoms = {10, 10};
-
-  // Check if floordivs can be computed when no other inequalities exist
-  // and floor divs do not depend on each other.
-  checkDivisionRepresentation(fac, divisions, denoms);
-}
-
-TEST(FlatAffineConstraintsTest, computeLocalReprConstantFloorDiv) {
-  FlatAffineConstraints fac(4);
-
-  fac.addInequality({1, 0, 3, 1, 2});
-  fac.addInequality({1, 2, -8, 1, 10});
-  fac.addEquality({1, 2, -4, 1, 10});
-
-  fac.addLocalFloorDiv({0, 0, 0, 0, 100}, 30);
-  fac.addLocalFloorDiv({0, 0, 0, 0, 0, 206}, 101);
-
-  std::vector<SmallVector<int64_t, 8>> divisions = {{0, 0, 0, 0, 0, 0, 3},
-                                                    {0, 0, 0, 0, 0, 0, 2}};
-  SmallVector<unsigned, 8> denoms = {1, 1};
-
-  // Check if floordivs with constant numerator can be computed.
-  checkDivisionRepresentation(fac, divisions, denoms);
-}
-
-TEST(FlatAffineConstraintsTest, computeLocalReprRecursive) {
-  FlatAffineConstraints fac(4);
-  fac.addInequality({1, 0, 3, 1, 2});
-  fac.addInequality({1, 2, -8, 1, 10});
-  fac.addEquality({1, 2, -4, 1, 10});
-
-  fac.addLocalFloorDiv({0, -2, 7, 2, 10}, 3);
-  fac.addLocalFloorDiv({3, 0, 9, 2, 2, 10}, 5);
-  fac.addLocalFloorDiv({0, 1, -123, 2, 0, -4, 10}, 3);
-
-  fac.addInequality({1, 2, -2, 1, -5, 0, 6, 100});
-  fac.addInequality({1, 2, -8, 1, 3, 7, 0, -9});
-
-  std::vector<SmallVector<int64_t, 8>> divisions = {
-      {0, -2, 7, 2, 0, 0, 0, 10},
-      {3, 0, 9, 2, 2, 0, 0, 10},
-      {0, 1, -123, 2, 0, -4, 0, 10}};
-
-  SmallVector<unsigned, 8> denoms = {3, 5, 3};
-
-  // Check if floordivs which may depend on other floordivs can be computed.
-  checkDivisionRepresentation(fac, divisions, denoms);
-}
-
-TEST(FlatAffineConstraintsTest, computeLocalReprTightUpperBound) {
-  MLIRContext context;
-
-  {
-    FlatAffineConstraints fac = parseFAC("(i) : (i mod 3 - 1 >= 0)", &context);
-
-    // The set formed by the fac is:
-    //        3q - i + 2 >= 0             <-- Division lower bound
-    //       -3q + i - 1 >= 0
-    //       -3q + i     >= 0             <-- Division upper bound
-    // We remove redundant constraints to get the set:
-    //        3q - i + 2 >= 0             <-- Division lower bound
-    //       -3q + i - 1 >= 0             <-- Tighter division upper bound
-    // thus, making the upper bound tighter.
-    fac.removeRedundantConstraints();
-
-    std::vector<SmallVector<int64_t, 8>> divisions = {{1, 0, 0}};
-    SmallVector<unsigned, 8> denoms = {3};
-
-    // Check if the divisions can be computed even with a tighter upper bound.
-    checkDivisionRepresentation(fac, divisions, denoms);
-  }
-
-  {
-    FlatAffineConstraints fac = parseFAC(
-        "(i, j, q) : (4*q - i - j + 2 >= 0, -4*q + i + j >= 0)", &context);
-    // Convert `q` to a local variable.
-    fac.convertDimToLocal(2, 3);
-
-    std::vector<SmallVector<int64_t, 8>> divisions = {{1, 1, 0, 1}};
-    SmallVector<unsigned, 8> denoms = {4};
-
-    // Check if the divisions can be computed even with a tighter upper bound.
-    checkDivisionRepresentation(fac, divisions, denoms);
-  }
-}
-
-TEST(FlatAffineConstraintsTest, computeLocalReprNoRepr) {
-  MLIRContext context;
-  FlatAffineConstraints fac =
-      parseFAC("(x, q) : (x - 3 * q >= 0, -x + 3 * q + 3 >= 0)", &context);
-  // Convert q to a local variable.
-  fac.convertDimToLocal(1, 2);
-
-  std::vector<SmallVector<int64_t, 8>> divisions = {{0, 0, 0}};
-  SmallVector<unsigned, 8> denoms = {0};
-
-  // Check that no division is computed.
-  checkDivisionRepresentation(fac, divisions, denoms);
-}
-
-TEST(FlatAffineConstraintsTest, simplifyLocalsTest) {
-  // (x) : (exists y: 2x + y = 1 and y = 2).
-  FlatAffineConstraints fac(1, 0, 1);
-  fac.addEquality({2, 1, -1});
-  fac.addEquality({0, 1, -2});
-
-  EXPECT_TRUE(fac.isEmpty());
-
-  // (x) : (exists y, z, w: 3x + y = 1 and 2y = z and 3y = w and z = w).
-  FlatAffineConstraints fac2(1, 0, 3);
-  fac2.addEquality({3, 1, 0, 0, -1});
-  fac2.addEquality({0, 2, -1, 0, 0});
-  fac2.addEquality({0, 3, 0, -1, 0});
-  fac2.addEquality({0, 0, 1, -1, 0});
-
-  EXPECT_TRUE(fac2.isEmpty());
-
-  // (x) : (exists y: x >= y + 1 and 2x + y = 0 and y >= -1).
-  FlatAffineConstraints fac3(1, 0, 1);
-  fac3.addInequality({1, -1, -1});
-  fac3.addInequality({0, 1, 1});
-  fac3.addEquality({2, 1, 0});
-
-  EXPECT_TRUE(fac3.isEmpty());
-}
-
-TEST(FlatAffineConstraintsTest, mergeDivisionsSimple) {
-  {
-    // (x) : (exists z, y  = [x / 2] : x = 3y and x + z + 1 >= 0).
-    FlatAffineConstraints fac1(1, 0, 1);
-    fac1.addLocalFloorDiv({1, 0, 0}, 2); // y = [x / 2].
-    fac1.addEquality({1, 0, -3, 0});     // x = 3y.
-    fac1.addInequality({1, 1, 0, 1});    // x + z + 1 >= 0.
-
-    // (x) : (exists y = [x / 2], z : x = 5y).
-    FlatAffineConstraints fac2(1);
-    fac2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2].
-    fac2.addEquality({1, -5, 0});     // x = 5y.
-    fac2.appendLocalId();             // Add local id z.
-
-    fac1.mergeLocalIds(fac2);
-
-    // Local space should be same.
-    EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds());
-
-    // 1 division should be matched + 2 unmatched local ids.
-    EXPECT_EQ(fac1.getNumLocalIds(), 3u);
-    EXPECT_EQ(fac2.getNumLocalIds(), 3u);
-  }
-
-  {
-    // (x) : (exists z = [x / 5], y = [x / 2] : x = 3y).
-    FlatAffineConstraints fac1(1);
-    fac1.addLocalFloorDiv({1, 0}, 5);    // z = [x / 5].
-    fac1.addLocalFloorDiv({1, 0, 0}, 2); // y = [x / 2].
-    fac1.addEquality({1, 0, -3, 0});     // x = 3y.
-
-    // (x) : (exists y = [x / 2], z = [x / 5]: x = 5z).
-    FlatAffineConstraints fac2(1);
-    fac2.addLocalFloorDiv({1, 0}, 2);    // y = [x / 2].
-    fac2.addLocalFloorDiv({1, 0, 0}, 5); // z = [x / 5].
-    fac2.addEquality({1, 0, -5, 0});     // x = 5z.
-
-    fac1.mergeLocalIds(fac2);
-
-    // Local space should be same.
-    EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds());
-
-    // 2 divisions should be matched.
-    EXPECT_EQ(fac1.getNumLocalIds(), 2u);
-    EXPECT_EQ(fac2.getNumLocalIds(), 2u);
-  }
-
-  {
-    // Division Normalization test.
-    // (x) : (exists z, y  = [x / 2] : x = 3y and x + z + 1 >= 0).
-    FlatAffineConstraints fac1(1, 0, 1);
-    // This division would be normalized.
-    fac1.addLocalFloorDiv({3, 0, 0}, 6); // y = [3x / 6] -> [x/2].
-    fac1.addEquality({1, 0, -3, 0});     // x = 3z.
-    fac1.addInequality({1, 1, 0, 1});    // x + y + 1 >= 0.
-
-    // (x) : (exists y = [x / 2], z : x = 5y).
-    FlatAffineConstraints fac2(1);
-    fac2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2].
-    fac2.addEquality({1, -5, 0});     // x = 5y.
-    fac2.appendLocalId();             // Add local id z.
-
-    fac1.mergeLocalIds(fac2);
-
-    // Local space should be same.
-    EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds());
-
-    // One division should be matched + 2 unmatched local ids.
-    EXPECT_EQ(fac1.getNumLocalIds(), 3u);
-    EXPECT_EQ(fac2.getNumLocalIds(), 3u);
-  }
-}
-
-TEST(FlatAffineConstraintsTest, mergeDivisionsNestedDivsions) {
-  {
-    // (x) : (exists y = [x / 2], z = [x + y / 3]: y + z >= x).
-    FlatAffineConstraints fac1(1);
-    fac1.addLocalFloorDiv({1, 0}, 2);    // y = [x / 2].
-    fac1.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3].
-    fac1.addInequality({-1, 1, 1, 0});   // y + z >= x.
-
-    // (x) : (exists y = [x / 2], z = [x + y / 3]: y + z <= x).
-    FlatAffineConstraints fac2(1);
-    fac2.addLocalFloorDiv({1, 0}, 2);    // y = [x / 2].
-    fac2.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3].
-    fac2.addInequality({1, -1, -1, 0});  // y + z <= x.
-
-    fac1.mergeLocalIds(fac2);
-
-    // Local space should be same.
-    EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds());
-
-    // 2 divisions should be matched.
-    EXPECT_EQ(fac1.getNumLocalIds(), 2u);
-    EXPECT_EQ(fac2.getNumLocalIds(), 2u);
-  }
-
-  {
-    // (x) : (exists y = [x / 2], z = [x + y / 3], w = [z + 1 / 5]: y + z >= x).
-    FlatAffineConstraints fac1(1);
-    fac1.addLocalFloorDiv({1, 0}, 2);       // y = [x / 2].
-    fac1.addLocalFloorDiv({1, 1, 0}, 3);    // z = [x + y / 3].
-    fac1.addLocalFloorDiv({0, 0, 1, 1}, 5); // w = [z + 1 / 5].
-    fac1.addInequality({-1, 1, 1, 0, 0});   // y + z >= x.
-
-    // (x) : (exists y = [x / 2], z = [x + y / 3], w = [z + 1 / 5]: y + z <= x).
-    FlatAffineConstraints fac2(1);
-    fac2.addLocalFloorDiv({1, 0}, 2);       // y = [x / 2].
-    fac2.addLocalFloorDiv({1, 1, 0}, 3);    // z = [x + y / 3].
-    fac2.addLocalFloorDiv({0, 0, 1, 1}, 5); // w = [z + 1 / 5].
-    fac2.addInequality({1, -1, -1, 0, 0});  // y + z <= x.
-
-    fac1.mergeLocalIds(fac2);
-
-    // Local space should be same.
-    EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds());
-
-    // 3 divisions should be matched.
-    EXPECT_EQ(fac1.getNumLocalIds(), 3u);
-    EXPECT_EQ(fac2.getNumLocalIds(), 3u);
-  }
-  {
-    // (x) : (exists y = [x / 2], z = [x + y / 3]: y + z >= x).
-    FlatAffineConstraints fac1(1);
-    fac1.addLocalFloorDiv({2, 0}, 4);    // y = [2x / 4] -> [x / 2].
-    fac1.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3].
-    fac1.addInequality({-1, 1, 1, 0});   // y + z >= x.
-
-    // (x) : (exists y = [x / 2], z = [x + y / 3]: y + z <= x).
-    FlatAffineConstraints fac2(1);
-    fac2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2].
-    // This division would be normalized.
-    fac2.addLocalFloorDiv({3, 3, 0}, 9); // z = [3x + 3y / 9] -> [x + y / 3].
-    fac2.addInequality({1, -1, -1, 0});  // y + z <= x.
-
-    fac1.mergeLocalIds(fac2);
-
-    // Local space should be same.
-    EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds());
-
-    // 2 divisions should be matched.
-    EXPECT_EQ(fac1.getNumLocalIds(), 2u);
-    EXPECT_EQ(fac2.getNumLocalIds(), 2u);
-  }
-}
-
-TEST(FlatAffineConstraintsTest, mergeDivisionsConstants) {
-  {
-    // (x) : (exists y = [x + 1 / 3], z = [x + 2 / 3]: y + z >= x).
-    FlatAffineConstraints fac1(1);
-    fac1.addLocalFloorDiv({1, 1}, 2);    // y = [x + 1 / 2].
-    fac1.addLocalFloorDiv({1, 0, 2}, 3); // z = [x + 2 / 3].
-    fac1.addInequality({-1, 1, 1, 0});   // y + z >= x.
-
-    // (x) : (exists y = [x + 1 / 3], z = [x + 2 / 3]: y + z <= x).
-    FlatAffineConstraints fac2(1);
-    fac2.addLocalFloorDiv({1, 1}, 2);    // y = [x + 1 / 2].
-    fac2.addLocalFloorDiv({1, 0, 2}, 3); // z = [x + 2 / 3].
-    fac2.addInequality({1, -1, -1, 0});  // y + z <= x.
-
-    fac1.mergeLocalIds(fac2);
-
-    // Local space should be same.
-    EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds());
-
-    // 2 divisions should be matched.
-    EXPECT_EQ(fac1.getNumLocalIds(), 2u);
-    EXPECT_EQ(fac2.getNumLocalIds(), 2u);
-  }
-  {
-    // (x) : (exists y = [x + 1 / 3], z = [x + 2 / 3]: y + z >= x).
-    FlatAffineConstraints fac1(1);
-    fac1.addLocalFloorDiv({1, 1}, 2); // y = [x + 1 / 2].
-    // Normalization test.
-    fac1.addLocalFloorDiv({3, 0, 6}, 9); // z = [3x + 6 / 9] -> [x + 2 / 3].
-    fac1.addInequality({-1, 1, 1, 0});   // y + z >= x.
-
-    // (x) : (exists y = [x + 1 / 3], z = [x + 2 / 3]: y + z <= x).
-    FlatAffineConstraints fac2(1);
-    // Normalization test.
-    fac2.addLocalFloorDiv({2, 2}, 4);    // y = [2x + 2 / 4] -> [x + 1 / 2].
-    fac2.addLocalFloorDiv({1, 0, 2}, 3); // z = [x + 2 / 3].
-    fac2.addInequality({1, -1, -1, 0});  // y + z <= x.
-
-    fac1.mergeLocalIds(fac2);
-
-    // Local space should be same.
-    EXPECT_EQ(fac1.getNumLocalIds(), fac2.getNumLocalIds());
-
-    // 2 divisions should be matched.
-    EXPECT_EQ(fac1.getNumLocalIds(), 2u);
-    EXPECT_EQ(fac2.getNumLocalIds(), 2u);
-  }
-}
-
-} // namespace mlir

diff  --git a/mlir/unittests/Analysis/CMakeLists.txt b/mlir/unittests/Analysis/CMakeLists.txt
index b70d6822ac87..c5e0beb3e5fe 100644
--- a/mlir/unittests/Analysis/CMakeLists.txt
+++ b/mlir/unittests/Analysis/CMakeLists.txt
@@ -1,7 +1,6 @@
 add_mlir_unittest(MLIRAnalysisTests
   AffineStructuresParser.cpp
   AffineStructuresParserTest.cpp
-  AffineStructuresTest.cpp
   PresburgerSetTest.cpp
 )
 

diff  --git a/mlir/unittests/Analysis/Presburger/IntegerPolyhedronTest.cpp b/mlir/unittests/Analysis/Presburger/IntegerPolyhedronTest.cpp
index 1598f051678f..15a9e1461ab2 100644
--- a/mlir/unittests/Analysis/Presburger/IntegerPolyhedronTest.cpp
+++ b/mlir/unittests/Analysis/Presburger/IntegerPolyhedronTest.cpp
@@ -7,14 +7,20 @@
 //===----------------------------------------------------------------------===//
 
 #include "mlir/Analysis/Presburger/IntegerPolyhedron.h"
+#include "../AffineStructuresParser.h"
+#include "mlir/IR/MLIRContext.h"
 
 #include <gmock/gmock.h>
 #include <gtest/gtest.h>
 
+#include <numeric>
+
 namespace mlir {
 
 using testing::ElementsAre;
 
+enum class TestFunction { Sample, Empty };
+
 /// Construct a IntegerPolyhedron from a set of inequality and
 /// equality constraints.
 static IntegerPolyhedron
@@ -30,6 +36,79 @@ makeSetFromConstraints(unsigned ids, ArrayRef<SmallVector<int64_t, 4>> ineqs,
   return set;
 }
 
+/// If fn is TestFunction::Sample (default):
+/// If hasSample is true, check that findIntegerSample returns a valid sample
+/// for the IntegerPolyhedron poly.
+/// If hasSample is false, check that findIntegerSample returns None.
+///
+/// If fn is TestFunction::Empty, check that isIntegerEmpty returns the
+/// opposite of hasSample.
+static void checkSample(bool hasSample, const IntegerPolyhedron &poly,
+                        TestFunction fn = TestFunction::Sample) {
+  Optional<SmallVector<int64_t, 8>> maybeSample;
+  switch (fn) {
+  case TestFunction::Sample:
+    maybeSample = poly.findIntegerSample();
+    if (!hasSample) {
+      EXPECT_FALSE(maybeSample.hasValue());
+      if (maybeSample.hasValue()) {
+        for (auto x : *maybeSample)
+          llvm::errs() << x << ' ';
+        llvm::errs() << '\n';
+      }
+    } else {
+      ASSERT_TRUE(maybeSample.hasValue());
+      EXPECT_TRUE(poly.containsPoint(*maybeSample));
+    }
+    break;
+  case TestFunction::Empty:
+    EXPECT_EQ(!hasSample, poly.isIntegerEmpty());
+    break;
+  }
+}
+
+/// Check sampling for all the permutations of the dimensions for the given
+/// constraint set. Since the GBR algorithm progresses dimension-wise, 
diff erent
+/// orderings may cause the algorithm to proceed 
diff erently. At least some of
+///.these permutations should make it past the heuristics and test the
+/// implementation of the GBR algorithm itself.
+/// Use TestFunction fn to test.
+static void checkPermutationsSample(bool hasSample, unsigned nDim,
+                                    ArrayRef<SmallVector<int64_t, 4>> ineqs,
+                                    ArrayRef<SmallVector<int64_t, 4>> eqs,
+                                    TestFunction fn = TestFunction::Sample) {
+  SmallVector<unsigned, 4> perm(nDim);
+  std::iota(perm.begin(), perm.end(), 0);
+  auto permute = [&perm](ArrayRef<int64_t> coeffs) {
+    SmallVector<int64_t, 4> permuted;
+    for (unsigned id : perm)
+      permuted.push_back(coeffs[id]);
+    permuted.push_back(coeffs.back());
+    return permuted;
+  };
+  do {
+    SmallVector<SmallVector<int64_t, 4>, 4> permutedIneqs, permutedEqs;
+    for (const auto &ineq : ineqs)
+      permutedIneqs.push_back(permute(ineq));
+    for (const auto &eq : eqs)
+      permutedEqs.push_back(permute(eq));
+
+    checkSample(hasSample,
+                makeSetFromConstraints(nDim, permutedIneqs, permutedEqs), fn);
+  } while (std::next_permutation(perm.begin(), perm.end()));
+}
+
+/// Parses a IntegerPolyhedron from a StringRef. It is expected that the
+/// string represents a valid IntegerSet, otherwise it will violate a gtest
+/// assertion.
+static IntegerPolyhedron parsePoly(StringRef str, MLIRContext *context) {
+  FailureOr<IntegerPolyhedron> poly = parseIntegerSetToFAC(str, context);
+
+  EXPECT_TRUE(succeeded(poly));
+
+  return *poly;
+}
+
 TEST(IntegerPolyhedronTest, removeInequality) {
   IntegerPolyhedron set =
       makeSetFromConstraints(1, {{0, 0}, {1, 1}, {2, 2}, {3, 3}, {4, 4}}, {});
@@ -100,4 +179,794 @@ TEST(IntegerPolyhedronTest, removeIdRange) {
   EXPECT_THAT(set.getInequality(0), testing::ElementsAre(12, 20, 40));
 }
 
+TEST(IntegerPolyhedronTest, FindSampleTest) {
+  // Bounded sets with only inequalities.
+
+  MLIRContext context;
+
+  // 0 <= 7x <= 5
+  checkSample(true, parsePoly("(x) : (7 * x >= 0, -7 * x + 5 >= 0)", &context));
+
+  // 1 <= 5x and 5x <= 4 (no solution).
+  checkSample(false,
+              parsePoly("(x) : (5 * x - 1 >= 0, -5 * x + 4 >= 0)", &context));
+
+  // 1 <= 5x and 5x <= 9 (solution: x = 1).
+  checkSample(true,
+              parsePoly("(x) : (5 * x - 1 >= 0, -5 * x + 9 >= 0)", &context));
+
+  // Bounded sets with equalities.
+  // x >= 8 and 40 >= y and x = y.
+  checkSample(true, parsePoly("(x,y) : (x - 8 >= 0, -y + 40 >= 0, x - y == 0)",
+                              &context));
+
+  // x <= 10 and y <= 10 and 10 <= z and x + 2y = 3z.
+  // solution: x = y = z = 10.
+  checkSample(true, parsePoly("(x,y,z) : (-x + 10 >= 0, -y + 10 >= 0, "
+                              "z - 10 >= 0, x + 2 * y - 3 * z == 0)",
+                              &context));
+
+  // x <= 10 and y <= 10 and 11 <= z and x + 2y = 3z.
+  // This implies x + 2y >= 33 and x + 2y <= 30, which has no solution.
+  checkSample(false, parsePoly("(x,y,z) : (-x + 10 >= 0, -y + 10 >= 0, "
+                               "z - 11 >= 0, x + 2 * y - 3 * z == 0)",
+                               &context));
+
+  // 0 <= r and r <= 3 and 4q + r = 7.
+  // Solution: q = 1, r = 3.
+  checkSample(
+      true,
+      parsePoly("(q,r) : (r >= 0, -r + 3 >= 0, 4 * q + r - 7 == 0)", &context));
+
+  // 4q + r = 7 and r = 0.
+  // Solution: q = 1, r = 3.
+  checkSample(false,
+              parsePoly("(q,r) : (4 * q + r - 7 == 0, r == 0)", &context));
+
+  // The next two sets are large sets that should take a long time to sample
+  // with a naive branch and bound algorithm but can be sampled efficiently with
+  // the GBR algorithm.
+  //
+  // This is a triangle with vertices at (1/3, 0), (2/3, 0) and (10000, 10000).
+  checkSample(true, parsePoly("(x,y) : (y >= 0, "
+                              "300000 * x - 299999 * y - 100000 >= 0, "
+                              "-300000 * x + 299998 * y + 200000 >= 0)",
+                              &context));
+
+  // This is a tetrahedron with vertices at
+  // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 10000), and (10000, 10000, 10000).
+  // The first three points form a triangular base on the xz plane with the
+  // apex at the fourth point, which is the only integer point.
+  checkPermutationsSample(
+      true, 3,
+      {
+          {0, 1, 0, 0},  // y >= 0
+          {0, -1, 1, 0}, // z >= y
+          {300000, -299998, -1,
+           -100000},                    // -300000x + 299998y + 100000 + z <= 0.
+          {-150000, 149999, 0, 100000}, // -150000x + 149999y + 100000 >= 0.
+      },
+      {});
+
+  // Same thing with some spurious extra dimensions equated to constants.
+  checkSample(
+      true,
+      parsePoly("(a,b,c,d,e) : (b + d - e >= 0, -b + c - d + e >= 0, "
+                "300000 * a - 299998 * b - c - 9 * d + 21 * e - 112000 >= 0, "
+                "-150000 * a + 149999 * b - 15 * d + 47 * e + 68000 >= 0, "
+                "d - e == 0, d + e - 2000 == 0)",
+                &context));
+
+  // This is a tetrahedron with vertices at
+  // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 100), (100, 100 - 1/3, 100).
+  checkPermutationsSample(false, 3,
+                          {
+                              {0, 1, 0, 0},
+                              {0, -300, 299, 0},
+                              {300 * 299, -89400, -299, -100 * 299},
+                              {-897, 894, 0, 598},
+                          },
+                          {});
+
+  // Two tests involving equalities that are integer empty but not rational
+  // empty.
+
+  // This is a line segment from (0, 1/3) to (100, 100 + 1/3).
+  checkSample(
+      false,
+      parsePoly("(x,y) : (x >= 0, -x + 100 >= 0, 3 * x - 3 * y + 1 == 0)",
+                &context));
+
+  // A thin parallelogram. 0 <= x <= 100 and x + 1/3 <= y <= x + 2/3.
+  checkSample(false,
+              parsePoly("(x,y) : (x >= 0, -x + 100 >= 0, "
+                        "3 * x - 3 * y + 2 >= 0, -3 * x + 3 * y - 1 >= 0)",
+                        &context));
+
+  checkSample(true, parsePoly("(x,y) : (2 * x >= 0, -2 * x + 99 >= 0, "
+                              "2 * y >= 0, -2 * y + 99 >= 0)",
+                              &context));
+
+  // 2D cone with apex at (10000, 10000) and
+  // edges passing through (1/3, 0) and (2/3, 0).
+  checkSample(true, parsePoly("(x,y) : (300000 * x - 299999 * y - 100000 >= 0, "
+                              "-300000 * x + 299998 * y + 200000 >= 0)",
+                              &context));
+
+  // Cartesian product of a tetrahedron and a 2D cone.
+  // The tetrahedron has vertices at
+  // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 10000), and (10000, 10000, 10000).
+  // The first three points form a triangular base on the xz plane with the
+  // apex at the fourth point, which is the only integer point.
+  // The cone has apex at (10000, 10000) and
+  // edges passing through (1/3, 0) and (2/3, 0).
+  checkPermutationsSample(
+      true /* not empty */, 5,
+      {
+          // Tetrahedron contraints:
+          {0, 1, 0, 0, 0, 0},  // y >= 0
+          {0, -1, 1, 0, 0, 0}, // z >= y
+                               // -300000x + 299998y + 100000 + z <= 0.
+          {300000, -299998, -1, 0, 0, -100000},
+          // -150000x + 149999y + 100000 >= 0.
+          {-150000, 149999, 0, 0, 0, 100000},
+
+          // Triangle constraints:
+          // 300000p - 299999q >= 100000
+          {0, 0, 0, 300000, -299999, -100000},
+          // -300000p + 299998q + 200000 >= 0
+          {0, 0, 0, -300000, 299998, 200000},
+      },
+      {});
+
+  // Cartesian product of same tetrahedron as above and {(p, q) : 1/3 <= p <=
+  // 2/3}. Since the second set is empty, the whole set is too.
+  checkPermutationsSample(
+      false /* empty */, 5,
+      {
+          // Tetrahedron contraints:
+          {0, 1, 0, 0, 0, 0},  // y >= 0
+          {0, -1, 1, 0, 0, 0}, // z >= y
+                               // -300000x + 299998y + 100000 + z <= 0.
+          {300000, -299998, -1, 0, 0, -100000},
+          // -150000x + 149999y + 100000 >= 0.
+          {-150000, 149999, 0, 0, 0, 100000},
+
+          // Second set constraints:
+          // 3p >= 1
+          {0, 0, 0, 3, 0, -1},
+          // 3p <= 2
+          {0, 0, 0, -3, 0, 2},
+      },
+      {});
+
+  // Cartesian product of same tetrahedron as above and
+  // {(p, q, r) : 1 <= p <= 2 and p = 3q + 3r}.
+  // Since the second set is empty, the whole set is too.
+  checkPermutationsSample(
+      false /* empty */, 5,
+      {
+          // Tetrahedron contraints:
+          {0, 1, 0, 0, 0, 0, 0},  // y >= 0
+          {0, -1, 1, 0, 0, 0, 0}, // z >= y
+                                  // -300000x + 299998y + 100000 + z <= 0.
+          {300000, -299998, -1, 0, 0, 0, -100000},
+          // -150000x + 149999y + 100000 >= 0.
+          {-150000, 149999, 0, 0, 0, 0, 100000},
+
+          // Second set constraints:
+          // p >= 1
+          {0, 0, 0, 1, 0, 0, -1},
+          // p <= 2
+          {0, 0, 0, -1, 0, 0, 2},
+      },
+      {
+          {0, 0, 0, 1, -3, -3, 0}, // p = 3q + 3r
+      });
+
+  // Cartesian product of a tetrahedron and a 2D cone.
+  // The tetrahedron is empty and has vertices at
+  // (1/3, 0, 0), (2/3, 0, 0), (2/3, 0, 100), and (100, 100 - 1/3, 100).
+  // The cone has apex at (10000, 10000) and
+  // edges passing through (1/3, 0) and (2/3, 0).
+  // Since the tetrahedron is empty, the Cartesian product is too.
+  checkPermutationsSample(false /* empty */, 5,
+                          {
+                              // Tetrahedron contraints:
+                              {0, 1, 0, 0, 0, 0},
+                              {0, -300, 299, 0, 0, 0},
+                              {300 * 299, -89400, -299, 0, 0, -100 * 299},
+                              {-897, 894, 0, 0, 0, 598},
+
+                              // Triangle constraints:
+                              // 300000p - 299999q >= 100000
+                              {0, 0, 0, 300000, -299999, -100000},
+                              // -300000p + 299998q + 200000 >= 0
+                              {0, 0, 0, -300000, 299998, 200000},
+                          },
+                          {});
+
+  // Cartesian product of same tetrahedron as above and
+  // {(p, q) : 1/3 <= p <= 2/3}.
+  checkPermutationsSample(false /* empty */, 5,
+                          {
+                              // Tetrahedron contraints:
+                              {0, 1, 0, 0, 0, 0},
+                              {0, -300, 299, 0, 0, 0},
+                              {300 * 299, -89400, -299, 0, 0, -100 * 299},
+                              {-897, 894, 0, 0, 0, 598},
+
+                              // Second set constraints:
+                              // 3p >= 1
+                              {0, 0, 0, 3, 0, -1},
+                              // 3p <= 2
+                              {0, 0, 0, -3, 0, 2},
+                          },
+                          {});
+
+  checkSample(true, parsePoly("(x, y, z) : (2 * x - 1 >= 0, x - y - 1 == 0, "
+                              "y - z == 0)",
+                              &context));
+
+  // Regression tests for the computation of dual coefficients.
+  checkSample(false, parsePoly("(x, y, z) : ("
+                               "6*x - 4*y + 9*z + 2 >= 0,"
+                               "x + 5*y + z + 5 >= 0,"
+                               "-4*x + y + 2*z - 1 >= 0,"
+                               "-3*x - 2*y - 7*z - 1 >= 0,"
+                               "-7*x - 5*y - 9*z - 1 >= 0)",
+                               &context));
+  checkSample(true, parsePoly("(x, y, z) : ("
+                              "3*x + 3*y + 3 >= 0,"
+                              "-4*x - 8*y - z + 4 >= 0,"
+                              "-7*x - 4*y + z + 1 >= 0,"
+                              "2*x - 7*y - 8*z - 7 >= 0,"
+                              "9*x + 8*y - 9*z - 7 >= 0)",
+                              &context));
+}
+
+TEST(IntegerPolyhedronTest, IsIntegerEmptyTest) {
+
+  MLIRContext context;
+
+  // 1 <= 5x and 5x <= 4 (no solution).
+  EXPECT_TRUE(parsePoly("(x) : (5 * x - 1 >= 0, -5 * x + 4 >= 0)", &context)
+                  .isIntegerEmpty());
+  // 1 <= 5x and 5x <= 9 (solution: x = 1).
+  EXPECT_FALSE(parsePoly("(x) : (5 * x - 1 >= 0, -5 * x + 9 >= 0)", &context)
+                   .isIntegerEmpty());
+
+  // Unbounded sets.
+  EXPECT_TRUE(parsePoly("(x,y,z) : (2 * y - 1 >= 0, -2 * y + 1 >= 0, "
+                        "2 * z - 1 >= 0, 2 * x - 1 == 0)",
+                        &context)
+                  .isIntegerEmpty());
+
+  EXPECT_FALSE(parsePoly("(x,y,z) : (2 * x - 1 >= 0, -3 * x + 3 >= 0, "
+                         "5 * z - 6 >= 0, -7 * z + 17 >= 0, 3 * y - 2 >= 0)",
+                         &context)
+                   .isIntegerEmpty());
+
+  EXPECT_FALSE(
+      parsePoly("(x,y,z) : (2 * x - 1 >= 0, x - y - 1 == 0, y - z == 0)",
+                &context)
+          .isIntegerEmpty());
+
+  // IntegerPolyhedron::isEmpty() does not detect the following sets to be
+  // empty.
+
+  // 3x + 7y = 1 and 0 <= x, y <= 10.
+  // Since x and y are non-negative, 3x + 7y can never be 1.
+  EXPECT_TRUE(parsePoly("(x,y) : (x >= 0, -x + 10 >= 0, y >= 0, -y + 10 >= 0, "
+                        "3 * x + 7 * y - 1 == 0)",
+                        &context)
+                  .isIntegerEmpty());
+
+  // 2x = 3y and y = x - 1 and x + y = 6z + 2 and 0 <= x, y <= 100.
+  // Substituting y = x - 1 in 3y = 2x, we obtain x = 3 and hence y = 2.
+  // Since x + y = 5 cannot be equal to 6z + 2 for any z, the set is empty.
+  EXPECT_TRUE(
+      parsePoly("(x,y,z) : (x >= 0, -x + 100 >= 0, y >= 0, -y + 100 >= 0, "
+                "2 * x - 3 * y == 0, x - y - 1 == 0, x + y - 6 * z - 2 == 0)",
+                &context)
+          .isIntegerEmpty());
+
+  // 2x = 3y and y = x - 1 + 6z and x + y = 6q + 2 and 0 <= x, y <= 100.
+  // 2x = 3y implies x is a multiple of 3 and y is even.
+  // Now y = x - 1 + 6z implies y = 2 mod 3. In fact, since y is even, we have
+  // y = 2 mod 6. Then since x = y + 1 + 6z, we have x = 3 mod 6, implying
+  // x + y = 5 mod 6, which contradicts x + y = 6q + 2, so the set is empty.
+  EXPECT_TRUE(
+      parsePoly(
+          "(x,y,z,q) : (x >= 0, -x + 100 >= 0, y >= 0, -y + 100 >= 0, "
+          "2 * x - 3 * y == 0, x - y + 6 * z - 1 == 0, x + y - 6 * q - 2 == 0)",
+          &context)
+          .isIntegerEmpty());
+
+  // Set with symbols.
+  EXPECT_FALSE(parsePoly("(x)[s] : (x + s >= 0, x - s == 0)", &context)
+                   .isIntegerEmpty());
+}
+
+TEST(IntegerPolyhedronTest, removeRedundantConstraintsTest) {
+  MLIRContext context;
+
+  IntegerPolyhedron poly =
+      parsePoly("(x) : (x - 2 >= 0, -x + 2 >= 0, x - 2 == 0)", &context);
+  poly.removeRedundantConstraints();
+
+  // Both inequalities are redundant given the equality. Both have been removed.
+  EXPECT_EQ(poly.getNumInequalities(), 0u);
+  EXPECT_EQ(poly.getNumEqualities(), 1u);
+
+  IntegerPolyhedron poly2 =
+      parsePoly("(x,y) : (x - 3 >= 0, y - 2 >= 0, x - y == 0)", &context);
+  poly2.removeRedundantConstraints();
+
+  // The second inequality is redundant and should have been removed. The
+  // remaining inequality should be the first one.
+  EXPECT_EQ(poly2.getNumInequalities(), 1u);
+  EXPECT_THAT(poly2.getInequality(0), ElementsAre(1, 0, -3));
+  EXPECT_EQ(poly2.getNumEqualities(), 1u);
+
+  IntegerPolyhedron poly3 =
+      parsePoly("(x,y,z) : (x - y == 0, x - z == 0, y - z == 0)", &context);
+  poly3.removeRedundantConstraints();
+
+  // One of the three equalities can be removed.
+  EXPECT_EQ(poly3.getNumInequalities(), 0u);
+  EXPECT_EQ(poly3.getNumEqualities(), 2u);
+
+  IntegerPolyhedron poly4 =
+      parsePoly("(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q) : ("
+                "b - 1 >= 0,"
+                "-b + 500 >= 0,"
+                "-16 * d + f >= 0,"
+                "f - 1 >= 0,"
+                "-f + 998 >= 0,"
+                "16 * d - f + 15 >= 0,"
+                "-16 * e + g >= 0,"
+                "g - 1 >= 0,"
+                "-g + 998 >= 0,"
+                "16 * e - g + 15 >= 0,"
+                "h >= 0,"
+                "-h + 1 >= 0,"
+                "j - 1 >= 0,"
+                "-j + 500 >= 0,"
+                "-f + 16 * l + 15 >= 0,"
+                "f - 16 * l >= 0,"
+                "-16 * m + o >= 0,"
+                "o - 1 >= 0,"
+                "-o + 998 >= 0,"
+                "16 * m - o + 15 >= 0,"
+                "p >= 0,"
+                "-p + 1 >= 0,"
+                "-g - h + 8 * q + 8 >= 0,"
+                "-o - p + 8 * q + 8 >= 0,"
+                "o + p - 8 * q - 1 >= 0,"
+                "g + h - 8 * q - 1 >= 0,"
+                "-f + n >= 0,"
+                "f - n >= 0,"
+                "k - 10 >= 0,"
+                "-k + 10 >= 0,"
+                "i - 13 >= 0,"
+                "-i + 13 >= 0,"
+                "c - 10 >= 0,"
+                "-c + 10 >= 0,"
+                "a - 13 >= 0,"
+                "-a + 13 >= 0"
+                ")",
+                &context);
+
+  // The above is a large set of constraints without any redundant constraints,
+  // as verified by the Fourier-Motzkin based removeRedundantInequalities.
+  unsigned nIneq = poly4.getNumInequalities();
+  unsigned nEq = poly4.getNumEqualities();
+  poly4.removeRedundantInequalities();
+  ASSERT_EQ(poly4.getNumInequalities(), nIneq);
+  ASSERT_EQ(poly4.getNumEqualities(), nEq);
+  // Now we test that removeRedundantConstraints does not find any constraints
+  // to be redundant either.
+  poly4.removeRedundantConstraints();
+  EXPECT_EQ(poly4.getNumInequalities(), nIneq);
+  EXPECT_EQ(poly4.getNumEqualities(), nEq);
+
+  IntegerPolyhedron poly5 = parsePoly(
+      "(x,y) : (128 * x + 127 >= 0, -x + 7 >= 0, -128 * x + y >= 0, y >= 0)",
+      &context);
+  // 128x + 127 >= 0  implies that 128x >= 0, since x has to be an integer.
+  // (This should be caught by GCDTightenInqualities().)
+  // So -128x + y >= 0 and 128x + 127 >= 0 imply y >= 0 since we have
+  // y >= 128x >= 0.
+  poly5.removeRedundantConstraints();
+  EXPECT_EQ(poly5.getNumInequalities(), 3u);
+  SmallVector<int64_t, 8> redundantConstraint = {0, 1, 0};
+  for (unsigned i = 0; i < 3; ++i) {
+    // Ensure that the removed constraint was the redundant constraint [3].
+    EXPECT_NE(poly5.getInequality(i), ArrayRef<int64_t>(redundantConstraint));
+  }
+}
+
+TEST(IntegerPolyhedronTest, addConstantUpperBound) {
+  IntegerPolyhedron poly(2);
+  poly.addBound(IntegerPolyhedron::UB, 0, 1);
+  EXPECT_EQ(poly.atIneq(0, 0), -1);
+  EXPECT_EQ(poly.atIneq(0, 1), 0);
+  EXPECT_EQ(poly.atIneq(0, 2), 1);
+
+  poly.addBound(IntegerPolyhedron::UB, {1, 2, 3}, 1);
+  EXPECT_EQ(poly.atIneq(1, 0), -1);
+  EXPECT_EQ(poly.atIneq(1, 1), -2);
+  EXPECT_EQ(poly.atIneq(1, 2), -2);
+}
+
+TEST(IntegerPolyhedronTest, addConstantLowerBound) {
+  IntegerPolyhedron poly(2);
+  poly.addBound(IntegerPolyhedron::LB, 0, 1);
+  EXPECT_EQ(poly.atIneq(0, 0), 1);
+  EXPECT_EQ(poly.atIneq(0, 1), 0);
+  EXPECT_EQ(poly.atIneq(0, 2), -1);
+
+  poly.addBound(IntegerPolyhedron::LB, {1, 2, 3}, 1);
+  EXPECT_EQ(poly.atIneq(1, 0), 1);
+  EXPECT_EQ(poly.atIneq(1, 1), 2);
+  EXPECT_EQ(poly.atIneq(1, 2), 2);
+}
+
+/// Check if the expected division representation of local variables matches the
+/// computed representation. The expected division representation is given as
+/// a vector of expressions set in `expectedDividends` and the corressponding
+/// denominator in `expectedDenominators`. The `denominators` and `dividends`
+/// obtained through `getLocalRepr` function is verified against the
+/// `expectedDenominators` and `expectedDividends` respectively.
+static void checkDivisionRepresentation(
+    IntegerPolyhedron &poly,
+    const std::vector<SmallVector<int64_t, 8>> &expectedDividends,
+    const SmallVectorImpl<unsigned> &expectedDenominators) {
+
+  std::vector<SmallVector<int64_t, 8>> dividends;
+  SmallVector<unsigned, 4> denominators;
+
+  poly.getLocalReprs(dividends, denominators);
+
+  // Check that the `denominators` and `expectedDenominators` match.
+  EXPECT_TRUE(expectedDenominators == denominators);
+
+  // Check that the `dividends` and `expectedDividends` match. If the
+  // denominator for a division is zero, we ignore its dividend.
+  EXPECT_TRUE(dividends.size() == expectedDividends.size());
+  for (unsigned i = 0, e = dividends.size(); i < e; ++i)
+    if (denominators[i] != 0)
+      EXPECT_TRUE(expectedDividends[i] == dividends[i]);
+}
+
+TEST(IntegerPolyhedronTest, computeLocalReprSimple) {
+  IntegerPolyhedron poly(1);
+
+  poly.addLocalFloorDiv({1, 4}, 10);
+  poly.addLocalFloorDiv({1, 0, 100}, 10);
+
+  std::vector<SmallVector<int64_t, 8>> divisions = {{1, 0, 0, 4},
+                                                    {1, 0, 0, 100}};
+  SmallVector<unsigned, 8> denoms = {10, 10};
+
+  // Check if floordivs can be computed when no other inequalities exist
+  // and floor divs do not depend on each other.
+  checkDivisionRepresentation(poly, divisions, denoms);
+}
+
+TEST(IntegerPolyhedronTest, computeLocalReprConstantFloorDiv) {
+  IntegerPolyhedron poly(4);
+
+  poly.addInequality({1, 0, 3, 1, 2});
+  poly.addInequality({1, 2, -8, 1, 10});
+  poly.addEquality({1, 2, -4, 1, 10});
+
+  poly.addLocalFloorDiv({0, 0, 0, 0, 100}, 30);
+  poly.addLocalFloorDiv({0, 0, 0, 0, 0, 206}, 101);
+
+  std::vector<SmallVector<int64_t, 8>> divisions = {{0, 0, 0, 0, 0, 0, 3},
+                                                    {0, 0, 0, 0, 0, 0, 2}};
+  SmallVector<unsigned, 8> denoms = {1, 1};
+
+  // Check if floordivs with constant numerator can be computed.
+  checkDivisionRepresentation(poly, divisions, denoms);
+}
+
+TEST(IntegerPolyhedronTest, computeLocalReprRecursive) {
+  IntegerPolyhedron poly(4);
+  poly.addInequality({1, 0, 3, 1, 2});
+  poly.addInequality({1, 2, -8, 1, 10});
+  poly.addEquality({1, 2, -4, 1, 10});
+
+  poly.addLocalFloorDiv({0, -2, 7, 2, 10}, 3);
+  poly.addLocalFloorDiv({3, 0, 9, 2, 2, 10}, 5);
+  poly.addLocalFloorDiv({0, 1, -123, 2, 0, -4, 10}, 3);
+
+  poly.addInequality({1, 2, -2, 1, -5, 0, 6, 100});
+  poly.addInequality({1, 2, -8, 1, 3, 7, 0, -9});
+
+  std::vector<SmallVector<int64_t, 8>> divisions = {
+      {0, -2, 7, 2, 0, 0, 0, 10},
+      {3, 0, 9, 2, 2, 0, 0, 10},
+      {0, 1, -123, 2, 0, -4, 0, 10}};
+
+  SmallVector<unsigned, 8> denoms = {3, 5, 3};
+
+  // Check if floordivs which may depend on other floordivs can be computed.
+  checkDivisionRepresentation(poly, divisions, denoms);
+}
+
+TEST(IntegerPolyhedronTest, computeLocalReprTightUpperBound) {
+  MLIRContext context;
+
+  {
+    IntegerPolyhedron poly = parsePoly("(i) : (i mod 3 - 1 >= 0)", &context);
+
+    // The set formed by the poly is:
+    //        3q - i + 2 >= 0             <-- Division lower bound
+    //       -3q + i - 1 >= 0
+    //       -3q + i     >= 0             <-- Division upper bound
+    // We remove redundant constraints to get the set:
+    //        3q - i + 2 >= 0             <-- Division lower bound
+    //       -3q + i - 1 >= 0             <-- Tighter division upper bound
+    // thus, making the upper bound tighter.
+    poly.removeRedundantConstraints();
+
+    std::vector<SmallVector<int64_t, 8>> divisions = {{1, 0, 0}};
+    SmallVector<unsigned, 8> denoms = {3};
+
+    // Check if the divisions can be computed even with a tighter upper bound.
+    checkDivisionRepresentation(poly, divisions, denoms);
+  }
+
+  {
+    IntegerPolyhedron poly = parsePoly(
+        "(i, j, q) : (4*q - i - j + 2 >= 0, -4*q + i + j >= 0)", &context);
+    // Convert `q` to a local variable.
+    poly.convertDimToLocal(2, 3);
+
+    std::vector<SmallVector<int64_t, 8>> divisions = {{1, 1, 0, 1}};
+    SmallVector<unsigned, 8> denoms = {4};
+
+    // Check if the divisions can be computed even with a tighter upper bound.
+    checkDivisionRepresentation(poly, divisions, denoms);
+  }
+}
+
+TEST(IntegerPolyhedronTest, computeLocalReprNoRepr) {
+  MLIRContext context;
+  IntegerPolyhedron poly =
+      parsePoly("(x, q) : (x - 3 * q >= 0, -x + 3 * q + 3 >= 0)", &context);
+  // Convert q to a local variable.
+  poly.convertDimToLocal(1, 2);
+
+  std::vector<SmallVector<int64_t, 8>> divisions = {{0, 0, 0}};
+  SmallVector<unsigned, 8> denoms = {0};
+
+  // Check that no division is computed.
+  checkDivisionRepresentation(poly, divisions, denoms);
+}
+
+TEST(IntegerPolyhedronTest, simplifyLocalsTest) {
+  // (x) : (exists y: 2x + y = 1 and y = 2).
+  IntegerPolyhedron poly(1, 0, 1);
+  poly.addEquality({2, 1, -1});
+  poly.addEquality({0, 1, -2});
+
+  EXPECT_TRUE(poly.isEmpty());
+
+  // (x) : (exists y, z, w: 3x + y = 1 and 2y = z and 3y = w and z = w).
+  IntegerPolyhedron poly2(1, 0, 3);
+  poly2.addEquality({3, 1, 0, 0, -1});
+  poly2.addEquality({0, 2, -1, 0, 0});
+  poly2.addEquality({0, 3, 0, -1, 0});
+  poly2.addEquality({0, 0, 1, -1, 0});
+
+  EXPECT_TRUE(poly2.isEmpty());
+
+  // (x) : (exists y: x >= y + 1 and 2x + y = 0 and y >= -1).
+  IntegerPolyhedron poly3(1, 0, 1);
+  poly3.addInequality({1, -1, -1});
+  poly3.addInequality({0, 1, 1});
+  poly3.addEquality({2, 1, 0});
+
+  EXPECT_TRUE(poly3.isEmpty());
+}
+
+TEST(IntegerPolyhedronTest, mergeDivisionsSimple) {
+  {
+    // (x) : (exists z, y  = [x / 2] : x = 3y and x + z + 1 >= 0).
+    IntegerPolyhedron poly1(1, 0, 1);
+    poly1.addLocalFloorDiv({1, 0, 0}, 2); // y = [x / 2].
+    poly1.addEquality({1, 0, -3, 0});     // x = 3y.
+    poly1.addInequality({1, 1, 0, 1});    // x + z + 1 >= 0.
+
+    // (x) : (exists y = [x / 2], z : x = 5y).
+    IntegerPolyhedron poly2(1);
+    poly2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2].
+    poly2.addEquality({1, -5, 0});     // x = 5y.
+    poly2.appendLocalId();             // Add local id z.
+
+    poly1.mergeLocalIds(poly2);
+
+    // Local space should be same.
+    EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds());
+
+    // 1 division should be matched + 2 unmatched local ids.
+    EXPECT_EQ(poly1.getNumLocalIds(), 3u);
+    EXPECT_EQ(poly2.getNumLocalIds(), 3u);
+  }
+
+  {
+    // (x) : (exists z = [x / 5], y = [x / 2] : x = 3y).
+    IntegerPolyhedron poly1(1);
+    poly1.addLocalFloorDiv({1, 0}, 5);    // z = [x / 5].
+    poly1.addLocalFloorDiv({1, 0, 0}, 2); // y = [x / 2].
+    poly1.addEquality({1, 0, -3, 0});     // x = 3y.
+
+    // (x) : (exists y = [x / 2], z = [x / 5]: x = 5z).
+    IntegerPolyhedron poly2(1);
+    poly2.addLocalFloorDiv({1, 0}, 2);    // y = [x / 2].
+    poly2.addLocalFloorDiv({1, 0, 0}, 5); // z = [x / 5].
+    poly2.addEquality({1, 0, -5, 0});     // x = 5z.
+
+    poly1.mergeLocalIds(poly2);
+
+    // Local space should be same.
+    EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds());
+
+    // 2 divisions should be matched.
+    EXPECT_EQ(poly1.getNumLocalIds(), 2u);
+    EXPECT_EQ(poly2.getNumLocalIds(), 2u);
+  }
+
+  {
+    // Division Normalization test.
+    // (x) : (exists z, y  = [x / 2] : x = 3y and x + z + 1 >= 0).
+    IntegerPolyhedron poly1(1, 0, 1);
+    // This division would be normalized.
+    poly1.addLocalFloorDiv({3, 0, 0}, 6); // y = [3x / 6] -> [x/2].
+    poly1.addEquality({1, 0, -3, 0});     // x = 3z.
+    poly1.addInequality({1, 1, 0, 1});    // x + y + 1 >= 0.
+
+    // (x) : (exists y = [x / 2], z : x = 5y).
+    IntegerPolyhedron poly2(1);
+    poly2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2].
+    poly2.addEquality({1, -5, 0});     // x = 5y.
+    poly2.appendLocalId();             // Add local id z.
+
+    poly1.mergeLocalIds(poly2);
+
+    // Local space should be same.
+    EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds());
+
+    // One division should be matched + 2 unmatched local ids.
+    EXPECT_EQ(poly1.getNumLocalIds(), 3u);
+    EXPECT_EQ(poly2.getNumLocalIds(), 3u);
+  }
+}
+
+TEST(IntegerPolyhedronTest, mergeDivisionsNestedDivsions) {
+  {
+    // (x) : (exists y = [x / 2], z = [x + y / 3]: y + z >= x).
+    IntegerPolyhedron poly1(1);
+    poly1.addLocalFloorDiv({1, 0}, 2);    // y = [x / 2].
+    poly1.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3].
+    poly1.addInequality({-1, 1, 1, 0});   // y + z >= x.
+
+    // (x) : (exists y = [x / 2], z = [x + y / 3]: y + z <= x).
+    IntegerPolyhedron poly2(1);
+    poly2.addLocalFloorDiv({1, 0}, 2);    // y = [x / 2].
+    poly2.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3].
+    poly2.addInequality({1, -1, -1, 0});  // y + z <= x.
+
+    poly1.mergeLocalIds(poly2);
+
+    // Local space should be same.
+    EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds());
+
+    // 2 divisions should be matched.
+    EXPECT_EQ(poly1.getNumLocalIds(), 2u);
+    EXPECT_EQ(poly2.getNumLocalIds(), 2u);
+  }
+
+  {
+    // (x) : (exists y = [x / 2], z = [x + y / 3], w = [z + 1 / 5]: y + z >= x).
+    IntegerPolyhedron poly1(1);
+    poly1.addLocalFloorDiv({1, 0}, 2);       // y = [x / 2].
+    poly1.addLocalFloorDiv({1, 1, 0}, 3);    // z = [x + y / 3].
+    poly1.addLocalFloorDiv({0, 0, 1, 1}, 5); // w = [z + 1 / 5].
+    poly1.addInequality({-1, 1, 1, 0, 0});   // y + z >= x.
+
+    // (x) : (exists y = [x / 2], z = [x + y / 3], w = [z + 1 / 5]: y + z <= x).
+    IntegerPolyhedron poly2(1);
+    poly2.addLocalFloorDiv({1, 0}, 2);       // y = [x / 2].
+    poly2.addLocalFloorDiv({1, 1, 0}, 3);    // z = [x + y / 3].
+    poly2.addLocalFloorDiv({0, 0, 1, 1}, 5); // w = [z + 1 / 5].
+    poly2.addInequality({1, -1, -1, 0, 0});  // y + z <= x.
+
+    poly1.mergeLocalIds(poly2);
+
+    // Local space should be same.
+    EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds());
+
+    // 3 divisions should be matched.
+    EXPECT_EQ(poly1.getNumLocalIds(), 3u);
+    EXPECT_EQ(poly2.getNumLocalIds(), 3u);
+  }
+  {
+    // (x) : (exists y = [x / 2], z = [x + y / 3]: y + z >= x).
+    IntegerPolyhedron poly1(1);
+    poly1.addLocalFloorDiv({2, 0}, 4);    // y = [2x / 4] -> [x / 2].
+    poly1.addLocalFloorDiv({1, 1, 0}, 3); // z = [x + y / 3].
+    poly1.addInequality({-1, 1, 1, 0});   // y + z >= x.
+
+    // (x) : (exists y = [x / 2], z = [x + y / 3]: y + z <= x).
+    IntegerPolyhedron poly2(1);
+    poly2.addLocalFloorDiv({1, 0}, 2); // y = [x / 2].
+    // This division would be normalized.
+    poly2.addLocalFloorDiv({3, 3, 0}, 9); // z = [3x + 3y / 9] -> [x + y / 3].
+    poly2.addInequality({1, -1, -1, 0});  // y + z <= x.
+
+    poly1.mergeLocalIds(poly2);
+
+    // Local space should be same.
+    EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds());
+
+    // 2 divisions should be matched.
+    EXPECT_EQ(poly1.getNumLocalIds(), 2u);
+    EXPECT_EQ(poly2.getNumLocalIds(), 2u);
+  }
+}
+
+TEST(IntegerPolyhedronTest, mergeDivisionsConstants) {
+  {
+    // (x) : (exists y = [x + 1 / 3], z = [x + 2 / 3]: y + z >= x).
+    IntegerPolyhedron poly1(1);
+    poly1.addLocalFloorDiv({1, 1}, 2);    // y = [x + 1 / 2].
+    poly1.addLocalFloorDiv({1, 0, 2}, 3); // z = [x + 2 / 3].
+    poly1.addInequality({-1, 1, 1, 0});   // y + z >= x.
+
+    // (x) : (exists y = [x + 1 / 3], z = [x + 2 / 3]: y + z <= x).
+    IntegerPolyhedron poly2(1);
+    poly2.addLocalFloorDiv({1, 1}, 2);    // y = [x + 1 / 2].
+    poly2.addLocalFloorDiv({1, 0, 2}, 3); // z = [x + 2 / 3].
+    poly2.addInequality({1, -1, -1, 0});  // y + z <= x.
+
+    poly1.mergeLocalIds(poly2);
+
+    // Local space should be same.
+    EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds());
+
+    // 2 divisions should be matched.
+    EXPECT_EQ(poly1.getNumLocalIds(), 2u);
+    EXPECT_EQ(poly2.getNumLocalIds(), 2u);
+  }
+  {
+    // (x) : (exists y = [x + 1 / 3], z = [x + 2 / 3]: y + z >= x).
+    IntegerPolyhedron poly1(1);
+    poly1.addLocalFloorDiv({1, 1}, 2); // y = [x + 1 / 2].
+    // Normalization test.
+    poly1.addLocalFloorDiv({3, 0, 6}, 9); // z = [3x + 6 / 9] -> [x + 2 / 3].
+    poly1.addInequality({-1, 1, 1, 0});   // y + z >= x.
+
+    // (x) : (exists y = [x + 1 / 3], z = [x + 2 / 3]: y + z <= x).
+    IntegerPolyhedron poly2(1);
+    // Normalization test.
+    poly2.addLocalFloorDiv({2, 2}, 4);    // y = [2x + 2 / 4] -> [x + 1 / 2].
+    poly2.addLocalFloorDiv({1, 0, 2}, 3); // z = [x + 2 / 3].
+    poly2.addInequality({1, -1, -1, 0});  // y + z <= x.
+
+    poly1.mergeLocalIds(poly2);
+
+    // Local space should be same.
+    EXPECT_EQ(poly1.getNumLocalIds(), poly2.getNumLocalIds());
+
+    // 2 divisions should be matched.
+    EXPECT_EQ(poly1.getNumLocalIds(), 2u);
+    EXPECT_EQ(poly2.getNumLocalIds(), 2u);
+  }
+}
+
 } // namespace mlir


        


More information about the Mlir-commits mailing list