Exact integer emptiness checks for FlatAffineConstraints

[arjunp]

Summary

This patch adds the capability to perform exact integer emptiness checks for FlatAffineConstraints using the General Basis Reduction algorithm (GBR). Previously, only a heuristic was available for emptiness checks, which was not guaranteed to always give a conclusive result.

This patch adds a Simplex class, which can be constructed using a FlatAffineConstraints, and can find an integer sample point (if one exists) using the GBR algorithm. Additionally, it adds two classes Matrix and Fraction, which are used by Simplex.

The integer emptiness check functionality can be accessed through the new FlatAffineConstraints::isIntegerEmpty() function, which runs the existing heuristic first and, if that proves to be inconclusive, runs the GBR algorithm to produce a conclusive result.

[Kayjukh]

It looks like there are clang-format errors (https://results.llvm-merge-guard.org/amd64_debian_testing_clang-1651/summary.html). Can you format your changes using git-clang-format HEAD^ and update the diff? This would allow us to have a clean build.

[arjunp]

Fix some whitespace and clang-format.

[ftynse]

First of all, thanks for your contribution! This is an important and hard piece of polyhedral machinery.

The sheer size of the patch is huge and this will take some time to review. I only looked at the +/-style and simple issues and did not go into details of the algorithm itself. It would be helpful to know if you ran some code coverage tools on it.

Here is a list of major observations, which I did not always repeat:

• I am concerned by 32 as default number of stack elements in SmallVector, do you have evidence for this specific number?
• Please drop most of llvm:: prefixes, common classes are re-exported in the mlir:: namespace.
• In general, all top-level entities as well as class methods must have documentation.
• Consider making the API surface, both public and private, of classes smaller and using auxiliary functions with static linkage instead; same goes for multiple in-line lambdas that capture the world by-reference. These functions only need to be methods if they access class state, everything else can take that state as (const when possible) arguments and be cleaner about mutating it.
• Some additional documentation on the practical implementation details is necessary, in particular, the ownership model on the tableau and the bitwise negation in indexing.

Please bear with us.

mlir/include/mlir/Analysis/AffineStructures.h

150

Multiple common LLVM ADT structures are re-exported in mlir namespace, see the full list here https://github.com/llvm/llvm-project/blob/master/mlir/include/mlir/Support/LLVM.h. We generally avoid prefixing these with llvm:: or mlir:: for brevity.

mlir/include/mlir/Analysis/Presburger/Fraction.h

18

Prefer <cstdint> in C++ code for namespacing reasons

22

This needs a doc

34

This may still be subject to overflows. I'm fine with documenting this behavior somewhere and leaving as is, but we need to be aware of the problem. Have you considered performance implications of basic Fraction on APInt instead of int64?

mlir/include/mlir/Analysis/Presburger/Matrix.h

36

All public methods should have at least some documentation.

66

This rationale should be documented in the code. The interplay between resizes and cache behavior in case of consecutive reads (which I suppose are much more common) is non-trivial, especially if we allow reserve-ing space. I'm also slightly concerned about using 8k bytes on stack by default, we rarely use more than 8 stack elements in small vectors in MLIR.

148

Please keep the print. It allows one to print objects into streams different than llvm::errs().

mlir/include/mlir/Analysis/Presburger/Simplex.h

37–41

This comment lost me about what are unknowns, variables and constraints, and what is the relation between them, even though I know the theory behind this.

47

Could you please elaborate a little about this bit-completent indexing, e.g. that -1 gets transformed into 0 and -2 into 1 and so on, sparing you an addition. This is a clever trick, but clever tricks are very expensive to maintain

56

Nit: suffix Op is pervasively used for IR Operations, please consider a different name

96

Nit: null value -> None

113

Nit: empty llvm::Optional -> None, here and below

126

Extra points: use a bit field or a uint8_t for flags

131

Nit: -1u (or UINT_MAX if you prefer)

133

Please, provide documentation about this class. Especially the ownership model and how it is maintained.

187

Prefer mlir::LogicalResult for success/failure conditions to avoid deadling with LLVM's changing mapping of boolean values to success status

203

Nit: true _if_ value is positive

218

General design question: how much of these functions need to be private methods, i.e. have access to class members, as opposed to be just free functions, static to the compilation unit, that actual methods can call ?

237

Why is this int when pretty much everything else uses either unsigned or int64_t?

mlir/lib/Analysis/Presburger/Matrix.cpp

18 ↗

(On Diff #267824)

unsigned is used everywhere else

76 ↗

(On Diff #267824)

nit: use single quotes for one character

mlir/lib/Analysis/Presburger/Simplex.cpp

2

Drop "C++" in .cpp files, it's only necessary for headers (.h) so that tooling can tell between C and C++ that share the extension

15

Is this include necessary?

417

Please use braces symmetrically, i.e. use them if at least one if/else block needs them

534

Why does this copy the first two columns from A?

563

Please use braces symmetrically; here and below

591

Use SmallVectorImpl with references so you don't need to hardcode the size

607

Nit: consider a ternary expression here

614

This class needs a better description of what are the invariants on the snapshot stack. This method rolls back one snapshot, which happens to be be before an equality was added, but without proper documentation, there is a risk this might be broken in future.

628

Don't evaluate the exit condition on every iteration, use

for (unsigned i = 0, e = dir.size(); i < e; ++i)

697

Please consider a lengthier name for this and all derived objects, e.g. fAtI

758

Nit: asserts are for things that should never happen, no need to mention it in the description

834

This function would be cleaner if it took opt as argument and did not capture the world by-reference. It could then also go outside as a helper with static linkage

867

MLIR generally discourages recursive functions on anything that runs on pieces of IR, stack overflows in the compiler are bad. Any chance of making this iterative? Otherwise, please provide some textual argument on how fast the recursion stops.

883

Nit: print may be used outside dump, so I'd abstain from refering to "Dumping" here

[arjunp]

Add two tests of "long" sets. Fix some bugs that this revealed.

[bondhugula]

I am happy to see this submitted for review. An exact integer emptiness check had been missing and having one would be a really useful contribution. Although expensive in theory, for many of the use cases here (which are typically low / fixed dimensionality ones), things may be practical and fast enough. It just has be carefully used so as to not be detrimental to compile time for "standard" enough use cases - since there may be alternative faster approaches. It's nevertheless useful to have this technique at our disposal and have this maintained and optimized for use cases of interest, including for experimental/research purposes.

I strongly agree with all of @ftynse's comments; they already cover pretty much all of what I may have wanted to say. Just one additional comment below on the storage for Matrix.

mlir/include/mlir/Analysis/Presburger/Matrix.h

66

Please avoid using vector of vector or SmallVector of SmallVector here. (A SmallVector isn't really optimized to hold other vectors or SmallVectors as its elements.) All of these would be bad enough for performance to warrant using a 1-d storage to start with. Using a single (contiguous) buffer for all the coefficients (a single allocation holding it) will give you cache locality (spatial), fewer/no conflict misses, better prefetching, and fewer TLB misses (just one in nearly all cases). You'll also get completely inline / stack allocations for small matrices. Your API won't be really impacted. FlatAffineConstraints uses such a single buffer of coefficients via a SmallVector (you can actually copy over its relevant methods - instead of equalities and inequalities, there's just coefficients here.) Resizing isn't really a downside with this. Adding rows is straightforward since they are at the tail, and adding columns will require movement within the contiguous buffer, which is fine. FAC reserves columns to reduce movement most of the time for column addition. With that, I don't see any additional resizing or movement penalty with a 1-d buffer over vector of vectors - you have too many allocations with the latter.

[bondhugula]

mlir/include/mlir/Analysis/Presburger/Matrix.h

66

Esp. for mid-size matrices like 10x10 or 20x20, all of the allocs and deallocs from a smallvector of smallvector will be a heap killer!

[mehdi_amini]

mlir/lib/Analysis/Presburger/Simplex.cpp

628

or better with range-style: for (int i : llvm::seq<int>(dir.size())) :)

[arjunp]

Thanks for your review. Apologies for the delayed response.

• Here is a code coverage report.
• The lengths of these SmallVectors corresponds to the dimensionality. I've changed the length parameter to 8 for vectors and 64 for the 2D matrix, as it is in FlatAffineConstraints.
• The lambdas in makeProduct and reduceBasis are very specific to these functions and access a lot of the local variables present in those functions. Does it still make sense to pull these out into free functions?
• Added more documentation for the Simplex class.

mlir/include/mlir/Analysis/Presburger/Fraction.h

34

Added documentation.

56

I didn't document these overloads as it's unclear to me what the documentation could add. Is it necessary to write something for these?

mlir/include/mlir/Analysis/Presburger/Simplex.h

218

I move flippedDirection and signMatchesDirection out into free functions. reduceBasis does need access to this, but only to construct the GBRSimplex.

237

Added documentation. It's not unsigned because it's convenient to use sign comparisons to determine if an index corresponds to a constraint or a variable.

mlir/lib/Analysis/Presburger/Simplex.cpp

498

Pulling this out would require making Simplex::Unknown public or templating on the store type.

534

Nice catch. Fixed.

607

(Not applicable anymore.)

697

Changed to width. It's not strictly speaking the width but the width under some constraints, but I guess this is a good enough approximation. (and it is already referred to as "width" in computeWidthAndDuals).

834

It needs access to this for computeOptimum. Factored out into a member function computeIntegerBounds.

867

The depth of the recursion is at most the dimensionality of the set, which is not expected to be greater than say 32. This could be made iterative, albeit with some loss in readability.

[arjunp]

This comment was removed by arjunp.

[mehdi_amini]

mlir/lib/Analysis/Presburger/Simplex.cpp

867

Please document clearly the bound on the recursion. What impacts the dimensionality of the set here? When you "is not expected to be greater than say 32", how can it grow larger?

[mehdi_amini]

mlir/include/mlir/Analysis/Presburger/Simplex.h

118

s/particlar/particular/

123

Nice doc for the overall class! Thanks :)

mlir/lib/Analysis/Presburger/Simplex.cpp

474

Please avoid auto (here and everywhere applicable) when it doesn't improve readability (unsigned here right?)

[bondhugula]

mlir/lib/Analysis/Presburger/Simplex.cpp

908

Please avoid auto here and similarly below.

945

Please do not return SmallVectors this way - return value optimization here will not happen I believe leading to too many copies. Instead pass SmallVectorImpl<int64_t> by reference as an output argument. @mehdi_amini or @rriddle can confirm on the style part and RVO.

mlir/unittests/Analysis/AffineStructuresTest.cpp

35

Avoid auto here for readability.

39

Avoid auto.

115

Please terminate comments with a period when it ends.

[bondhugula]

mlir/lib/Analysis/Presburger/Simplex.cpp

74

Nit: You could leave the first three lines as the doc comment and move the ones below inline into the method as regular comments.

[arjunp]

Address comments. Made findIntegerSample iterative.

[arjunp]

mlir/include/mlir/Analysis/Presburger/Simplex.h

123

:)

mlir/lib/Analysis/Presburger/Simplex.cpp

867

It could theoretically grow arbitrarily large, although it may not happen in practice. I've made findIntegerSample iterative.

945

(I think this is no longer applicable since findIntegerSample is iterative now)

[bondhugula]

Thanks very much for all the changes. A few more requests - a couple on test cases and other minor ones.

mlir/unittests/Analysis/AffineStructuresTest.cpp

179

I think you are missing a test case involving equalities that is integer empty but not rational empty - a flat one along the lines of the test case added with purely inequalities. It could just be a 1-d line segment in a 2-d space that passes in between and thus misses all integer points - it could be potentially slanted as well.

mlir/unittests/Analysis/Presburger/SimplexTest.cpp

17

It may be good to add a line on what the test is testing.

87–88

Missing a test case here that involves equalities?

[bondhugula]

mlir/lib/Analysis/Presburger/Simplex.cpp

455–459

A short comment here please.

576

Comment here please.

mlir/unittests/Analysis/Presburger/SimplexTest.cpp

23

Avoid auto here please.

[bondhugula]

Not for this patch, but if you want to just quickly try and test it on more cases, you could just replace mlir::boundCheckLoadOrStoreOp's isEmpty() with isIntegerEmpty() and see if -test-memref-bound-check test cases pass. The -test-memref-bound-check pass tests for out of bound memref accesses and actually needs an integer exact method (whenever it determines non-empty) to be accurate. This way you would also have a speed comparison, and it would also be a practical application of this technique for an analysis pass.

tests/Transforms/memref-bound-check.mlir

[arjunp]

Addressed comments.

[arjunp]

I've added some more tests and documentation as requested.

[bondhugula]

mlir/unittests/Analysis/AffineStructuresTest.cpp

187

I'm sorry this would be caught by the GCD test itself. Instead, could you please replace the equality with something to the effect:

3x - 3y >= -2
3x - 3y  <= -1

This should be empty as well but would exercise your method given how isIntegerEmpty() is structured now.

[bondhugula]

mlir/include/mlir/Analysis/Presburger/Fraction.h

66

Are you assuming f.den is always positive? If both num and den are both negative, you'd want to add one.

[arjunp]

Added the requested test.

[arjunp]

Fix a comment in Fraction.

[arjunp]

mlir/include/mlir/Analysis/Presburger/Fraction.h

66

Yes, den is always positive. There was an incomplete comment for the member declaration that I've fixed. (This is also documented in the top-level class documentation)

mlir/unittests/Analysis/AffineStructuresTest.cpp

187

Since checkSample directly calls findIntegerSample() rather than isIntegerEmpty(), this shouldn't be an issue. Nonetheless, I've also added the test you requested.

[ftynse]

Thanks for providing more detailed documentation!

I did not extensively analyze the algorithms. At a glance, the implementation looks plausible. Again, given the complexity, I would encourage you to get 100% test coverage on the core functions. Another way to stress-test this would be to generate a bunch of arbitrary sets and run this + gcd emptiness checks and analyze the differences manually to make sure they only happen when gcd is not supposed to be correct.

I don't have any further comments.

@andydavis1 had a bug which might have been related to integer exactness, but I don't know if it is still active.

mlir/include/mlir/Analysis/Presburger/Simplex.h

38

Nit: SImplex -> simplex

46

Nit: 2x - 3y >= 0 (currently uses z)

[bondhugula]

mlir/include/mlir/Analysis/Presburger/Fraction.h

66

If den is positive, this can be simplified to avoid an extra check.

f.num % f.den > 0 ? f.num / f.den + 1 : f.num / f.den;

(one comparison instead of two comparisons)

Also, if there is a doc comment here, please mention that den is positive. Actually, you don't even need this logic; mlir::ceilDiv does exactly this.

[bondhugula]

In D80860#2106414, @ftynse wrote:

Thanks for providing more detailed documentation!

I did not extensively analyze the algorithms. At a glance, the implementation looks plausible. Again, given the complexity, I would encourage you to get 100% test coverage

+1 I don't think it's practical to have a review checking each of the steps of the algorithm. I would just rely on (1) good test case coverage, and (2) sufficient documentation/comments for now. Other implementation aspects can be brought up and revised later when there are other users with a specialized interest in digging deeper.

[andydavis1]

Thanks for your contribution!

It would be great to have tests on integer sets which resulting in an over approximation using the existing FM algorithm, but which returned integer exact results using GBR from this change.

Also wondering what the plan is to check for integer overflow of coefficient values for some of these algorithms?

[arjunp]

Improved test coverage. There is still one uncovered line; it seems quite difficult to find a case that covers this. I ran hundreds of tests with randomly generated sets intersected with a 1000 x 1000 x 1000 cube. Whenever GBR and the FM/GCD-based algorithm disagreed, I checked against a brute-force that checks every point in this cube. All these tests passed, as well as tests from some other types of generated cases.

Added some tests that GBR detects as empty but the existing FM/GCD-based algorithm does not.

Regarding overflow issues, it does not seem to be something specific to this patch. It looks to me like there already exists other functionality in FlatAffineConstraints, like FM, that performs operations without complete overflow checks. I feel that all these need to be addressed together by, say, introducing arbitrary-precision arithmetic, which is beyond the scope of this patch.

Thanks again for all your help with improving this patch!

[bondhugula]

mlir/include/mlir/Analysis/Presburger/Fraction.h

39

This has to be a /// comment.

[bondhugula]

In D80860#2115848, @arjunp wrote:

Added some tests that GBR detects as empty but the existing FM/GCD-based algorithm does not.

It'll be useful if you marked these tests with a comment to that effect. I might have missed it if you already did.

[bondhugula]

A bunch of mostly minor comments.

mlir/lib/Analysis/Presburger/Simplex.cpp

707

Can you rephrase this? <dir, x-y> is unclear.

711

Please use coeffs.reserve since you know big it's going to be.

712–714

Why emplace_back instead of push_back for a POD?

719–723

Triple comments please.

866–867

Reflow - use the width?

875

Please make this a complete sentence.

961

You can use llvm::to_vector.

1011

llvm::to_vector

mlir/unittests/Analysis/AffineStructuresTest.cpp

30

Nit: hasValue -> 'hasValue'

56

SmallVector -> ArrayRef

57

Likewise.

68

Typo: least

72–73

Please pass ArrayRef instead of SmallVector for inputs. (If these were output args or inouts, SmallVectorImpl is to be used.)

mlir/unittests/Analysis/Presburger/SimplexTest.cpp

54

SmallVector 4?

60

Likewise.

[arjunp]

Addressed comments.

[arjunp]

In D80860#2118522, @bondhugula wrote:

In D80860#2115848, @arjunp wrote:

Added some tests that GBR detects as empty but the existing FM/GCD-based algorithm does not.

It'll be useful if you marked these tests with a comment to that effect. I might have missed it if you already did.

This is marked at AffineStructuresTest.cpp:247.

mlir/lib/Analysis/Presburger/Simplex.cpp

707

I've replaced it with the more explicit dotProduct everywhere now. (this was originally explained in the top-level class comment)

mlir/unittests/Analysis/AffineStructuresTest.cpp

65

Is there a better way? This is needed since we moved to ArrayRefs for arguments but in checkPermutationsSample below we need to pass SmallVectors of SmallVectors to makeFACFromConstraints.

[bondhugula]

mlir/unittests/Analysis/AffineStructuresTest.cpp

65

We don't need to switch to ArrayRef of ArrayRef if you have SmallVector of SmallVector, but only to ArrayRef of SmallVector.

[bondhugula]

mlir/unittests/Analysis/AffineStructuresTest.cpp

65

I actually meant ArrayRef<SmallVector,...> as the function argument instead of SmallVector<SmallVector<.... Make sure you don't run into underlying storage lifetime issues. ArrayRef<SmallVector..> should be fine.

[arjunp]

Switched from ArrayRef of ArrayRef to ArrayRef of SmallVectors.

[ftynse]

@Kayjukh your review is blocking this revision, have your concerns been addressed?

[Kayjukh]

Yes, sorry about blocking the revision. Note that there is still a clang-format issue in one of the files.

mlir/unittests/Analysis/AffineStructuresTest.cpp

89

Please fix formatting

[arjunp]

clang-format.

[arjunp]

I don't have commit access yet, so I think someone else will need to land this.

[bondhugula]

In D80860#2127981, @arjunp wrote:

I don't have commit access yet, so I think someone else will need to land this.

I've committed this.

[mehdi_amini]

mlir/lib/Analysis/Presburger/Simplex.cpp

132

What guarantees that gcd can't be 0 and produce an invalid division here?
Should we guard this loop?

(flagged by Coverity)

[mehdi_amini]

mlir/lib/Analysis/Presburger/Simplex.cpp

326

If retRow is true, then retConst is used uninitialized right? Can you look into this?
(flagged by Coverity)

326

(Same with retElem)

[arjunp]

mlir/lib/Analysis/Presburger/Simplex.cpp

132

One of the elements being GCD'd is a "denominator" value that is always non-zero. I will add a comment noting this.

326

If retRow has a value then lines 319-322 were executed in the iteration that that values was set in, so the values of retElem and retConst were also set in that iteration already.

[grosser]

Hi @arjunp, I am wondering if there are ways to prevent future coverity warnings on this code?

mlir/lib/Analysis/Presburger/Simplex.cpp

132

Given that these are coverity errors, it might be useful to document this with an assert to ensure that coverity does not warn about this anymore.

326

Again, as these are coverity errors the question is how we can avoid coverity to flag these issues. I am almost tempted to just suggest to set retElement and retConst to 0 and add a comment stating that this is done to avoid the coverity warning.

[arjunp]

mlir/lib/Analysis/Presburger/Simplex.cpp

132

Done: https://reviews.llvm.org/rG0f925339e52e361946a1e894cfa6dbe494b39290

326

Done: https://reviews.llvm.org/rGef8351598ef33b906fd1962420c8b464812cea85

mlir/include/mlir/Analysis/AffineStructures.h

public:
   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 has no integer solution, and False
+  // if it does have a solution.
+  bool isIntegerEmpty() const;
+
+  // Find a sample point satisfying the constraints. This uses a branch and
+  // bound algorithm with general basis reduction.
+  //
+  // Returns such a point if one exists, or an empty llvm::Optional otherwise.
+  llvm::Optional<std::vector<int64_t>> findIntegerSample() const;
+
   // Clones this object.
   std::unique_ptr<FlatAffineConstraints> clone() const;
 
   /// Returns the value at the specified equality row and column.
   inline int64_t atEq(unsigned i, unsigned j) const {
     return equalities[i * numReservedCols + j];
   }
   inline int64_t &atEq(unsigned i, unsigned j) {

mlir/include/mlir/Analysis/Presburger/Fraction.h

+//===- Simplex.h - MLIR Fraction Class --------------------------*- 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
+//
+//===----------------------------------------------------------------------===//
+//
+// This is a simple class to represent fractions. It supports multiplication,
+// comparison, floor, and ceiling operations.
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef MLIR_ANALYSIS_PRESBURGER_FRACTION_H
+#define MLIR_ANALYSIS_PRESBURGER_FRACTION_H
+
+namespace mlir {
+
+template <typename T>
+struct Fraction {
+  Fraction() : num(0), den(1) {}
+  Fraction(T oNum, T oDen) : num(oNum), den(oDen) {
+    if (den < 0) {
+      num = -num;
+      den = -den;
+    }
+  }
+  T num, den;
+};
+
+template <typename T>
+int compare(const Fraction<T> &x, const Fraction<T> &y) {
+  auto diff = x.num * y.den - y.num * x.den;
+  if (diff > 0)
+    return +1;
+  if (diff < 0)
+    return -1;
+  return 0;
+}
+
+// Division rounds towards zero, so if the result is already non-negative
+// then num/den is the floor. If the result is negative and the division
+// leaves a remainder then we need to subtract one to get the floor.
+template <typename T>
+T floor(Fraction<T> f) {
+  return f.num / f.den - (f.num < 0 && f.num % f.den != 0);
+}
+
+// Division rounds towards zero, so if the result is already non-positive
+// then num/den is the ceiling. If the result is positive and the division
+// leaves a remainder then we need to add one to get the ceiling.
+template <typename T>
+T ceil(Fraction<T> f) {
+  return f.num / f.den + (f.num > 0 && f.num % f.den != 0);
+}
+
+template <typename T>
+int compare(const T &x, const Fraction<T> &y) {
+  return sign(x * y.den - y.num);
+}
+
+template <typename T>
+int compare(const Fraction<T> &x, const T &y) {
+  return -compare(y, x);
+}
+
+template <typename T>
+Fraction<T> operator-(const Fraction<T> &x) {
+  return Fraction<T>(-x.num, x.den);
+}
+
+template <typename T, typename U>
+bool operator<(const T &x, const Fraction<U> &y) {
+  return compare(x, y) < 0;
+}
+
+template <typename T, typename U>
+bool operator<=(const T &x, const Fraction<U> &y) {
+  return compare(x, y) <= 0;
+}
+
+template <typename T, typename U>
+bool operator==(const T &x, const Fraction<U> &y) {
+  return compare(x, y) == 0;
+}
+template <typename T, typename U>
+bool operator>(const T &x, const Fraction<U> &y) {
+  return compare(x, y) > 0;
+}
+
+template <typename T, typename U>
+bool operator>=(const T &x, const Fraction<U> &y) {
+  return compare(x, y) >= 0;
+}
+
+template <typename T>
+Fraction<T> operator*(const Fraction<T> &x, const Fraction<T> &y) {
+  return Fraction<T>(x.num * y.num, x.den * y.den);
+}
+} // namespace mlir
+#endif // MLIR_ANALYSIS_PRESBURGER_FRACTION_H

mlir/include/mlir/Analysis/Presburger/Matrix.h

+//===- Matrix.h - MLIR Matrix Class -----------------------------*- 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
+//
+//===----------------------------------------------------------------------===//
+//
+// This a simplex 2D matrix class that supports reading, writing, resizing,
+// and swapping rows, and swapping columns.
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef MLIR_ANALYSIS_PRESBURGER_MATRIX_H
+#define MLIR_ANALYSIS_PRESBURGER_MATRIX_H
+
+#include "llvm/Support/raw_ostream.h"
+
+#include <cassert>
+#include <vector>
+
+namespace mlir {
+
+/// This is a simple class to represent a resizable matrix.
+///
+/// The data is stored in the form of a vector of vectors.
+template <typename INT>
+class Matrix {
+public:
+  Matrix() = delete;
+
+  /// Construct a matrix with the specified number of rows and columns.
+  /// Initially, the values are default initialized.
+  Matrix(unsigned rows, unsigned columns);
+
+  static Matrix getIdentityMatrix(unsigned dimension);
+
+  /// Access the element at the specified row and column.
+  INT &operator()(unsigned row, unsigned column);
+  INT operator()(unsigned row, unsigned column) const;
+
+  /// Swap the given columns.
+  void swapColumns(unsigned column, unsigned otherColumn);
+
+  /// Swap the given rows.
+  void swapRows(unsigned row, unsigned otherRow);
+
+  unsigned getNumRows() const;
+
+  unsigned getNumColumns() const;
+
+  const std::vector<INT> getRow(unsigned row) const;
+
+  void addToRow(unsigned sourceRow, unsigned targetRow, INT scale);
+
+  /// Resize the matrix to the specified dimensions. If a dimension is smaller,
+  /// the values are truncated; if it is bigger, the new values are default
+  /// initialized.
+  void resize(unsigned newNRows, unsigned newNColumns);
+
+  void print(llvm::raw_ostream &os) const;
+  void dump() const;
+
+private:
+  unsigned nColumns;
+  std::vector<std::vector<INT>> data;
+};
+
+template <typename INT>
+Matrix<INT>::Matrix(unsigned rows, unsigned columns)
+    : nColumns(columns), data(rows, std::vector<INT>(columns)) {}
+
+template <typename INT>
+Matrix<INT> Matrix<INT>::getIdentityMatrix(unsigned dimension) {
+  Matrix<INT> matrix(dimension, dimension);
+  for (size_t i = 0; i < dimension; ++i)
+    matrix(i, i) = 1;
+  return matrix;
+}
+
+template <typename INT>
+INT &Matrix<INT>::operator()(unsigned row, unsigned column) {
+  assert(row < getNumRows() && "Row outside of range");
+  assert(column < getNumColumns() && "Column outside of range");
+  return data[row][column];
+}
+
+template <typename INT>
+INT Matrix<INT>::operator()(unsigned row, unsigned column) const {
+  assert(row < getNumRows() && "Row outside of range");
+  assert(column < getNumColumns() && "Column outside of range");
+  return data[row][column];
+}
+
+template <typename INT>
+unsigned Matrix<INT>::getNumRows() const {
+  return data.size();
+}
+
+template <typename INT>
+unsigned Matrix<INT>::getNumColumns() const {
+  return nColumns;
+}
+
+template <typename INT>
+void Matrix<INT>::resize(unsigned newNRows, unsigned newNColumns) {
+  nColumns = newNColumns;
+  for (auto &row : data)
+    row.resize(nColumns);
+  data.resize(newNRows, std::vector<INT>(nColumns));
+}
+
+template <typename INT>
+void Matrix<INT>::swapRows(unsigned row, unsigned otherRow) {
+  assert((row < getNumRows() && otherRow < getNumRows()) &&
+         "Given row out of bounds");
+  if (row == otherRow)
+    return;
+  swap(data[row], data[otherRow]);
+}
+
+template <typename INT>
+void Matrix<INT>::swapColumns(unsigned column, unsigned otherColumn) {
+  assert((column < getNumColumns() && otherColumn < getNumColumns()) &&
+         "Given column out of bounds");
+  if (column == otherColumn)
+    return;
+  for (auto &row : data)
+    std::swap(row[column], row[otherColumn]);
+}
+
+template <typename INT>
+const std::vector<INT> Matrix<INT>::getRow(unsigned row) const {
+  return data[row];
+}
+
+template <typename INT>
+void Matrix<INT>::addToRow(unsigned sourceRow, unsigned targetRow, INT scale) {
+  if (scale == 0)
+    return;
+  for (unsigned col = 0; col < getNumColumns(); ++col)
+    data[targetRow][col] += scale * data[sourceRow][col];
+  return;
+}
+
+template <typename INT>
+void Matrix<INT>::print(llvm::raw_ostream &os) const {
+  os << "Dumping matrix, rows = " << getNumRows()
+     << ", columns: " << getNumColumns() << '\n';
+  os << "r/c  ";
+  for (unsigned column = 0; column < getNumColumns(); ++column)
+    os << "| " << column << " ";
+  os << '\n';
+  os << std::string(5 + getNumColumns() * 5, '-') << '\n';
+  for (unsigned row = 0; row < getNumRows(); ++row) {
+    os << row << " | ";
+    for (unsigned column = 0; column < getNumColumns(); ++column)
+      os << data[row][column] << " ";
+    os << '\n';
+  }
+}
+
+template <typename INT>
+void Matrix<INT>::dump() const {
+  print(llvm::errs());
+}
+
+} // namespace mlir
+#endif // MLIR_ANALYSIS_PRESBURGER_MATRIX_H

mlir/include/mlir/Analysis/Presburger/Simplex.h

+//===- Simplex.h - MLIR Simplex Class ---------------------------*- 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
+//
+//===----------------------------------------------------------------------===//
+//
+// This class can perform analysis on FlatAffineConstraints. In particular,
+// it can be used to perform emptiness checks.
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef MLIR_ANALYSIS_PRESBURGER_SIMPLEX_H
+#define MLIR_ANALYSIS_PRESBURGER_SIMPLEX_H
+
+#include "mlir/Analysis/AffineStructures.h"
+#include "mlir/Analysis/Presburger/Fraction.h"
+#include "mlir/Analysis/Presburger/Matrix.h"
+#include "llvm/ADT/ArrayRef.h"
+#include "llvm/ADT/Optional.h"
+#include "llvm/ADT/SmallVector.h"
+#include "llvm/Support/raw_ostream.h"
+
+namespace mlir {
+
+class GBRSimplex;
+
+/// This class implements the Simplex algorithm. It supports adding affine
+/// equalities and inequalities, and can perform emptiness checks, i.e., it can
+/// find a solution to the set of constraints if one exists, or say that the
+/// set is empty if no solution exists.
+///
+/// An unknown is addressed by its index. If the index i is non-negative, then
+/// the ith variable is the Unknown being addressed. If the index is negative,
+/// then a constraint is being addressed, having index ~i.
+///
+/// The unknown corresponding to each row r (resp. column c) has index rowVar[r]
+/// (resp. colVar[c]). If at some point it is detected that the set of
+/// constraints are mutually contradictory and have no solution, then empty will
+/// be set to true.
+class Simplex {
+public:
+  enum class Direction { UP, DOWN };
+  enum class UndoOp {
+    DEALLOCATE,
+    UNMARK_EMPTY,
+  };
+
+  Simplex() = delete;
+  explicit Simplex(unsigned nVar);
+  explicit Simplex(const FlatAffineConstraints &constraints);
+  unsigned getNumRows() const;
+  unsigned getNumColumns() const;
+
+  /// \returns True if the tableau is empty (has conflicting constraints),
+  /// False otherwise.
+  bool isEmpty() const;
+
+  /// Add an inequality to the tableau. If coeffs is c_0, c_1, ... c_n, where n
+  /// is the current number of variables, then the corresponding inequality is
+  /// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} >= 0.
+  void addInequality(ArrayRef<int64_t> coeffs);
+
+  /// \returns the number of variables in the tableau.
+  unsigned numberVariables() const;
+
+  /// \returns the number of constraints in the tableau.
+  unsigned numberConstraints() const;
+
+  /// Add an equality to the tableau. If coeffs is c_0, c_1, ... c_n, where n
+  /// is the current number of variables, then the corresponding equality is
+  /// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} == 0.
+  void addEquality(ArrayRef<int64_t> coeffs);
+
+  /// Mark the tableau as being empty.
+  void markEmpty();
+
+  /// Get a snapshot of the current state. This is used for rolling back.
+  unsigned getSnapshot() const;
+
+  /// Rollback to a snapshot. This invalidates all later snapshots.
+  void rollback(unsigned snapshot);
+
+  /// Compute the maximum or minimum value of the given row, depending on
+  /// \p direction.
+  ///
+  /// \returns a {num, den} pair denoting the optimum, or a null value if no
+  /// optimum exists, i.e., if the expression is unbounded in this direction.
+  llvm::Optional<Fraction<int64_t>> computeRowOptimum(Direction direction,
+                                                      unsigned row);
+
+  /// Compute the maximum or minimum value of the given expression, depending on
+  /// \p direction.
+  ///
+  /// \returns a {num, den} pair denoting the optimum, or a null value if no
+  /// optimum exists, i.e., if the expression is unbounded in this direction.
+  llvm::Optional<Fraction<int64_t>> computeOptimum(Direction direction,
+                                                   ArrayRef<int64_t> coeffs);
+
+  /// Make a tableau to represent a pair of points in the given tableaus, one in
+  /// tableau \p A and one in \p B.
+  static Simplex makeProduct(const Simplex &a, const Simplex &b);
+
+  /// \returns the current sample point if it is integral. Otherwise, returns an
+  /// empty llvm::Optional.
+  llvm::Optional<std::vector<int64_t>> getSamplePointIfIntegral() const;
+
+  /// \returns an integral sample point if one exists, or an empty
+  /// llvm::Optional otherwise.
+  llvm::Optional<std::vector<int64_t>> findIntegerSample();
+
+  // Print the tableau's internal state.
+  void print(llvm::raw_ostream &os) const;
+  void dump() const;
+
+private:
+  friend class GBRSimplex;
+
+  struct Unknown {
+    Unknown(bool oOwnsRow, bool oRestricted, unsigned oPos)
+        : ownsRow(oOwnsRow), restricted(oRestricted), pos(oPos) {}
+    Unknown() : Unknown(false, false, -1) {}
+    bool ownsRow;
+    bool restricted;
+    unsigned pos;
+  };
+
+  // Dump the internal state of the unknown.
+  void printUnknown(llvm::raw_ostream &os, const Unknown &u) const;
+
+  /// Find a pivot to change the sample value of \p row in the specified
+  /// direction. The returned pivot row will be \p row if and only
+  /// if the unknown is unbounded in the specified direction.
+  ///
+  /// \returns a [row, col] pair denoting a pivot, or an empty llvm::Optional if
+  /// no valid pivot exists.
+  llvm::Optional<std::pair<unsigned, unsigned>>
+  findPivot(int row, Direction direction) const;
+
+  /// Swap the row with the column in the tableau's data structures but not the
+  /// tableau itself. This is used by pivot.
+  void swapRowWithCol(unsigned row, unsigned col);
+
+  /// Pivot the row with the column.
+  void pivot(unsigned row, unsigned col);
+  void pivot(const std::pair<unsigned, unsigned> &p);
+
+  /// \returns the unknown associated with \p index.
+  const Unknown &unknownFromIndex(int index) const;
+  /// \returns the unknown associated with \p col.
+  const Unknown &unknownFromColumn(unsigned col) const;
+  /// \returns the unknown associated with \p row.
+  const Unknown &unknownFromRow(unsigned row) const;
+  /// \returns the unknown associated with \p index.
+  Unknown &unknownFromIndex(int index);
+  /// \returns the unknown associated with \p col.
+  Unknown &unknownFromColumn(unsigned col);
+  /// \returns the unknown associated with \p row.
+  Unknown &unknownFromRow(unsigned row);
+
+  /// Add a new row to the tableau and the associated data structures.
+  unsigned addRow(ArrayRef<int64_t> coeffs);
+
+  /// Normalize the given row by removing common factors between the numerator
+  /// and the denominator.
+  void normalizeRow(unsigned row);
+
+  /// Swap the two rows in the tableau and associated data structures.
+  void swapRows(unsigned i, unsigned j);
+
+  /// Restore the unknown to a non-negative sample value.
+  ///
+  /// \returns True if the unknown was successfully restored to a non-negative
+  /// sample value, False otherwise.
+  bool restoreRow(Unknown &u);
+
+  void undoOp(UndoOp op, llvm::Optional<int> index);
+
+  /// Find a row that can be used to pivot the column in the specified
+  /// direction. If \p skipRow is not null, then this row is excluded
+  /// from consideration. The returned pivot will maintain all constraints
+  /// except the column itself and \p skipRow, if it is set. (if these unknowns
+  /// are restricted).
+  ///
+  /// \returns the row to pivot to, or an empty llvm::Optional if the column
+  /// is unbounded in the specified direction.
+  llvm::Optional<unsigned> findPivotRow(llvm::Optional<unsigned> skipRow,
+                                        Direction direction,
+                                        unsigned col) const;
+
+  /// \returns True \p value is positive and direction is Direction::UP, or if
+  /// \p value is negative and direction is Direction::DOWN. Returns False
+  /// otherwise.
+  bool signMatchesDirection(int64_t value, Direction direction) const;
+
+  /// \returns Direction::UP if \p direction is Direction::DOWN and vice versa.
+  Direction flippedDirection(Direction direction) const;
+
+  /// Searches for an integer sample point recursively using a branch and bound
+  /// algorithm and general basis reduction.
+  llvm::Optional<std::vector<int64_t>>
+  findIntegerSampleRecursively(Matrix<int64_t> &basis, unsigned level);
+
+  /// Reduce the given basis, starting at the specified level, using general
+  /// basis reduction.
+  void reduceBasis(Matrix<int64_t> &basis, unsigned level);
+
+  /// The number of rows in the tableau.
+  unsigned nRow;
+
+  /// The number of columns in the tableau, including the common denominator
+  /// and the constant column.
+  unsigned nCol;
+
+  /// The matrix representing the tableau.
+  Matrix<int64_t> tableau;
+
+  /// True if the tableau has been detected to be empty, False otherwise.
+  bool empty;
+
+  /// Holds a log of operations, used for rolling back to a previoous state.
+  std::vector<std::pair<UndoOp, llvm::Optional<int>>> undoLog;
+
+  /// These hold the indexes of the unknown at a given row or column position.
+  std::vector<int> rowVar, colVar;
+
+  /// These hold information about each unknown.
+  std::vector<Unknown> con, var;
+};
+} // namespace mlir
+#endif // MLIR_ANALYSIS_PRESBURGER_SIMPLEX_H

mlir/lib/Analysis/AffineStructures.cpp

 //===- AffineStructures.cpp - MLIR Affine Structures Class-----------------===//
 //
 // 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
 //
 //===----------------------------------------------------------------------===//
 //
 // Structures for affine/polyhedral analysis of affine dialect ops.
 //
 //===----------------------------------------------------------------------===//
 
 #include "mlir/Analysis/AffineStructures.h"
+#include "mlir/Analysis/Presburger/Simplex.h"
 #include "mlir/Dialect/Affine/IR/AffineOps.h"
 #include "mlir/Dialect/Affine/IR/AffineValueMap.h"
 #include "mlir/Dialect/StandardOps/IR/Ops.h"
 #include "mlir/IR/AffineExprVisitor.h"
 #include "mlir/IR/IntegerSet.h"
 #include "mlir/Support/LLVM.h"
 #include "mlir/Support/MathExtras.h"
 #include "llvm/ADT/SmallPtrSet.h"
for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
     int64_t v = std::abs(atEq(i, numCols - 1));
     if (gcd > 0 && (v % gcd != 0)) {
       return true;
     }
   }
   return false;
 }
 
+// First, try the GCD test heuristic. If that doesn't find the set empty,
+// use the complete emptiness check provided by findIntegerSample, which gives
+// a definitive answer.
+bool FlatAffineConstraints::isIntegerEmpty() const {
+  if (isEmptyByGCDTest())
+    return true;
+  return !findIntegerSample().hasValue();
+}
+
+llvm::Optional<std::vector<int64_t>>
+FlatAffineConstraints::findIntegerSample() const {
+  Simplex simplex(*this);
+  return simplex.findIntegerSample();
+}
+
 /// 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.

mlir/lib/Analysis/CMakeLists.txt

add_mlir_library(MLIRAnalysis
   DEPENDS
   mlir-headers
 
   LINK_LIBS PUBLIC
   MLIRAffineOps
   MLIRCallInterfaces
   MLIRControlFlowInterfaces
   MLIRInferTypeOpInterface
+  MLIRPresburger
   MLIRSCF
   )
 
 add_mlir_library(MLIRLoopAnalysis
   AffineAnalysis.cpp
   AffineStructures.cpp
   LoopAnalysis.cpp
   NestedMatcher.cpp
   Utils.cpp
 
   ADDITIONAL_HEADER_DIRS
   ${MLIR_MAIN_INCLUDE_DIR}/mlir/Analysis
 
   DEPENDS
   mlir-headers
 
   LINK_LIBS PUBLIC
   MLIRAffineOps
   MLIRCallInterfaces
   MLIRControlFlowInterfaces
   MLIRInferTypeOpInterface
+  MLIRPresburger
   MLIRSCF
   )
+  
+add_subdirectory(Presburger)

mlir/lib/Analysis/Presburger/CMakeLists.txt

+add_mlir_library(MLIRPresburger
+  Simplex.cpp
+  )
 No newline at end of file

mlir/lib/Analysis/Presburger/Simplex.cpp

+//===- Simplex.cpp - MLIR Simplex Class -------------------------*- 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/Presburger/Simplex.h"
+#include "mlir/Analysis/Presburger/Matrix.h"
+#include "mlir/Support/MathExtras.h"
+
+#include <algorithm>
+#include <cstdio>
+
+namespace mlir {
+const int nullIndex = std::numeric_limits<int>::max();
+
+// Construct a Simplex object with `nVar` variables.
+Simplex::Simplex(unsigned nVar)
+    : nRow(0), nCol(2), tableau(0, 2 + nVar), empty(false) {
+  colVar.push_back(nullIndex);
+  colVar.push_back(nullIndex);
+  for (unsigned i = 0; i < nVar; ++i) {
+    var.push_back(
+        Unknown(/*ownsRow=*/false, /*restricted=*/false, /*pos=*/nCol));
+    colVar.push_back(i);
+    nCol++;
+  }
+}
+
+Simplex::Simplex(const FlatAffineConstraints &constraints)
+    : Simplex(constraints.getNumIds()) {
+  for (unsigned i = 0; i < constraints.getNumInequalities(); ++i)
+    addInequality(constraints.getInequality(i));
+  for (unsigned i = 0; i < constraints.getNumEqualities(); ++i)
+    addEquality(constraints.getEquality(i));
+}
+
+const Simplex::Unknown &Simplex::unknownFromIndex(int index) const {
+  assert(index != nullIndex && "nullIndex passed to unknownFromIndex");
+  return index >= 0 ? var[index] : con[~index];
+}
+
+const Simplex::Unknown &Simplex::unknownFromColumn(unsigned col) const {
+  assert(col < nCol && "Invalid column");
+  return unknownFromIndex(colVar[col]);
+}
+
+const Simplex::Unknown &Simplex::unknownFromRow(unsigned row) const {
+  assert(row < nRow && "Invalid row");
+  return unknownFromIndex(rowVar[row]);
+}
+
+Simplex::Unknown &Simplex::unknownFromIndex(int index) {
+  assert(index != nullIndex && "nullIndex passed to unknownFromIndex");
+  return index >= 0 ? var[index] : con[~index];
+}
+
+Simplex::Unknown &Simplex::unknownFromColumn(unsigned col) {
+  assert(col < nCol && "Invalid column");
+  return unknownFromIndex(colVar[col]);
+}
+
+Simplex::Unknown &Simplex::unknownFromRow(unsigned row) {
+  assert(row < nRow && "Invalid row");
+  return unknownFromIndex(rowVar[row]);
+}
+
+// Add a new row to the tableau corresponding to the given constant term and
+// list of coefficients. The coefficients are specified as a vector of
+// (variable index, coefficient) pairs.
+//
+// If the tableau is not big enough to accomodate the extra row, we extend it.
+//
+// We then process each given variable coefficient.
+// If a variable is in column position at column col, then we just add the
+// coefficient for that variable (scaled by the common row denominator) to
+// the corresponding entry in the new row.
+// If the variable is in row position, we need to add that row to the new row,
+// scaled by the coefficient for the variable, accounting for the two rows
+// potentially having different denominators. The new denominator is the
+// lcm of the two.
+unsigned Simplex::addRow(ArrayRef<int64_t> coeffs) {
+  assert(coeffs.size() == 1 + var.size() &&
+         "Incorrect number of coefficients!");
+
+  nRow++;
+  if (nRow >= tableau.getNumRows())
+    tableau.resize(nRow, tableau.getNumColumns());
+  rowVar.push_back(~con.size());
+  con.emplace_back(true, false, nRow - 1);
+
+  tableau(nRow - 1, 0) = 1;
+  tableau(nRow - 1, 1) = coeffs.back();
+  for (unsigned col = 2; col < nCol; ++col)
+    tableau(nRow - 1, col) = 0;
+
+  for (unsigned i = 0; i < var.size(); ++i) {
+    unsigned pos = var[i].pos;
+    if (coeffs[i] == 0)
+      continue;
+
+    if (!var[i].ownsRow) {
+      tableau(nRow - 1, pos) += coeffs[i] * tableau(nRow - 1, 0);
+      continue;
+    }
+
+    int64_t lcm = mlir::lcm(tableau(nRow - 1, 0), tableau(pos, 0));
+    int64_t nRowCoeff = lcm / tableau(nRow - 1, 0);
+    int64_t idxRowCoeff = coeffs[i] * (lcm / tableau(pos, 0));
+    tableau(nRow - 1, 0) = lcm;
+    for (unsigned col = 1; col < nCol; ++col)
+      tableau(nRow - 1, col) =
+          nRowCoeff * tableau(nRow - 1, col) + idxRowCoeff * tableau(pos, col);
+  }
+
+  normalizeRow(nRow - 1);
+  // Push to undo stack along with the index of the new constraint.
+  undoLog.emplace_back(UndoOp::DEALLOCATE, ~(con.size() - 1));
+  return con.size() - 1;
+}
+
+unsigned Simplex::getNumColumns() const { return nCol; }
+
+unsigned Simplex::getNumRows() const { return nRow; }
+
+// Normalize the row by removing common factors that are common between the
+// denominator and all the numerator coefficients.
+void Simplex::normalizeRow(unsigned row) {
+  int64_t gcd = 0;
+  for (unsigned col = 0; col < nCol; ++col) {
+    if (gcd == 1)
+      break;
+    gcd = llvm::greatestCommonDivisor(gcd, std::abs(tableau(row, col)));
+  }
+  for (unsigned col = 0; col < nCol; ++col)
+    tableau(row, col) /= gcd;
+}
+
+bool Simplex::signMatchesDirection(int64_t elem, Direction direction) const {
+  assert(elem != 0 && "elem should not be 0");
+  return direction == Direction::UP ? elem > 0 : elem < 0;
+}
+
+Simplex::Direction Simplex::flippedDirection(Direction direction) const {
+  return direction == Direction::UP ? Direction::DOWN : Direction::UP;
+}
+
+// Find a pivot to change the sample value of `row` in the specified direction.
+// The returned pivot row will involve `row` if and only if the unknown is
+// unbounded in the specified direction.
+//
+// To increase (resp. decrease) the value of a row, we need to find a live
+// column with a non-zero coefficient. If the coefficient is positive, we need
+// to increase (decrease) the value of the column, and if the coefficient is
+// negative, we need to decrease (increase) the value of the column. Also,
+// we cannot decrease the sample value of restricted columns.
+//
+// If multiple columns are valid, we break ties by considering a lexicographic
+// ordering where we prefer unknowns with lower index.
+llvm::Optional<std::pair<unsigned, unsigned>>
+Simplex::findPivot(int row, Direction direction) const {
+  llvm::Optional<unsigned> col;
+  for (unsigned j = 2; j < nCol; ++j) {
+    int64_t elem = tableau(row, j);
+    if (elem == 0)
+      continue;
+
+    if (unknownFromColumn(j).restricted &&
+        !signMatchesDirection(elem, direction))
+      continue;
+    if (!col || colVar[j] < colVar[*col])
+      col = j;
+  }
+
+  if (!col)
+    return {};
+
+  Direction newDirection =
+      tableau(row, *col) < 0 ? flippedDirection(direction) : direction;
+  auto opt = findPivotRow(row, newDirection, *col);
+  return std::pair<unsigned, unsigned>{opt.getValueOr(row), *col};
+}
+
+void Simplex::swapRowWithCol(unsigned row, unsigned col) {
+  std::swap(rowVar[row], colVar[col]);
+  Unknown *uCol = &unknownFromColumn(col);
+  Unknown *uRow = &unknownFromRow(row);
+  uCol->ownsRow = false;
+  uRow->ownsRow = true;
+  uCol->pos = col;
+  uRow->pos = row;
+}
+
+void Simplex::pivot(const std::pair<unsigned, unsigned> &p) {
+  pivot(p.first, p.second);
+}
+
+// Pivot pivotRow and pivotCol.
+//
+// Let R be the pivot row unknown and let C be the pivot col unknown.
+// Since initially R = a*C + sum b_i * X_i
+// (where the sum is over the other column's unknowns, x_i)
+// C = (R - (sum b_i * X_i))/a
+//
+// Let u be some other row unknown.
+// u = c*C + sum d_i * X_i
+// So u = c*(R - sum b_i * X_i)/a + sum d_i * X_i
+//
+// This results in the following transform:
+//            pivot col    other col                   pivot col    other col
+// pivot row     a             b       ->   pivot row     1/a         -b/a
+// other row     c             d            other row     c/a        d - bc/a
+//
+// Taking into account the common denominators p and q:
+//
+//            pivot col    other col                    pivot col   other col
+// pivot row     a/p          b/p     ->   pivot row      p/a         -b/a
+// other row     c/q          d/q          other row     cp/aq    (da - bc)/aq
+//
+// The pivot row transform is accomplished be swapping a with the pivot row's
+// common denominator and negating the pivot row except for the pivot column
+// element.
+void Simplex::pivot(unsigned pivotRow, unsigned pivotCol) {
+  assert(pivotCol >= 2 && "Refusing to pivot invalid column");
+
+  swapRowWithCol(pivotRow, pivotCol);
+  std::swap(tableau(pivotRow, 0), tableau(pivotRow, pivotCol));
+  // We need to negate the whole pivot row except for the pivot column.
+  if (tableau(pivotRow, 0) < 0) {
+    // If the denominator is negative, we negate the row by simply negating the
+    // denominator.
+    tableau(pivotRow, 0) = -tableau(pivotRow, 0);
+    tableau(pivotRow, pivotCol) = -tableau(pivotRow, pivotCol);
+  } else {
+    for (unsigned col = 1; col < nCol; ++col) {
+      if (col == pivotCol)
+        continue;
+      tableau(pivotRow, col) = -tableau(pivotRow, col);
+    }
+  }
+  normalizeRow(pivotRow);
+
+  for (unsigned row = 0; row < nRow; ++row) {
+    if (row == pivotRow)
+      continue;
+    if (tableau(row, pivotCol) == 0) // Nothing to do.
+      continue;
+    tableau(row, 0) *= tableau(pivotRow, 0);
+    for (unsigned j = 1; j < nCol; ++j) {
+      if (j == pivotCol)
+        continue;
+      // Add rather than subtract because the pivot row has been negated.
+      tableau(row, j) = tableau(row, j) * tableau(pivotRow, 0) +
+                        tableau(row, pivotCol) * tableau(pivotRow, j);
+    }
+    tableau(row, pivotCol) *= tableau(pivotRow, pivotCol);
+    normalizeRow(row);
+  }
+}
+
+// Perform pivots until the unknown has a non-negative sample value or until
+// no more upward pivots can be performed. Return the sign of the final sample
+// value.
+bool Simplex::restoreRow(Unknown &u) {
+  assert(u.ownsRow && "unknown should be in row position");
+
+  while (tableau(u.pos, 1) < 0) {
+    auto p = findPivot(u.pos, Direction::UP);
+    if (!p)
+      break;
+
+    pivot(*p);
+    if (!u.ownsRow)
+      return 1; // the unknown is unbounded above.
+  }
+  return tableau(u.pos, 1) >= 0;
+}
+
+// Find a row that can be used to pivot the column in the specified direction.
+// This returns an empty optional if and only if the column is unbounded in the
+// specified direction (ignoring skipRow, if skipRow is set).
+//
+// If skipRow is set, this row is not considered, and (if it is restricted) its
+// restriction may be violated by the returned pivot. Usually, skipRow is set
+// because we don't want to move it to column position unless it is unbounded,
+// and we are either trying to increase the value of skipRow or explicitly
+// trying to make skipRow negative, so we are not concerned about this.
+//
+// If the direction is up (resp. down) and a restricted row has a negative
+// (positive) coefficient for the column, then this row imposes a bound on how
+// much the sample value of the column can change. Such a row with constant term
+// c and coefficient f for the column imposes a bound of c/|f| on the change in
+// sample value (in the specified direction). (note that c is non-negative here
+// since the row is restricted and the tableau is consistent)
+//
+// We iterate through the rows and pick the row which imposes the most stringent
+// bound, since pivoting with a row changes the row's sample value to 0 and
+// hence saturates the bound it imposes. We break ties between rows that impose
+// the same bound by considering a lexicographic ordering where we prefer
+// unknowns with lower index value.
+llvm::Optional<unsigned> Simplex::findPivotRow(llvm::Optional<unsigned> skipRow,
+                                               Direction direction,
+                                               unsigned col) const {
+  llvm::Optional<unsigned> retRow;
+  int64_t retElem, retConst;
+  for (unsigned row = 0; row < nRow; ++row) {
+    if (skipRow && row == *skipRow)
+      continue;
+    auto elem = tableau(row, col);
+    if (elem == 0)
+      continue;
+    if (!unknownFromRow(row).restricted)
+      continue;
+    if (signMatchesDirection(elem, direction))
+      continue;
+    auto constTerm = tableau(row, 1);
+
+    if (!retRow) {
+      retRow = row;
+      retElem = elem;
+      retConst = constTerm;
+      continue;
+    }
+
+    auto diff = retConst * elem - constTerm * retElem;
+    if ((diff == 0 && rowVar[row] < rowVar[*retRow]) ||
+        (diff != 0 && !signMatchesDirection(diff, direction))) {
+      retRow = row;
+      retElem = elem;
+      retConst = constTerm;
+    }
+  }
+  return retRow;
+}
+
+bool Simplex::isEmpty() const { return empty; }
+
+void Simplex::swapRows(unsigned i, unsigned j) {
+  if (i == j)
+    return;
+  tableau.swapRows(i, j);
+  std::swap(rowVar[i], rowVar[j]);
+  unknownFromRow(i).pos = i;
+  unknownFromRow(j).pos = j;
+}
+
+// Mark this tableau empty and push an entry to the undo stack.
+void Simplex::markEmpty() {
+  undoLog.emplace_back(UndoOp::UNMARK_EMPTY, llvm::Optional<int>());
+  empty = true;
+}
+
+// Add an inequality to the tableau. If coeffs is c_0, c_1, ... c_n, where n
+// is the curent number of variables, then the corresponding inequality is
+// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} >= 0.
+//
+// We add the inequality and mark it as restricted. We then try to make its
+// sample value non-negative. If this is not possible, the tableau has become
+// empty and we mark it as such.
+void Simplex::addInequality(ArrayRef<int64_t> coeffs) {
+  unsigned conIndex = addRow(coeffs);
+  Unknown &u = con[conIndex];
+  u.restricted = true;
+  bool success = restoreRow(u);
+  if (!success)
+    markEmpty();
+}
+
+// Add an equality to the tableau. If coeffs is c_0, c_1, ... c_n, where n
+// is the curent number of variables, then the corresponding equality is
+// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} == 0.
+//
+// We simply add two opposing inequalities, which force the expression to
+// be zero.
+void Simplex::addEquality(ArrayRef<int64_t> coeffs) {
+  addInequality(coeffs);
+  SmallVector<int64_t, 64> negatedCoeffs;
+  for (auto coeff : coeffs)
+    negatedCoeffs.emplace_back(-coeff);
+  addInequality(negatedCoeffs);
+}
+
+unsigned Simplex::numberVariables() const { return var.size(); }
+unsigned Simplex::numberConstraints() const { return con.size(); }
+
+// Return a snapshot of the curent state. This is just the current size of the
+// undo log.
+unsigned Simplex::getSnapshot() const { return undoLog.size(); }
+
+void Simplex::undoOp(UndoOp op, llvm::Optional<int> index) {
+  if (op == UndoOp::DEALLOCATE) {
+    assert(index.hasValue() &&
+           "DEALLOCATE undo entry must be accompanied by an index");
+
+    assert(index < 0 && "Unknown to be deallocated must be a constraint");
+    Unknown &unknown = unknownFromIndex(*index);
+
+    if (!unknown.ownsRow) {
+      unsigned column = unknown.pos;
+      llvm::Optional<unsigned> row;
+
+      // Try to find any pivot row for this column that preserves tableau
+      // consistency (except possibly the column itself, which is going to be
+      // deallocated anyway).
+      //
+      // If no pivot row is found in either direction,
+      // then the column is unbounded in both directions and we are free to
+      // perform any pivot at all. To do this, we just need to find any row with
+      // a non-zero coefficient for the column.
+      if (auto maybeRow = findPivotRow({}, Direction::UP, column))
+        row = *maybeRow;
+      else if (auto maybeRow = findPivotRow({}, Direction::DOWN, column))
+        row = *maybeRow;
+      else {
+        // The loop doesn't find a pivot row only if the column has zero
+        // coefficients for every row. But the unknown is a constraint,
+        // so it was added initially as a row. Such a row could never have been
+        // pivoted to a column. So a pivot row will always be found.
+        for (unsigned i = 0; i < nRow; ++i) {
+          if (tableau(i, column) != 0) {
+            row = i;
+            break;
+          }
+        }
+      }
+      assert(row.hasValue() && "No pivot row found!");
+      pivot(*row, column);
+    }
+
+    // Move this unknown to the last row and remove the last row from the
+    // tableau.
+    swapRows(unknown.pos, nRow - 1);
+    // It is not strictly necessary to resize tableau, but for now we maintain
+    // the invariant that the actual size of the tableau is (nRow, nCol).
+    tableau.resize(nRow - 1, nCol);
+    nRow--;
+    rowVar.pop_back();
+    con.pop_back();
+  } else if (op == UndoOp::UNMARK_EMPTY) {
+    empty = false;
+  }
+}
+
+// Rollback to the specified snapshot.
+//
+// We undo all the log entries until the log size when the snapshot was taken
+// is reached.
+void Simplex::rollback(unsigned snapshot) {
+  while (undoLog.size() > snapshot) {
+    auto entry = undoLog.back();
+    undoLog.pop_back();
+    undoOp(entry.first, entry.second);
+  }
+}
+
+llvm::Optional<Fraction<int64_t>>
+Simplex::computeRowOptimum(Direction direction, unsigned row) {
+  while (auto maybePivot = findPivot(row, direction)) {
+    if (maybePivot->first == row)
+      return {};
+    pivot(*maybePivot);
+  }
+  return Fraction<int64_t>(tableau(row, 1), tableau(row, 0));
+}
+
+llvm::Optional<Fraction<int64_t>>
+Simplex::computeOptimum(Direction direction, ArrayRef<int64_t> coeffs) {
+  assert(!empty && "Tableau should not be empty");
+
+  auto snap = getSnapshot();
+  unsigned conIndex = addRow(coeffs);
+  unsigned row = con[conIndex].pos;
+  auto optimum = computeRowOptimum(direction, row);
+  rollback(snap);
+  return optimum;
+}
+
+template <typename T>
+std::vector<T> concat(const std::vector<T> &v, const std::vector<T> &w) {
+  auto result = v;
+  result.insert(result.end(), w.begin(), w.end());
+  return std::move(result);
+}
+
+// Make a tableau to represent a pair of points in the original tableau.
+//
+// The product constraints and variables are stored as: first A's, then B's.
+//
+// The product tableau has row layout:
+//   A's rows, B's rows.
+//
+// It has column layout:
+//   denominator, constant, A's columns, B's columns.
+Simplex Simplex::makeProduct(const Simplex &a, const Simplex &b) {
+  unsigned numVar = a.numberVariables() + b.numberVariables();
+  unsigned numCon = a.numberConstraints() + b.numberConstraints();
+  Simplex result(numVar);
+
+  result.tableau.resize(numCon, 2 + numVar);
+  result.empty = a.empty || b.empty;
+
+  auto indexFromBIndex = [&](int index) {
+    return index >= 0 ? a.numberVariables() + index
+                      : ~(a.numberConstraints() + ~index);
+  };
+
+  result.con = concat(a.con, b.con);
+  result.var = concat(a.var, b.var);
+  for (unsigned i = 2; i < a.nCol; i++) {
+    result.colVar.push_back(a.colVar[i]);
+    result.unknownFromIndex(result.colVar.back()).pos =
+        result.colVar.size() - 1;
+  }
+  for (unsigned i = 2; i < b.nCol; i++) {
+    result.colVar.push_back(indexFromBIndex(b.colVar[i]));
+    result.unknownFromIndex(result.colVar.back()).pos =
+        result.colVar.size() - 1;
+  }
+
+  auto appendRowFromA = [&](unsigned row) {
+    for (unsigned col = 0; col < a.nCol; col++)
+      result.tableau(result.nRow, col) = a.tableau(row, col);
+    result.rowVar.push_back(a.rowVar[row]);
+    result.unknownFromIndex(result.rowVar.back()).pos =
+        result.rowVar.size() - 1;
+    result.nRow++;
+  };
+
+  // Also fixes the corresponding entry in rowVar and var/con.
+  auto appendRowFromB = [&](unsigned row) {
+    result.tableau(result.nRow, 0) = a.tableau(row, 0);
+    result.tableau(result.nRow, 1) = a.tableau(row, 1);
+
+    unsigned offset = a.nCol - 2;
+    for (unsigned col = 2; col < b.nCol; col++)
+      result.tableau(result.nRow, offset + col) = b.tableau(row, col);
+    result.rowVar.push_back(indexFromBIndex(b.rowVar[row]));
+    result.unknownFromIndex(result.rowVar.back()).pos =
+        result.rowVar.size() - 1;
+    result.nRow++;
+  };
+
+  for (unsigned row = 0; row < a.nRow; row++)
+    appendRowFromA(row);
+  for (unsigned row = 0; row < b.nRow; row++)
+    appendRowFromB(row);
+
+  return result;
+}
+
+llvm::Optional<std::vector<int64_t>> Simplex::getSamplePointIfIntegral() const {
+  if (empty)
+    return {};
+
+  std::vector<int64_t> sample;
+  for (const Unknown &u : var) {
+    if (!u.ownsRow)
+      sample.push_back(0);
+    else {
+      if (tableau(u.pos, 1) % tableau(u.pos, 0) != 0)
+        return {};
+      sample.push_back(tableau(u.pos, 1) / tableau(u.pos, 0));
+    }
+  }
+  return sample;
+}
+
+// Given a simplex for a polytope, construct a new simplex whose variables are
+// identified with a pair of points (x, y) in the original polytope. Supports
+// some operations needed for general basis reduction. In what follows, <x, y>
+// denotes the dot product of the vectors x and y.
+class GBRSimplex {
+public:
+  GBRSimplex(const Simplex &originalSimplex)
+      : simplex(Simplex::makeProduct(originalSimplex, originalSimplex)),
+        simplexConstraintOffset(simplex.numberConstraints()) {}
+
+  // Add an equality <dir, x - y> = 0.
+  void addEqualityForDirection(const std::vector<int64_t> &dir) {
+    snapshotStack.push_back(simplex.getSnapshot());
+    simplex.addEquality(getCoeffsForDirection(dir));
+  }
+
+  // Compute max(<dir, x - y>) and save the dual variables for only the
+  // direction constraints to `dual`.
+  Fraction<int64_t> computeWidthAndDuals(const std::vector<int64_t> &dir,
+                                         std::vector<int64_t> &dual,
+                                         int64_t &dualDenom) {
+    unsigned snap = simplex.getSnapshot();
+    unsigned conIndex = simplex.addRow(getCoeffsForDirection(dir));
+    unsigned row = simplex.con[conIndex].pos;
+    auto width = simplex.computeRowOptimum(Simplex::Direction::UP, row);
+    assert(width && "Width should not be unbounded!");
+    dualDenom = simplex.tableau(row, 0);
+    dual.clear();
+    // The increment is i += 2 because equalities are added as two inequalities,
+    // one positive and one negative. We only want to process the positive ones.
+    for (unsigned i = simplexConstraintOffset; i < conIndex; i += 2) {
+      if (simplex.con[i].ownsRow)
+        dual.push_back(0);
+      else {
+        // The dual variable is the negative of the row coefficient.
+        dual.push_back(-simplex.tableau(row, simplex.con[i].pos));
+      }
+    }
+    simplex.rollback(snap);
+    return *width;
+  }
+
+  // Remove the last equality that was added through addEqualityForDirection.
+  void removeLastEquality() {
+    assert(!snapshotStack.empty() && "Snapshot stack is empty!");
+    simplex.rollback(snapshotStack.back());
+    snapshotStack.pop_back();
+  }
+
+private:
+  // Returns coefficients for the expression <dir, x - y>.
+  std::vector<int64_t> getCoeffsForDirection(ArrayRef<int64_t> dir) {
+    assert(2 * dir.size() == simplex.numberVariables() &&
+           "Direction vector has wrong dimensionality");
+    std::vector<int64_t> coeffs;
+    for (unsigned i = 0; i < dir.size(); i++)
+      coeffs.emplace_back(dir[i]);
+    for (unsigned i = 0; i < dir.size(); i++)
+      coeffs.emplace_back(-dir[i]);
+    coeffs.emplace_back(0); // constant term
+    return coeffs;
+  }
+
+  Simplex simplex;
+  // The first index of the equality constraints, the index immediately after
+  // the last constraint in the initial product simplex.
+  unsigned simplexConstraintOffset;
+  // A stack of snapshots, used for rolling back.
+  std::vector<unsigned> snapshotStack;
+};
+
+// Let b_{level}, b_{level + 1}, ... b_n be the current basis.
+// Let f_i(v) = max <v, x - y> where x and y are points in the original polytope
+// and <b_j, x - y> = 0 is satisfied for all j < i. (here <u, v> denotes the
+// inner product)
+//
+// In every iteration, we first replace b_{i+1} with b_{i+1} + u*b_i, where u is
+// the integer such that f_i(b_{i+1} + u*b_i) minimized. Let alpha be the dual
+// variable associated with the constraint <b_i, x - y> = 0 when computing
+// F_{i+1}(b_{i+1}). alpha must be the minimizing value of u, if it were allowed
+// to be rational. Due to convexity, the minimizing integer value is either
+// floor(alpha) or ceil(alpha), so we just need to check which of these gives a
+// lower F_{i+1} value. If alpha turned out to be an integer, then u = alpha.
+//
+// Now if f_i(b_{i+1}) < eps * f_i(b_i), we swap b_i and (the new) b_{i + 1}
+// and decrement i (unless i = level, in which case we stay at the same i).
+// Otherwise, we increment i. We use eps = 0.75.
+//
+// In an iteration we need to compute:
+//
+// Some of the required F and alpha values may already be known. We cache the
+// known values and reuse them if possible. In particular:
+//
+// When we set b_{i+1} to b_{i+1} + u*b_i, no F values are changed since we only
+// added a multiple of b_i to b_{i+1}. In particular F_{i+1}(b_{i+1})
+// = min f_i(b_{i+1} + alpha * b_i) so adding u*b_i to b_{i+1} does not
+// change this. Also <b_i, x - y> = 0 and <b_{i+1}, x - y> = 0 already
+// imply <b_{i+1} + u*b_i, x - y> = 0, so the constraints are unchanged.
+//
+// When we decrement i, we swap b_i and b_{i+1}. In the following paragraphs we
+// always refer to the final vector b_{i+1} after updating). But note that when
+// we come to the end of an iteration we always know f_i(b_{i+1}), so we just
+// need to update the cached value to reflect this. However the value of
+// F_{i+1}(b_i) is not already known, so if there was a stale value of F[i+1] in
+// the cache we remove this. Moreover, the iteration after decrementing will
+// want the dual variables from this computation so we cache this when we
+// compute the minimizing value of u.
+//
+// If alpha turns out to be an integer, then we never compute f_i(b_{i+1}) in
+// this iteration. But in this case, F_{i+1}(b_{i+1}) = F_{i+1}(b'_{i+1}) where
+// b'_{i+1} is the vector before updating. Therefore, we can update the cache
+// with this value. Furthermore, we can just inherit the dual variables from
+// this computation.
+//
+// When incrementing i we do not make any changes to the basis so no
+// invalidation occurs.
+void Simplex::reduceBasis(Matrix<int64_t> &basis, unsigned level) {
+  const Fraction<int64_t> epsilon(3, 4);
+
+  if (level == basis.getNumRows() - 1)
+    return;
+
+  GBRSimplex gbrSimplex(*this);
+  std::vector<Fraction<int64_t>> f;
+  std::vector<int64_t> alpha;
+  int64_t alphaDenom;
+  auto findUAndGetFCandidate = [&](unsigned i) -> Fraction<int64_t> {
+    assert(i < level + alpha.size() && "alpha_i is not known!");
+
+    int64_t u = floorDiv(alpha[i - level], alphaDenom);
+    basis.addToRow(i, i + 1, u);
+    if (alpha[i - level] % alphaDenom != 0) {
+      std::vector<int64_t> uAlpha[2];
+      int64_t uAlphaDenom[2];
+      Fraction<int64_t> fAtI[2];
+
+      // Initially u is floor(alpha) and basis reflects this.
+      fAtI[0] = gbrSimplex.computeWidthAndDuals(basis.getRow(i + 1), uAlpha[0],
+                                                uAlphaDenom[0]);
+
+      // Now try ceil(alpha), i.e. floor(alpha) + 1.
+      ++u;
+      basis.addToRow(i, i + 1, 1);
+      fAtI[1] = gbrSimplex.computeWidthAndDuals(basis.getRow(i + 1), uAlpha[1],
+                                                uAlphaDenom[1]);
+
+      int j = fAtI[0] < fAtI[1] ? 0 : 1;
+      if (j == 0)
+        // Subtract 1 to go from u = ceil(alpha) back to floor(alpha).
+        basis.addToRow(i, i + 1, -1);
+
+      alpha = std::move(uAlpha[j]);
+      alphaDenom = uAlphaDenom[j];
+      return fAtI[j];
+    }
+    assert(i + 1 - level < f.size() && "F_{i+1} wasn't saved");
+    // When alpha minimizes f_i(b_{i+1} + alpha*b_i), this is equal to
+    // F_{i+1}(b_{i+1}).
+    return f[i + 1 - level];
+  };
+
+  // In the ith iteration of the loop, gbrSimplex has constraints for directions
+  // from `level` to i - 1.
+  unsigned i = level;
+  while (i < basis.getNumRows() - 1) {
+    Fraction<int64_t> fAtICandidate; // f_i(b_{i+1} + u*b_i)
+    if (i >= level + f.size()) {
+      // We don't even know the value of f_i(b_i), so let's find that first.
+      // We have to do this first since later we assume that f already contains
+      // values up to and including i.
+
+      assert((i == 0 || i - 1 < level + f.size()) &&
+             "We are at level i but we don't know the value of F_{i-1}");
+
+      // We don't actually use these duals at all, but it doesn't matter
+      // because this case should only occur when i is level, and there are no
+      // duals in that case anyway.
+      assert(i == level && "This case should only occur when i == level");
+      f.push_back(
+          gbrSimplex.computeWidthAndDuals(basis.getRow(i), alpha, alphaDenom));
+    }
+
+    if (i >= level + alpha.size()) {
+      assert(i + 1 >= level + f.size() && "We don't know alpha_i but we know "
+                                          "F_{i+1}, this should never happen");
+      // We don't know alpha for our level, so let's find it.
+      gbrSimplex.addEqualityForDirection(basis.getRow(i));
+      f.push_back(gbrSimplex.computeWidthAndDuals(basis.getRow(i + 1), alpha,
+                                                  alphaDenom));
+      gbrSimplex.removeLastEquality();
+    }
+
+    fAtICandidate = findUAndGetFCandidate(i);
+
+    if (fAtICandidate < epsilon * f[i - level]) {
+      basis.swapRows(i, i + 1);
+      f[i - level] = fAtICandidate;
+      // The values of f_{i+1}(b_{i+1}) and higher may change after the swap,
+      // so we remove the cached values here.
+      f.resize(i - level + 1);
+      if (i == level) {
+        alpha.clear();
+        continue;
+      }
+
+      gbrSimplex.removeLastEquality();
+      i--;
+      continue;
+    }
+
+    alpha.clear();
+    gbrSimplex.addEqualityForDirection(basis.getRow(i));
+    i++;
+  }
+}
+
+// Try to find an integer sample point in the polytope.
+//
+// If such a point exists, this function returns it. Otherwise, it returns and
+// empty llvm::Optional.
+llvm::Optional<std::vector<int64_t>> Simplex::findIntegerSample() {
+  if (empty)
+    return {};
+
+  unsigned nDims = var.size();
+  Matrix<int64_t> basis = Matrix<int64_t>::getIdentityMatrix(nDims);
+  return findIntegerSampleRecursively(basis, 0);
+}
+
+// Search for an integer sample point using a branch and bound algorithm.
+//
+// Each row in the basis matrix is a vector, and the set of basis vectors should
+// span the space. Initially this is called with the identity matrix, i.e., the
+// basis vectors are just the variables.
+//
+// In every level, a value is assigned to the level-th basis vector, as follows.
+// Compute the minimum and maximum rational values of this direction.
+// If only one integer point lies in this range, constrain the variable to
+// have this value and recurse to the next variable.
+//
+// If the range has multiple values, perform general basis reduction via
+// reduceBasis and then compute the bounds again. Now we can't do any better
+// than this, so we just recurse on every integer value in this range.
+//
+// If the range contains no integer value, then of course the polytope is empty
+// for the current assignment of the values in previous levels, so return to
+// the previous level.
+//
+// If we reach the last level where all the variables have been assigned values
+// already, then we simply return the current sample point if it is integral, or
+// an empty llvm::Optional otherwise.
+llvm::Optional<std::vector<int64_t>>
+Simplex::findIntegerSampleRecursively(Matrix<int64_t> &basis, unsigned level) {
+  if (level == basis.getNumRows())
+    return getSamplePointIfIntegral();
+
+  std::vector<int64_t> basisCoeffVector = basis.getRow(level);
+  basisCoeffVector.emplace_back(0); // constant term
+
+  auto getBounds = [&]() -> std::pair<int64_t, int64_t> {
+    int64_t minRoundedUp;
+    if (auto opt = computeOptimum(Direction::DOWN, basisCoeffVector))
+      minRoundedUp = ceil(*opt);
+    else
+      llvm_unreachable("Tableau should not be unbounded");
+
+    int64_t maxRoundedDown;
+    if (auto opt = computeOptimum(Direction::UP, basisCoeffVector))
+      maxRoundedDown = floor(*opt);
+    else
+      llvm_unreachable("Tableau should not be unbounded");
+
+    return {minRoundedUp, maxRoundedDown};
+  };
+
+  int64_t minRoundedUp, maxRoundedDown;
+  std::tie(minRoundedUp, maxRoundedDown) = getBounds();
+
+  // Heuristic: if the sample point is integral at this point, just return it.
+  if (auto opt = getSamplePointIfIntegral())
+    return *opt;
+
+  if (minRoundedUp < maxRoundedDown) {
+    reduceBasis(basis, level);
+    std::tie(minRoundedUp, maxRoundedDown) = getBounds();
+  }
+
+  for (int64_t i = minRoundedUp; i <= maxRoundedDown; ++i) {
+    auto snapshot = getSnapshot();
+    basisCoeffVector.back() = -i;
+    // Add the constraint `basisCoeffVector = i`.
+    addEquality(basisCoeffVector);
+    if (auto opt = findIntegerSampleRecursively(basis, level + 1))
+      return *opt;
+    rollback(snapshot);
+  }
+
+  return {};
+}
+
+void Simplex::printUnknown(llvm::raw_ostream &os, const Unknown &u) const {
+  os << (u.ownsRow ? "r" : "c");
+  os << u.pos;
+  if (u.restricted)
+    os << " [>=0]";
+}
+
+void Simplex::print(llvm::raw_ostream &os) const {
+  os << "Dumping Simplex, rows = " << nRow << ", columns = " << nCol << "\n";
+  if (empty)
+    os << "Simplex marked empty!\n";
+  os << "var: ";
+  for (unsigned i = 0; i < var.size(); ++i) {
+    if (i > 0)
+      os << ", ";
+    printUnknown(os, var[i]);
+  }
+  os << "\ncon: ";
+  for (unsigned i = 0; i < con.size(); ++i) {
+    if (i > 0)
+      os << ", ";
+    printUnknown(os, con[i]);
+  }
+  os << '\n';
+  for (unsigned row = 0; row < nRow; ++row) {
+    if (row > 0)
+      os << ", ";
+    os << "r" << row << ": " << rowVar[row];
+  }
+  os << '\n';
+  os << "c0: denom, c1: const";
+  for (unsigned col = 2; col < nCol; ++col)
+    os << ", c" << col << ": " << colVar[col];
+  os << '\n';
+  for (unsigned row = 0; row < nRow; ++row) {
+    for (unsigned col = 0; col < nCol; ++col)
+      os << tableau(row, col) << '\t';
+    os << '\n';
+  }
+  os << '\n';
+}
+
+void Simplex::dump() const { print(llvm::errs()); }
+
+} // namespace mlir

mlir/unittests/Analysis/AffineStructuresTest.cpp

+//===- 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 <gmock/gmock.h>
+#include <gtest/gtest.h>
+
+namespace mlir {
+
+int64_t valueAt(ArrayRef<int64_t> expr, const std::vector<int64_t> &point) {
+  int64_t value = expr.back();
+  assert(expr.size() == 1 + point.size());
+  for (unsigned i = 0; i < point.size(); ++i)
+    value += expr[i] * point[i];
+  return value;
+}
+
+void checkSample(bool hasValue, const FlatAffineConstraints &fac) {
+  auto maybeSample = fac.findIntegerSample();
+  if (!hasValue)
+    EXPECT_FALSE(maybeSample.hasValue());
+  else {
+    ASSERT_TRUE(maybeSample.hasValue());
+    for (unsigned i = 0; i < fac.getNumEqualities(); ++i)
+      EXPECT_EQ(valueAt(fac.getEquality(i), *maybeSample), 0);
+    for (unsigned i = 0; i < fac.getNumInequalities(); ++i)
+      EXPECT_GE(valueAt(fac.getInequality(i), *maybeSample), 0);
+  }
+}
+
+FlatAffineConstraints
+makeFACFromConstraints(unsigned dims, std::vector<std::vector<int64_t>> ineqs,
+                       std::vector<std::vector<int64_t>> eqs) {
+  FlatAffineConstraints fac(ineqs.size(), eqs.size(), dims + 1, dims);
+  for (const auto &eq : eqs)
+    fac.addEquality(eq);
+  for (const auto &ineq : ineqs)
+    fac.addInequality(ineq);
+  return fac;
+}
+
+TEST(FlatAffineConstraintsTest, FindSampleTest) {
+  // bounded sets, only inequalities
+
+  // 1 <= 5x and 5x <= 4 (no solution)
+  checkSample(false, makeFACFromConstraints(1, {{5, -1}, {-5, 4}}, {}));
+
+  // 1 <= 5x and 5x <= 9 (solution: x = 1)
+  checkSample(true, makeFACFromConstraints(1, {{5, -1}, {-5, 9}}, {}));
+
+  // bounded sets with equalities
+
+  // x >= 8 and 40 >= y and x = y
+  checkSample(
+      true, makeFACFromConstraints(2, {{1, 0, -8}, {0, -1, 40}}, {{1, -1, 0}}));
+
+  // x <= 10 and y <= 10 and 10 <= z and x + 2y = 3z
+  // solution: x = y = z = 10.
+  checkSample(true, makeFACFromConstraints(
+                        3, {{-1, 0, 0, 10}, {0, -1, 0, 10}, {0, 0, 1, -10}},
+                        {{1, 2, -3, 0}}));
+
+  // 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, makeFACFromConstraints(
+                         3, {{-1, 0, 0, 10}, {0, -1, 0, 10}, {0, 0, 1, -11}},
+                         {{1, 2, -3, 0}}));
+
+  // 0 <= r and r <= 3 and 4q + r = 7
+  // Solution: q = 1, r = 3
+  checkSample(true,
+              makeFACFromConstraints(2, {{0, 1, 0}, {0, -1, 3}}, {{4, 1, -7}}));
+
+  // 4q + r = 7 and r = 0
+  // Solution: q = 1, r = 3
+  checkSample(false, makeFACFromConstraints(2, {}, {{4, 1, -7}, {0, 1, 0}}));
+}
+
+TEST(FlatAffineConstraintsTest, IsIntegerEmptyTest) {
+  // 1 <= 5x and 5x <= 4 (no solution)
+  EXPECT_TRUE(
+      makeFACFromConstraints(1, {{5, -1}, {-5, 4}}, {}).isIntegerEmpty());
+  // 1 <= 5x and 5x <= 9 (solution: x = 1)
+  EXPECT_FALSE(
+      makeFACFromConstraints(1, {{5, -1}, {-5, 9}}, {}).isIntegerEmpty());
+}
+
+} // namespace mlir

mlir/unittests/Analysis/CMakeLists.txt

+add_mlir_unittest(MLIRAnalysisTests
+  AffineStructuresTest.cpp
+)
+
+target_link_libraries(MLIRAnalysisTests
+  PRIVATE MLIRLoopAnalysis)
+
+add_subdirectory(Presburger)

mlir/unittests/Analysis/Presburger/CMakeLists.txt

+add_mlir_unittest(MLIRPresburgerTests
+  MatrixTest.cpp
+  SimplexTest.cpp
+)
+
+target_link_libraries(MLIRPresburgerTests
+  PRIVATE MLIRPresburger)

mlir/unittests/Analysis/Presburger/MatrixTest.cpp

+//===- MatrixTest.cpp - Tests for Matrix ------------------------*- 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/Presburger/Matrix.h"
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+
+namespace mlir {
+
+TEST(MatrixTest, ReadWrite) {
+  Matrix<int> mat(5, 5);
+  for (unsigned row = 0; row < 5; ++row)
+    for (unsigned col = 0; col < 5; ++col)
+      mat(row, col) = 10 * row + col;
+  for (unsigned row = 0; row < 5; ++row)
+    for (unsigned col = 0; col < 5; ++col)
+      EXPECT_EQ(mat(row, col), int(10 * row + col));
+}
+
+TEST(MatrixTest, SwapColumns) {
+  Matrix<int> mat(5, 5);
+  for (unsigned row = 0; row < 5; ++row)
+    for (unsigned col = 0; col < 5; ++col)
+      mat(row, col) = col == 3 ? 1 : 0;
+  mat.swapColumns(3, 1);
+  for (unsigned row = 0; row < 5; ++row)
+    for (unsigned col = 0; col < 5; ++col)
+      EXPECT_EQ(mat(row, col), col == 1 ? 1 : 0);
+
+  // swap around all the other columns, swap (1, 3) twice for no effect.
+  mat.swapColumns(3, 1);
+  mat.swapColumns(2, 4);
+  mat.swapColumns(1, 3);
+  mat.swapColumns(0, 4);
+
+  for (unsigned row = 0; row < 5; ++row)
+    for (unsigned col = 0; col < 5; ++col)
+      EXPECT_EQ(mat(row, col), col == 1 ? 1 : 0);
+}
+
+TEST(MatrixTest, SwapRows) {
+  Matrix<int> mat(5, 5);
+  for (unsigned row = 0; row < 5; ++row)
+    for (unsigned col = 0; col < 5; ++col)
+      mat(row, col) = row == 2 ? 1 : 0;
+  mat.swapRows(2, 0);
+  for (unsigned row = 0; row < 5; ++row)
+    for (unsigned col = 0; col < 5; ++col)
+      EXPECT_EQ(mat(row, col), row == 0 ? 1 : 0);
+
+  // swap around all the other rows, swap (2, 0) twice for no effect.
+  mat.swapRows(3, 4);
+  mat.swapRows(1, 4);
+  mat.swapRows(2, 0);
+  mat.swapRows(0, 2);
+
+  for (unsigned row = 0; row < 5; ++row)
+    for (unsigned col = 0; col < 5; ++col)
+      EXPECT_EQ(mat(row, col), row == 0 ? 1 : 0);
+}
+
+TEST(MatrixTest, Resize) {
+  Matrix<int> mat(5, 5);
+  EXPECT_EQ(mat.getNumRows(), 5u);
+  EXPECT_EQ(mat.getNumColumns(), 5u);
+  for (unsigned row = 0; row < 5; ++row)
+    for (unsigned col = 0; col < 5; ++col)
+      mat(row, col) = 10 * row + col;
+
+  mat.resize(3, 3);
+  EXPECT_EQ(mat.getNumRows(), 3u);
+  EXPECT_EQ(mat.getNumColumns(), 3u);
+  for (unsigned row = 0; row < 3; ++row)
+    for (unsigned col = 0; col < 3; ++col)
+      EXPECT_EQ(mat(row, col), int(10 * row + col));
+
+  mat.resize(5, 5);
+  EXPECT_EQ(mat.getNumRows(), 5u);
+  EXPECT_EQ(mat.getNumColumns(), 5u);
+  for (unsigned row = 0; row < 5; ++row)
+    for (unsigned col = 0; col < 5; ++col)
+      EXPECT_EQ(mat(row, col), row >= 3 || col >= 3 ? 0 : int(10 * row + col));
+}
+} // namespace mlir

mlir/unittests/Analysis/Presburger/SimplexTest.cpp

+//===- SimplexTest.cpp - Tests for Simplex ----------------------*- 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/Presburger/Simplex.h"
+
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+
+namespace mlir {
+
+TEST(SimplexTest, emptyRollback) {
+  Simplex tab(2);
+  // (u - v) >= 0
+  tab.addInequality({1, -1, 0});
+  EXPECT_FALSE(tab.isEmpty());
+
+  auto snap = tab.getSnapshot();
+  // (u - v) <= -1
+  tab.addInequality({-1, 1, -1});
+  EXPECT_TRUE(tab.isEmpty());
+  tab.rollback(snap);
+  EXPECT_FALSE(tab.isEmpty());
+}
+
+TEST(SimplexTest, addEquality_separate) {
+  Simplex tab(1);
+  tab.addInequality({1, -1}); // x >= 1
+  ASSERT_FALSE(tab.isEmpty());
+  tab.addEquality({1, 0}); // x == 0
+  EXPECT_TRUE(tab.isEmpty());
+}
+
+void expectInequalityMakesTabEmpty(Simplex &tab, ArrayRef<int64_t> coeffs,
+                                   bool expect) {
+  ASSERT_FALSE(tab.isEmpty());
+  auto snapshot = tab.getSnapshot();
+  tab.addInequality(coeffs);
+  EXPECT_EQ(tab.isEmpty(), expect);
+  tab.rollback(snapshot);
+}
+
+TEST(SimplexTest, addInequality_rollback) {
+  Simplex tab(3);
+
+  std::vector<int64_t> coeffs[]{{1, 0, 0, 0},   // u >= 0
+                                {-1, 0, 0, 0},  // u <= 0
+                                {1, -1, 1, 0},  // u - v + w >= 0
+                                {1, 1, -1, 0}}; // u + v - w >= 0
+  // The above constraints force u = 0 and v = w.
+  // The constraints below violate v = w.
+  std::vector<int64_t> checkCoeffs[]{{0, 1, -1, -1},  // v - w >= 1
+                                     {0, -1, 1, -1}}; // v - w <= -1
+  for (int run = 0; run < 4; run++) {
+    auto snapshot = tab.getSnapshot();
+
+    expectInequalityMakesTabEmpty(tab, checkCoeffs[0], false);
+    expectInequalityMakesTabEmpty(tab, checkCoeffs[1], false);
+
+    for (int i = 0; i < 4; i++)
+      tab.addInequality(coeffs[(run + i) % 4]);
+
+    expectInequalityMakesTabEmpty(tab, checkCoeffs[0], true);
+    expectInequalityMakesTabEmpty(tab, checkCoeffs[1], true);
+
+    tab.rollback(snapshot);
+    EXPECT_EQ(tab.numberConstraints(), 0u);
+
+    expectInequalityMakesTabEmpty(tab, checkCoeffs[0], false);
+    expectInequalityMakesTabEmpty(tab, checkCoeffs[1], false);
+  }
+}
+} // namespace mlir

mlir/unittests/CMakeLists.txt

 add_custom_target(MLIRUnitTests)
 set_target_properties(MLIRUnitTests PROPERTIES FOLDER "MLIR Tests")
 
 function(add_mlir_unittest test_dirname)
   add_unittest(MLIRUnitTests ${test_dirname} ${ARGN})
 endfunction()
 
+add_subdirectory(Analysis)
 add_subdirectory(Dialect)
 add_subdirectory(IR)
 add_subdirectory(Pass)
 add_subdirectory(SDBM)
 add_subdirectory(TableGen)