1	//===- Simplex.cpp - MLIR Simplex Class -----------------------------------===//
2	//
3	// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4	// See https://llvm.org/LICENSE.txt for license information.
5	// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6	//
7	//===----------------------------------------------------------------------===//
8
9	#include "mlir/Analysis/Presburger/Simplex.h"
10	#include "mlir/Analysis/Presburger/Fraction.h"
11	#include "mlir/Analysis/Presburger/IntegerRelation.h"
12	#include "mlir/Analysis/Presburger/Matrix.h"
13	#include "mlir/Analysis/Presburger/PresburgerSpace.h"
14	#include "mlir/Analysis/Presburger/Utils.h"
15	#include "llvm/ADT/DynamicAPInt.h"
16	#include "llvm/ADT/STLExtras.h"
17	#include "llvm/ADT/SmallBitVector.h"
18	#include "llvm/ADT/SmallVector.h"
19	#include "llvm/Support/Compiler.h"
20	#include "llvm/Support/ErrorHandling.h"
21	#include "llvm/Support/LogicalResult.h"
22	#include "llvm/Support/raw_ostream.h"
23	#include <cassert>
24	#include <functional>
25	#include <limits>
26	#include <optional>
27	#include <tuple>
28	#include <utility>
29
30	using namespace mlir;
31	using namespace presburger;
32
33	using Direction = Simplex::Direction;
34
35	const int nullIndex = std::numeric_limits<int>::max();
36
37	// Return a + scale*b;
38	LLVM_ATTRIBUTE_UNUSED
39	static SmallVector<DynamicAPInt, 8>
40	scaleAndAddForAssert(ArrayRef<DynamicAPInt> a, const DynamicAPInt &scale,
41	ArrayRef<DynamicAPInt> b) {
42	assert(a.size() == b.size());
43	SmallVector<DynamicAPInt, 8> res;
44	res.reserve(N: a.size());
45	for (unsigned i = 0, e = a.size(); i < e; ++i)
46	res.emplace_back(Args: a[i] + scale * b[i]);
47	return res;
48	}
49
50	SimplexBase::SimplexBase(unsigned nVar, bool mustUseBigM)
51	: usingBigM(mustUseBigM), nRedundant(0), nSymbol(0),
52	tableau(0, getNumFixedCols() + nVar), empty(false) {
53	var.reserve(N: nVar);
54	colUnknown.reserve(N: nVar + 1);
55	colUnknown.insert(I: colUnknown.begin(), NumToInsert: getNumFixedCols(), Elt: nullIndex);
56	for (unsigned i = 0; i < nVar; ++i) {
57	var.emplace_back(Args: Orientation::Column, /*restricted=*/Args: false,
58	/*pos=*/Args: getNumFixedCols() + i);
59	colUnknown.emplace_back(Args&: i);
60	}
61	}
62
63	SimplexBase::SimplexBase(unsigned nVar, bool mustUseBigM,
64	const llvm::SmallBitVector &isSymbol)
65	: SimplexBase(nVar, mustUseBigM) {
66	assert(isSymbol.size() == nVar && "invalid bitmask!");
67	// Invariant: nSymbol is the number of symbols that have been marked
68	// already and these occupy the columns
69	// [getNumFixedCols(), getNumFixedCols() + nSymbol).
70	for (unsigned symbolIdx : isSymbol.set_bits()) {
71	var[symbolIdx].isSymbol = true;
72	swapColumns(i: var[symbolIdx].pos, j: getNumFixedCols() + nSymbol);
73	++nSymbol;
74	}
75	}
76
77	const Simplex::Unknown &SimplexBase::unknownFromIndex(int index) const {
78	assert(index != nullIndex && "nullIndex passed to unknownFromIndex");
79	return index >= 0 ? var[index] : con[~index];
80	}
81
82	const Simplex::Unknown &SimplexBase::unknownFromColumn(unsigned col) const {
83	assert(col < getNumColumns() && "Invalid column");
84	return unknownFromIndex(index: colUnknown[col]);
85	}
86
87	const Simplex::Unknown &SimplexBase::unknownFromRow(unsigned row) const {
88	assert(row < getNumRows() && "Invalid row");
89	return unknownFromIndex(index: rowUnknown[row]);
90	}
91
92	Simplex::Unknown &SimplexBase::unknownFromIndex(int index) {
93	assert(index != nullIndex && "nullIndex passed to unknownFromIndex");
94	return index >= 0 ? var[index] : con[~index];
95	}
96
97	Simplex::Unknown &SimplexBase::unknownFromColumn(unsigned col) {
98	assert(col < getNumColumns() && "Invalid column");
99	return unknownFromIndex(index: colUnknown[col]);
100	}
101
102	Simplex::Unknown &SimplexBase::unknownFromRow(unsigned row) {
103	assert(row < getNumRows() && "Invalid row");
104	return unknownFromIndex(index: rowUnknown[row]);
105	}
106
107	unsigned SimplexBase::addZeroRow(bool makeRestricted) {
108	// Resize the tableau to accommodate the extra row.
109	unsigned newRow = tableau.appendExtraRow();
110	assert(getNumRows() == getNumRows() && "Inconsistent tableau size");
111	rowUnknown.emplace_back(Args: ~con.size());
112	con.emplace_back(Args: Orientation::Row, Args&: makeRestricted, Args&: newRow);
113	undoLog.emplace_back(Args: UndoLogEntry::RemoveLastConstraint);
114	tableau(newRow, 0) = 1;
115	return newRow;
116	}
117
118	/// Add a new row to the tableau corresponding to the given constant term and
119	/// list of coefficients. The coefficients are specified as a vector of
120	/// (variable index, coefficient) pairs.
121	unsigned SimplexBase::addRow(ArrayRef<DynamicAPInt> coeffs,
122	bool makeRestricted) {
123	assert(coeffs.size() == var.size() + 1 &&
124	"Incorrect number of coefficients!");
125	assert(var.size() + getNumFixedCols() == getNumColumns() &&
126	"inconsistent column count!");
127
128	unsigned newRow = addZeroRow(makeRestricted);
129	tableau(newRow, 1) = coeffs.back();
130	if (usingBigM) {
131	// When the lexicographic pivot rule is used, instead of the variables
132	//
133	// x, y, z ...
134	//
135	// we internally use the variables
136	//
137	// M, M + x, M + y, M + z, ...
138	//
139	// where M is the big M parameter. As such, when the user tries to add
140	// a row ax + by + cz + d, we express it in terms of our internal variables
141	// as -(a + b + c)M + a(M + x) + b(M + y) + c(M + z) + d.
142	//
143	// Symbols don't use the big M parameter since they do not get lex
144	// optimized.
145	DynamicAPInt bigMCoeff(0);
146	for (unsigned i = 0; i < coeffs.size() - 1; ++i)
147	if (!var[i].isSymbol)
148	bigMCoeff -= coeffs[i];
149	// The coefficient to the big M parameter is stored in column 2.
150	tableau(newRow, 2) = bigMCoeff;
151	}
152
153	// Process each given variable coefficient.
154	for (unsigned i = 0; i < var.size(); ++i) {
155	unsigned pos = var[i].pos;
156	if (coeffs[i] == 0)
157	continue;
158
159	if (var[i].orientation == Orientation::Column) {
160	// If a variable is in column position at column col, then we just add the
161	// coefficient for that variable (scaled by the common row denominator) to
162	// the corresponding entry in the new row.
163	tableau(newRow, pos) += coeffs[i] * tableau(newRow, 0);
164	continue;
165	}
166
167	// If the variable is in row position, we need to add that row to the new
168	// row, scaled by the coefficient for the variable, accounting for the two
169	// rows potentially having different denominators. The new denominator is
170	// the lcm of the two.
171	DynamicAPInt lcm = llvm::lcm(A: tableau(newRow, 0), B: tableau(pos, 0));
172	DynamicAPInt nRowCoeff = lcm / tableau(newRow, 0);
173	DynamicAPInt idxRowCoeff = coeffs[i] * (lcm / tableau(pos, 0));
174	tableau(newRow, 0) = lcm;
175	for (unsigned col = 1, e = getNumColumns(); col < e; ++col)
176	tableau(newRow, col) =
177	nRowCoeff * tableau(newRow, col) + idxRowCoeff * tableau(pos, col);
178	}
179
180	tableau.normalizeRow(row: newRow);
181	// Push to undo log along with the index of the new constraint.
182	return con.size() - 1;
183	}
184
185	namespace {
186	bool signMatchesDirection(const DynamicAPInt &elem, Direction direction) {
187	assert(elem != 0 && "elem should not be 0");
188	return direction == Direction::Up ? elem > 0 : elem < 0;
189	}
190
191	Direction flippedDirection(Direction direction) {
192	return direction == Direction::Up ? Direction::Down : Simplex::Direction::Up;
193	}
194	} // namespace
195
196	/// We simply make the tableau consistent while maintaining a lexicopositive
197	/// basis transform, and then return the sample value. If the tableau becomes
198	/// empty, we return empty.
199	///
200	/// Let the variables be x = (x_1, ... x_n).
201	/// Let the basis unknowns be y = (y_1, ... y_n).
202	/// We have that x = A*y + b for some n x n matrix A and n x 1 column vector b.
203	///
204	/// As we will show below, A*y is either zero or lexicopositive.
205	/// Adding a lexicopositive vector to b will make it lexicographically
206	/// greater, so A*y + b is always equal to or lexicographically greater than b.
207	/// Thus, since we can attain x = b, that is the lexicographic minimum.
208	///
209	/// We have that every column in A is lexicopositive, i.e., has at least
210	/// one non-zero element, with the first such element being positive. Since for
211	/// the tableau to be consistent we must have non-negative sample values not
212	/// only for the constraints but also for the variables, we also have x >= 0 and
213	/// y >= 0, by which we mean every element in these vectors is non-negative.
214	///
215	/// Proof that if every column in A is lexicopositive, and y >= 0, then
216	/// A*y is zero or lexicopositive. Begin by considering A_1, the first row of A.
217	/// If this row is all zeros, then (A*y)_1 = (A_1)*y = 0; proceed to the next
218	/// row. If we run out of rows, A*y is zero and we are done; otherwise, we
219	/// encounter some row A_i that has a non-zero element. Every column is
220	/// lexicopositive and so has some positive element before any negative elements
221	/// occur, so the element in this row for any column, if non-zero, must be
222	/// positive. Consider (A*y)_i = (A_i)*y. All the elements in both vectors are
223	/// non-negative, so if this is non-zero then it must be positive. Then the
224	/// first non-zero element of A*y is positive so A*y is lexicopositive.
225	///
226	/// Otherwise, if (A_i)*y is zero, then for every column j that had a non-zero
227	/// element in A_i, y_j is zero. Thus these columns have no contribution to A*y
228	/// and we can completely ignore these columns of A. We now continue downwards,
229	/// looking for rows of A that have a non-zero element other than in the ignored
230	/// columns. If we find one, say A_k, once again these elements must be positive
231	/// since they are the first non-zero element in each of these columns, so if
232	/// (A_k)*y is not zero then we have that A*y is lexicopositive and if not we
233	/// add these to the set of ignored columns and continue to the next row. If we
234	/// run out of rows, then A*y is zero and we are done.
235	MaybeOptimum<SmallVector<Fraction, 8>> LexSimplex::findRationalLexMin() {
236	if (restoreRationalConsistency().failed()) {
237	markEmpty();
238	return OptimumKind::Empty;
239	}
240	return getRationalSample();
241	}
242
243	/// Given a row that has a non-integer sample value, add an inequality such
244	/// that this fractional sample value is cut away from the polytope. The added
245	/// inequality will be such that no integer points are removed. i.e., the
246	/// integer lexmin, if it exists, is the same with and without this constraint.
247	///
248	/// Let the row be
249	/// (c + coeffM*M + a_1*s_1 + ... + a_m*s_m + b_1*y_1 + ... + b_n*y_n)/d,
250	/// where s_1, ... s_m are the symbols and
251	/// y_1, ... y_n are the other basis unknowns.
252	///
253	/// For this to be an integer, we want
254	/// coeffM*M + a_1*s_1 + ... + a_m*s_m + b_1*y_1 + ... + b_n*y_n = -c (mod d)
255	/// Note that this constraint must always hold, independent of the basis,
256	/// becuse the row unknown's value always equals this expression, even if *we*
257	/// later compute the sample value from a different expression based on a
258	/// different basis.
259	///
260	/// Let us assume that M has a factor of d in it. Imposing this constraint on M
261	/// does not in any way hinder us from finding a value of M that is big enough.
262	/// Moreover, this function is only called when the symbolic part of the sample,
263	/// a_1*s_1 + ... + a_m*s_m, is known to be an integer.
264	///
265	/// Also, we can safely reduce the coefficients modulo d, so we have:
266	///
267	/// (b_1%d)y_1 + ... + (b_n%d)y_n = (-c%d) + k*d for some integer `k`
268	///
269	/// Note that all coefficient modulos here are non-negative. Also, all the
270	/// unknowns are non-negative here as both constraints and variables are
271	/// non-negative in LexSimplexBase. (We used the big M trick to make the
272	/// variables non-negative). Therefore, the LHS here is non-negative.
273	/// Since 0 <= (-c%d) < d, k is the quotient of dividing the LHS by d and
274	/// is therefore non-negative as well.
275	///
276	/// So we have
277	/// ((b_1%d)y_1 + ... + (b_n%d)y_n - (-c%d))/d >= 0.
278	///
279	/// The constraint is violated when added (it would be useless otherwise)
280	/// so we immediately try to move it to a column.
281	LogicalResult LexSimplexBase::addCut(unsigned row) {
282	DynamicAPInt d = tableau(row, 0);
283	unsigned cutRow = addZeroRow(/*makeRestricted=*/true);
284	tableau(cutRow, 0) = d;
285	tableau(cutRow, 1) = -mod(LHS: -tableau(row, 1), RHS: d); // -c%d.
286	tableau(cutRow, 2) = 0;
287	for (unsigned col = 3 + nSymbol, e = getNumColumns(); col < e; ++col)
288	tableau(cutRow, col) = mod(LHS: tableau(row, col), RHS: d); // b_i%d.
289	return moveRowUnknownToColumn(row: cutRow);
290	}
291
292	std::optional<unsigned> LexSimplex::maybeGetNonIntegralVarRow() const {
293	for (const Unknown &u : var) {
294	if (u.orientation == Orientation::Column)
295	continue;
296	// If the sample value is of the form (a/d)M + b/d, we need b to be
297	// divisible by d. We assume M contains all possible
298	// factors and is divisible by everything.
299	unsigned row = u.pos;
300	if (tableau(row, 1) % tableau(row, 0) != 0)
301	return row;
302	}
303	return {};
304	}
305
306	MaybeOptimum<SmallVector<DynamicAPInt, 8>> LexSimplex::findIntegerLexMin() {
307	// We first try to make the tableau consistent.
308	if (restoreRationalConsistency().failed())
309	return OptimumKind::Empty;
310
311	// Then, if the sample value is integral, we are done.
312	while (std::optional<unsigned> maybeRow = maybeGetNonIntegralVarRow()) {
313	// Otherwise, for the variable whose row has a non-integral sample value,
314	// we add a cut, a constraint that remove this rational point
315	// while preserving all integer points, thus keeping the lexmin the same.
316	// We then again try to make the tableau with the new constraint
317	// consistent. This continues until the tableau becomes empty, in which
318	// case there is no integer point, or until there are no variables with
319	// non-integral sample values.
320	//
321	// Failure indicates that the tableau became empty, which occurs when the
322	// polytope is integer empty.
323	if (addCut(row: *maybeRow).failed())
324	return OptimumKind::Empty;
325	if (restoreRationalConsistency().failed())
326	return OptimumKind::Empty;
327	}
328
329	MaybeOptimum<SmallVector<Fraction, 8>> sample = getRationalSample();
330	assert(!sample.isEmpty() && "If we reached here the sample should exist!");
331	if (sample.isUnbounded())
332	return OptimumKind::Unbounded;
333	return llvm::to_vector<8>(
334	Range: llvm::map_range(C&: *sample, F: std::mem_fn(pm: &Fraction::getAsInteger)));
335	}
336
337	bool LexSimplex::isSeparateInequality(ArrayRef<DynamicAPInt> coeffs) {
338	SimplexRollbackScopeExit scopeExit(*this);
339	addInequality(coeffs);
340	return findIntegerLexMin().isEmpty();
341	}
342
343	bool LexSimplex::isRedundantInequality(ArrayRef<DynamicAPInt> coeffs) {
344	return isSeparateInequality(coeffs: getComplementIneq(ineq: coeffs));
345	}
346
347	SmallVector<DynamicAPInt, 8>
348	SymbolicLexSimplex::getSymbolicSampleNumerator(unsigned row) const {
349	SmallVector<DynamicAPInt, 8> sample;
350	sample.reserve(N: nSymbol + 1);
351	for (unsigned col = 3; col < 3 + nSymbol; ++col)
352	sample.emplace_back(Args: tableau(row, col));
353	sample.emplace_back(Args: tableau(row, 1));
354	return sample;
355	}
356
357	SmallVector<DynamicAPInt, 8>
358	SymbolicLexSimplex::getSymbolicSampleIneq(unsigned row) const {
359	SmallVector<DynamicAPInt, 8> sample = getSymbolicSampleNumerator(row);
360	// The inequality is equivalent to the GCD-normalized one.
361	normalizeRange(range: sample);
362	return sample;
363	}
364
365	void LexSimplexBase::appendSymbol() {
366	appendVariable();
367	swapColumns(i: 3 + nSymbol, j: getNumColumns() - 1);
368	var.back().isSymbol = true;
369	nSymbol++;
370	}
371
372	static bool isRangeDivisibleBy(ArrayRef<DynamicAPInt> range,
373	const DynamicAPInt &divisor) {
374	assert(divisor > 0 && "divisor must be positive!");
375	return llvm::all_of(
376	Range&: range, P: [divisor](const DynamicAPInt &x) { return x % divisor == 0; });
377	}
378
379	bool SymbolicLexSimplex::isSymbolicSampleIntegral(unsigned row) const {
380	DynamicAPInt denom = tableau(row, 0);
381	return tableau(row, 1) % denom == 0 &&
382	isRangeDivisibleBy(range: tableau.getRow(row).slice(N: 3, M: nSymbol), divisor: denom);
383	}
384
385	/// This proceeds similarly to LexSimplexBase::addCut(). We are given a row that
386	/// has a symbolic sample value with fractional coefficients.
387	///
388	/// Let the row be
389	/// (c + coeffM*M + sum_i a_i*s_i + sum_j b_j*y_j)/d,
390	/// where s_1, ... s_m are the symbols and
391	/// y_1, ... y_n are the other basis unknowns.
392	///
393	/// As in LexSimplex::addCut, for this to be an integer, we want
394	///
395	/// coeffM*M + sum_j b_j*y_j = -c + sum_i (-a_i*s_i) (mod d)
396	///
397	/// This time, a_1*s_1 + ... + a_m*s_m may not be an integer. We find that
398	///
399	/// sum_i (b_i%d)y_i = ((-c%d) + sum_i (-a_i%d)s_i)%d + k*d for some integer k
400	///
401	/// where we take a modulo of the whole symbolic expression on the right to
402	/// bring it into the range [0, d - 1]. Therefore, as in addCut(),
403	/// k is the quotient on dividing the LHS by d, and since LHS >= 0, we have
404	/// k >= 0 as well. If all the a_i are divisible by d, then we can add the
405	/// constraint directly. Otherwise, we realize the modulo of the symbolic
406	/// expression by adding a division variable
407	///
408	/// q = ((-c%d) + sum_i (-a_i%d)s_i)/d
409	///
410	/// to the symbol domain, so the equality becomes
411	///
412	/// sum_i (b_i%d)y_i = (-c%d) + sum_i (-a_i%d)s_i - q*d + k*d for some integer k
413	///
414	/// So the cut is
415	/// (sum_i (b_i%d)y_i - (-c%d) - sum_i (-a_i%d)s_i + q*d)/d >= 0
416	/// This constraint is violated when added so we immediately try to move it to a
417	/// column.
418	LogicalResult SymbolicLexSimplex::addSymbolicCut(unsigned row) {
419	DynamicAPInt d = tableau(row, 0);
420	if (isRangeDivisibleBy(range: tableau.getRow(row).slice(N: 3, M: nSymbol), divisor: d)) {
421	// The coefficients of symbols in the symbol numerator are divisible
422	// by the denominator, so we can add the constraint directly,
423	// i.e., ignore the symbols and add a regular cut as in addCut().
424	return addCut(row);
425	}
426
427	// Construct the division variable `q = ((-c%d) + sum_i (-a_i%d)s_i)/d`.
428	SmallVector<DynamicAPInt, 8> divCoeffs;
429	divCoeffs.reserve(N: nSymbol + 1);
430	DynamicAPInt divDenom = d;
431	for (unsigned col = 3; col < 3 + nSymbol; ++col)
432	divCoeffs.emplace_back(Args: mod(LHS: -tableau(row, col), RHS: divDenom)); // (-a_i%d)s_i
433	divCoeffs.emplace_back(Args: mod(LHS: -tableau(row, 1), RHS: divDenom)); // -c%d.
434	normalizeDiv(num: divCoeffs, denom&: divDenom);
435
436	domainSimplex.addDivisionVariable(coeffs: divCoeffs, denom: divDenom);
437	domainPoly.addLocalFloorDiv(dividend: divCoeffs, divisor: divDenom);
438
439	// Update `this` to account for the additional symbol we just added.
440	appendSymbol();
441
442	// Add the cut (sum_i (b_i%d)y_i - (-c%d) + sum_i -(-a_i%d)s_i + q*d)/d >= 0.
443	unsigned cutRow = addZeroRow(/*makeRestricted=*/true);
444	tableau(cutRow, 0) = d;
445	tableau(cutRow, 2) = 0;
446
447	tableau(cutRow, 1) = -mod(LHS: -tableau(row, 1), RHS: d); // -(-c%d).
448	for (unsigned col = 3; col < 3 + nSymbol - 1; ++col)
449	tableau(cutRow, col) = -mod(LHS: -tableau(row, col), RHS: d); // -(-a_i%d)s_i.
450	tableau(cutRow, 3 + nSymbol - 1) = d; // q*d.
451
452	for (unsigned col = 3 + nSymbol, e = getNumColumns(); col < e; ++col)
453	tableau(cutRow, col) = mod(LHS: tableau(row, col), RHS: d); // (b_i%d)y_i.
454	return moveRowUnknownToColumn(row: cutRow);
455	}
456
457	void SymbolicLexSimplex::recordOutput(SymbolicLexOpt &result) const {
458	IntMatrix output(0, domainPoly.getNumVars() + 1);
459	output.reserveRows(rows: result.lexopt.getNumOutputs());
460	for (const Unknown &u : var) {
461	if (u.isSymbol)
462	continue;
463
464	if (u.orientation == Orientation::Column) {
465	// M + u has a sample value of zero so u has a sample value of -M, i.e,
466	// unbounded.
467	result.unboundedDomain.unionInPlace(disjunct: domainPoly);
468	return;
469	}
470
471	DynamicAPInt denom = tableau(u.pos, 0);
472	if (tableau(u.pos, 2) < denom) {
473	// M + u has a sample value of fM + something, where f < 1, so
474	// u = (f - 1)M + something, which has a negative coefficient for M,
475	// and so is unbounded.
476	result.unboundedDomain.unionInPlace(disjunct: domainPoly);
477	return;
478	}
479	assert(tableau(u.pos, 2) == denom &&
480	"Coefficient of M should not be greater than 1!");
481
482	SmallVector<DynamicAPInt, 8> sample = getSymbolicSampleNumerator(row: u.pos);
483	for (DynamicAPInt &elem : sample) {
484	assert(elem % denom == 0 && "coefficients must be integral!");
485	elem /= denom;
486	}
487	output.appendExtraRow(elems: sample);
488	}
489
490	// Store the output in a MultiAffineFunction and add it the result.
491	PresburgerSpace funcSpace = result.lexopt.getSpace();
492	funcSpace.insertVar(kind: VarKind::Local, pos: 0, num: domainPoly.getNumLocalVars());
493
494	result.lexopt.addPiece(
495	piece: {.domain: PresburgerSet(domainPoly),
496	.output: MultiAffineFunction(funcSpace, output, domainPoly.getLocalReprs())});
497	}
498
499	std::optional<unsigned> SymbolicLexSimplex::maybeGetAlwaysViolatedRow() {
500	// First look for rows that are clearly violated just from the big M
501	// coefficient, without needing to perform any simplex queries on the domain.
502	for (unsigned row = 0, e = getNumRows(); row < e; ++row)
503	if (tableau(row, 2) < 0)
504	return row;
505
506	for (unsigned row = 0, e = getNumRows(); row < e; ++row) {
507	if (tableau(row, 2) > 0)
508	continue;
509	if (domainSimplex.isSeparateInequality(coeffs: getSymbolicSampleIneq(row))) {
510	// Sample numerator always takes negative values in the symbol domain.
511	return row;
512	}
513	}
514	return {};
515	}
516
517	std::optional<unsigned> SymbolicLexSimplex::maybeGetNonIntegralVarRow() {
518	for (const Unknown &u : var) {
519	if (u.orientation == Orientation::Column)
520	continue;
521	assert(!u.isSymbol && "Symbol should not be in row orientation!");
522	if (!isSymbolicSampleIntegral(row: u.pos))
523	return u.pos;
524	}
525	return {};
526	}
527
528	/// The non-branching pivots are just the ones moving the rows
529	/// that are always violated in the symbol domain.
530	LogicalResult SymbolicLexSimplex::doNonBranchingPivots() {
531	while (std::optional<unsigned> row = maybeGetAlwaysViolatedRow())
532	if (moveRowUnknownToColumn(row: *row).failed())
533	return failure();
534	return success();
535	}
536
537	SymbolicLexOpt SymbolicLexSimplex::computeSymbolicIntegerLexMin() {
538	SymbolicLexOpt result(PresburgerSpace::getRelationSpace(
539	/*numDomain=*/domainPoly.getNumDimVars(),
540	/*numRange=*/var.size() - nSymbol,
541	/*numSymbols=*/domainPoly.getNumSymbolVars()));
542
543	/// The algorithm is more naturally expressed recursively, but we implement
544	/// it iteratively here to avoid potential issues with stack overflows in the
545	/// compiler. We explicitly maintain the stack frames in a vector.
546	///
547	/// To "recurse", we store the current "stack frame", i.e., state variables
548	/// that we will need when we "return", into `stack`, increment `level`, and
549	/// `continue`. To "tail recurse", we just `continue`.
550	/// To "return", we decrement `level` and `continue`.
551	///
552	/// When there is no stack frame for the current `level`, this indicates that
553	/// we have just "recursed" or "tail recursed". When there does exist one,
554	/// this indicates that we have just "returned" from recursing. There is only
555	/// one point at which non-tail calls occur so we always "return" there.
556	unsigned level = 1;
557	struct StackFrame {
558	int splitIndex;
559	unsigned snapshot;
560	unsigned domainSnapshot;
561	IntegerRelation::CountsSnapshot domainPolyCounts;
562	};
563	SmallVector<StackFrame, 8> stack;
564
565	while (level > 0) {
566	assert(level >= stack.size());
567	if (level > stack.size()) {
568	if (empty || domainSimplex.findIntegerLexMin().isEmpty()) {
569	// No integer points; return.
570	--level;
571	continue;
572	}
573
574	if (doNonBranchingPivots().failed()) {
575	// Could not find pivots for violated constraints; return.
576	--level;
577	continue;
578	}
579
580	SmallVector<DynamicAPInt, 8> symbolicSample;
581	unsigned splitRow = 0;
582	for (unsigned e = getNumRows(); splitRow < e; ++splitRow) {
583	if (tableau(splitRow, 2) > 0)
584	continue;
585	assert(tableau(splitRow, 2) == 0 &&
586	"Non-branching pivots should have been handled already!");
587
588	symbolicSample = getSymbolicSampleIneq(row: splitRow);
589	if (domainSimplex.isRedundantInequality(coeffs: symbolicSample))
590	continue;
591
592	// It's neither redundant nor separate, so it takes both positive and
593	// negative values, and hence constitutes a row for which we need to
594	// split the domain and separately run each case.
595	assert(!domainSimplex.isSeparateInequality(symbolicSample) &&
596	"Non-branching pivots should have been handled already!");
597	break;
598	}
599
600	if (splitRow < getNumRows()) {
601	unsigned domainSnapshot = domainSimplex.getSnapshot();
602	IntegerRelation::CountsSnapshot domainPolyCounts =
603	domainPoly.getCounts();
604
605	// First, we consider the part of the domain where the row is not
606	// violated. We don't have to do any pivots for the row in this case,
607	// but we record the additional constraint that defines this part of
608	// the domain.
609	domainSimplex.addInequality(coeffs: symbolicSample);
610	domainPoly.addInequality(inEq: symbolicSample);
611
612	// Recurse.
613	//
614	// On return, the basis as a set is preserved but not the internal
615	// ordering within rows or columns. Thus, we take note of the index of
616	// the Unknown that caused the split, which may be in a different
617	// row when we come back from recursing. We will need this to recurse
618	// on the other part of the split domain, where the row is violated.
619	//
620	// Note that we have to capture the index above and not a reference to
621	// the Unknown itself, since the array it lives in might get
622	// reallocated.
623	int splitIndex = rowUnknown[splitRow];
624	unsigned snapshot = getSnapshot();
625	stack.emplace_back(
626	Args: StackFrame{.splitIndex: splitIndex, .snapshot: snapshot, .domainSnapshot: domainSnapshot, .domainPolyCounts: domainPolyCounts});
627	++level;
628	continue;
629	}
630
631	// The tableau is rationally consistent for the current domain.
632	// Now we look for non-integral sample values and add cuts for them.
633	if (std::optional<unsigned> row = maybeGetNonIntegralVarRow()) {
634	if (addSymbolicCut(row: *row).failed()) {
635	// No integral points; return.
636	--level;
637	continue;
638	}
639
640	// Rerun this level with the added cut constraint (tail recurse).
641	continue;
642	}
643
644	// Record output and return.
645	recordOutput(result);
646	--level;
647	continue;
648	}
649
650	if (level == stack.size()) {
651	// We have "returned" from "recursing".
652	const StackFrame &frame = stack.back();
653	domainPoly.truncate(counts: frame.domainPolyCounts);
654	domainSimplex.rollback(snapshot: frame.domainSnapshot);
655	rollback(snapshot: frame.snapshot);
656	const Unknown &u = unknownFromIndex(index: frame.splitIndex);
657
658	// Drop the frame. We don't need it anymore.
659	stack.pop_back();
660
661	// Now we consider the part of the domain where the unknown `splitIndex`
662	// was negative.
663	assert(u.orientation == Orientation::Row &&
664	"The split row should have been returned to row orientation!");
665	SmallVector<DynamicAPInt, 8> splitIneq =
666	getComplementIneq(ineq: getSymbolicSampleIneq(row: u.pos));
667	normalizeRange(range: splitIneq);
668	if (moveRowUnknownToColumn(row: u.pos).failed()) {
669	// The unknown can't be made non-negative; return.
670	--level;
671	continue;
672	}
673
674	// The unknown can be made negative; recurse with the corresponding domain
675	// constraints.
676	domainSimplex.addInequality(coeffs: splitIneq);
677	domainPoly.addInequality(inEq: splitIneq);
678
679	// We are now taking care of the second half of the domain and we don't
680	// need to do anything else here after returning, so it's a tail recurse.
681	continue;
682	}
683	}
684
685	return result;
686	}
687
688	bool LexSimplex::rowIsViolated(unsigned row) const {
689	if (tableau(row, 2) < 0)
690	return true;
691	if (tableau(row, 2) == 0 && tableau(row, 1) < 0)
692	return true;
693	return false;
694	}
695
696	std::optional<unsigned> LexSimplex::maybeGetViolatedRow() const {
697	for (unsigned row = 0, e = getNumRows(); row < e; ++row)
698	if (rowIsViolated(row))
699	return row;
700	return {};
701	}
702
703	/// We simply look for violated rows and keep trying to move them to column
704	/// orientation, which always succeeds unless the constraints have no solution
705	/// in which case we just give up and return.
706	LogicalResult LexSimplex::restoreRationalConsistency() {
707	if (empty)
708	return failure();
709	while (std::optional<unsigned> maybeViolatedRow = maybeGetViolatedRow())
710	if (moveRowUnknownToColumn(row: *maybeViolatedRow).failed())
711	return failure();
712	return success();
713	}
714
715	// Move the row unknown to column orientation while preserving lexicopositivity
716	// of the basis transform. The sample value of the row must be non-positive.
717	//
718	// We only consider pivots where the pivot element is positive. Suppose no such
719	// pivot exists, i.e., some violated row has no positive coefficient for any
720	// basis unknown. The row can be represented as (s + c_1*u_1 + ... + c_n*u_n)/d,
721	// where d is the denominator, s is the sample value and the c_i are the basis
722	// coefficients. If s != 0, then since any feasible assignment of the basis
723	// satisfies u_i >= 0 for all i, and we have s < 0 as well as c_i < 0 for all i,
724	// any feasible assignment would violate this row and therefore the constraints
725	// have no solution.
726	//
727	// We can preserve lexicopositivity by picking the pivot column with positive
728	// pivot element that makes the lexicographically smallest change to the sample
729	// point.
730	//
731	// Proof. Let
732	// x = (x_1, ... x_n) be the variables,
733	// z = (z_1, ... z_m) be the constraints,
734	// y = (y_1, ... y_n) be the current basis, and
735	// define w = (x_1, ... x_n, z_1, ... z_m) = B*y + s.
736	// B is basically the simplex tableau of our implementation except that instead
737	// of only describing the transform to get back the non-basis unknowns, it
738	// defines the values of all the unknowns in terms of the basis unknowns.
739	// Similarly, s is the column for the sample value.
740	//
741	// Our goal is to show that each column in B, restricted to the first n
742	// rows, is lexicopositive after the pivot if it is so before. This is
743	// equivalent to saying the columns in the whole matrix are lexicopositive;
744	// there must be some non-zero element in every column in the first n rows since
745	// the n variables cannot be spanned without using all the n basis unknowns.
746	//
747	// Consider a pivot where z_i replaces y_j in the basis. Recall the pivot
748	// transform for the tableau derived for SimplexBase::pivot:
749	//
750	// pivot col other col pivot col other col
751	// pivot row a b -> pivot row 1/a -b/a
752	// other row c d other row c/a d - bc/a
753	//
754	// Similarly, a pivot results in B changing to B' and c to c'; the difference
755	// between the tableau and these matrices B and B' is that there is no special
756	// case for the pivot row, since it continues to represent the same unknown. The
757	// same formula applies for all rows:
758	//
759	// B'.col(j) = B.col(j) / B(i,j)
760	// B'.col(k) = B.col(k) - B(i,k) * B.col(j) / B(i,j) for k != j
761	// and similarly, s' = s - s_i * B.col(j) / B(i,j).
762	//
763	// If s_i == 0, then the sample value remains unchanged. Otherwise, if s_i < 0,
764	// the change in sample value when pivoting with column a is lexicographically
765	// smaller than that when pivoting with column b iff B.col(a) / B(i, a) is
766	// lexicographically smaller than B.col(b) / B(i, b).
767	//
768	// Since B(i, j) > 0, column j remains lexicopositive.
769	//
770	// For the other columns, suppose C.col(k) is not lexicopositive.
771	// This means that for some p, for all t < p,
772	// C(t,k) = 0 => B(t,k) = B(t,j) * B(i,k) / B(i,j) and
773	// C(t,k) < 0 => B(p,k) < B(t,j) * B(i,k) / B(i,j),
774	// which is in contradiction to the fact that B.col(j) / B(i,j) must be
775	// lexicographically smaller than B.col(k) / B(i,k), since it lexicographically
776	// minimizes the change in sample value.
777	LogicalResult LexSimplexBase::moveRowUnknownToColumn(unsigned row) {
778	std::optional<unsigned> maybeColumn;
779	for (unsigned col = 3 + nSymbol, e = getNumColumns(); col < e; ++col) {
780	if (tableau(row, col) <= 0)
781	continue;
782	maybeColumn =
783	!maybeColumn ? col : getLexMinPivotColumn(row, colA: *maybeColumn, colB: col);
784	}
785
786	if (!maybeColumn)
787	return failure();
788
789	pivot(row, col: *maybeColumn);
790	return success();
791	}
792
793	unsigned LexSimplexBase::getLexMinPivotColumn(unsigned row, unsigned colA,
794	unsigned colB) const {
795	// First, let's consider the non-symbolic case.
796	// A pivot causes the following change. (in the diagram the matrix elements
797	// are shown as rationals and there is no common denominator used)
798	//
799	// pivot col big M col const col
800	// pivot row a p b
801	// other row c q d
802	// |
803	// v
804	//
805	// pivot col big M col const col
806	// pivot row 1/a -p/a -b/a
807	// other row c/a q - pc/a d - bc/a
808	//
809	// Let the sample value of the pivot row be s = pM + b before the pivot. Since
810	// the pivot row represents a violated constraint we know that s < 0.
811	//
812	// If the variable is a non-pivot column, its sample value is zero before and
813	// after the pivot.
814	//
815	// If the variable is the pivot column, then its sample value goes from 0 to
816	// (-p/a)M + (-b/a), i.e. 0 to -(pM + b)/a. Thus the change in the sample
817	// value is -s/a.
818	//
819	// If the variable is the pivot row, its sample value goes from s to 0, for a
820	// change of -s.
821	//
822	// If the variable is a non-pivot row, its sample value changes from
823	// qM + d to qM + d + (-pc/a)M + (-bc/a). Thus the change in sample value
824	// is -(pM + b)(c/a) = -sc/a.
825	//
826	// Thus the change in sample value is either 0, -s/a, -s, or -sc/a. Here -s is
827	// fixed for all calls to this function since the row and tableau are fixed.
828	// The callee just wants to compare the return values with the return value of
829	// other invocations of the same function. So the -s is common for all
830	// comparisons involved and can be ignored, since -s is strictly positive.
831	//
832	// Thus we take away this common factor and just return 0, 1/a, 1, or c/a as
833	// appropriate. This allows us to run the entire algorithm treating M
834	// symbolically, as the pivot to be performed does not depend on the value
835	// of M, so long as the sample value s is negative. Note that this is not
836	// because of any special feature of M; by the same argument, we ignore the
837	// symbols too. The caller ensure that the sample value s is negative for
838	// all possible values of the symbols.
839	auto getSampleChangeCoeffForVar = [this, row](unsigned col,
840	const Unknown &u) -> Fraction {
841	DynamicAPInt a = tableau(row, col);
842	if (u.orientation == Orientation::Column) {
843	// Pivot column case.
844	if (u.pos == col)
845	return {1, a};
846
847	// Non-pivot column case.
848	return {0, 1};
849	}
850
851	// Pivot row case.
852	if (u.pos == row)
853	return {1, 1};
854
855	// Non-pivot row case.
856	DynamicAPInt c = tableau(u.pos, col);
857	return {c, a};
858	};
859
860	for (const Unknown &u : var) {
861	Fraction changeA = getSampleChangeCoeffForVar(colA, u);
862	Fraction changeB = getSampleChangeCoeffForVar(colB, u);
863	if (changeA < changeB)
864	return colA;
865	if (changeA > changeB)
866	return colB;
867	}
868
869	// If we reached here, both result in exactly the same changes, so it
870	// doesn't matter which we return.
871	return colA;
872	}
873
874	/// Find a pivot to change the sample value of the row in the specified
875	/// direction. The returned pivot row will involve `row` if and only if the
876	/// unknown is unbounded in the specified direction.
877	///
878	/// To increase (resp. decrease) the value of a row, we need to find a live
879	/// column with a non-zero coefficient. If the coefficient is positive, we need
880	/// to increase (decrease) the value of the column, and if the coefficient is
881	/// negative, we need to decrease (increase) the value of the column. Also,
882	/// we cannot decrease the sample value of restricted columns.
883	///
884	/// If multiple columns are valid, we break ties by considering a lexicographic
885	/// ordering where we prefer unknowns with lower index.
886	std::optional<SimplexBase::Pivot>
887	Simplex::findPivot(int row, Direction direction) const {
888	std::optional<unsigned> col;
889	for (unsigned j = 2, e = getNumColumns(); j < e; ++j) {
890	DynamicAPInt elem = tableau(row, j);
891	if (elem == 0)
892	continue;
893
894	if (unknownFromColumn(col: j).restricted &&
895	!signMatchesDirection(elem, direction))
896	continue;
897	if (!col || colUnknown[j] < colUnknown[*col])
898	col = j;
899	}
900
901	if (!col)
902	return {};
903
904	Direction newDirection =
905	tableau(row, *col) < 0 ? flippedDirection(direction) : direction;
906	std::optional<unsigned> maybePivotRow = findPivotRow(skipRow: row, direction: newDirection, col: *col);
907	return Pivot{.row: maybePivotRow.value_or(u&: row), .column: *col};
908	}
909
910	/// Swap the associated unknowns for the row and the column.
911	///
912	/// First we swap the index associated with the row and column. Then we update
913	/// the unknowns to reflect their new position and orientation.
914	void SimplexBase::swapRowWithCol(unsigned row, unsigned col) {
915	std::swap(a&: rowUnknown[row], b&: colUnknown[col]);
916	Unknown &uCol = unknownFromColumn(col);
917	Unknown &uRow = unknownFromRow(row);
918	uCol.orientation = Orientation::Column;
919	uRow.orientation = Orientation::Row;
920	uCol.pos = col;
921	uRow.pos = row;
922	}
923
924	void SimplexBase::pivot(Pivot pair) { pivot(row: pair.row, col: pair.column); }
925
926	/// Pivot pivotRow and pivotCol.
927	///
928	/// Let R be the pivot row unknown and let C be the pivot col unknown.
929	/// Since initially R = a*C + sum b_i * X_i
930	/// (where the sum is over the other column's unknowns, x_i)
931	/// C = (R - (sum b_i * X_i))/a
932	///
933	/// Let u be some other row unknown.
934	/// u = c*C + sum d_i * X_i
935	/// So u = c*(R - sum b_i * X_i)/a + sum d_i * X_i
936	///
937	/// This results in the following transform:
938	/// pivot col other col pivot col other col
939	/// pivot row a b -> pivot row 1/a -b/a
940	/// other row c d other row c/a d - bc/a
941	///
942	/// Taking into account the common denominators p and q:
943	///
944	/// pivot col other col pivot col other col
945	/// pivot row a/p b/p -> pivot row p/a -b/a
946	/// other row c/q d/q other row cp/aq (da - bc)/aq
947	///
948	/// The pivot row transform is accomplished be swapping a with the pivot row's
949	/// common denominator and negating the pivot row except for the pivot column
950	/// element.
951	void SimplexBase::pivot(unsigned pivotRow, unsigned pivotCol) {
952	assert(pivotCol >= getNumFixedCols() && "Refusing to pivot invalid column");
953	assert(!unknownFromColumn(pivotCol).isSymbol);
954
955	swapRowWithCol(row: pivotRow, col: pivotCol);
956	std::swap(a&: tableau(pivotRow, 0), b&: tableau(pivotRow, pivotCol));
957	// We need to negate the whole pivot row except for the pivot column.
958	if (tableau(pivotRow, 0) < 0) {
959	// If the denominator is negative, we negate the row by simply negating the
960	// denominator.
961	tableau(pivotRow, 0) = -tableau(pivotRow, 0);
962	tableau(pivotRow, pivotCol) = -tableau(pivotRow, pivotCol);
963	} else {
964	for (unsigned col = 1, e = getNumColumns(); col < e; ++col) {
965	if (col == pivotCol)
966	continue;
967	tableau(pivotRow, col) = -tableau(pivotRow, col);
968	}
969	}
970	tableau.normalizeRow(row: pivotRow);
971
972	for (unsigned row = 0, numRows = getNumRows(); row < numRows; ++row) {
973	if (row == pivotRow)
974	continue;
975	if (tableau(row, pivotCol) == 0) // Nothing to do.
976	continue;
977	tableau(row, 0) *= tableau(pivotRow, 0);
978	for (unsigned col = 1, numCols = getNumColumns(); col < numCols; ++col) {
979	if (col == pivotCol)
980	continue;
981	// Add rather than subtract because the pivot row has been negated.
982	tableau(row, col) = tableau(row, col) * tableau(pivotRow, 0) +
983	tableau(row, pivotCol) * tableau(pivotRow, col);
984	}
985	tableau(row, pivotCol) *= tableau(pivotRow, pivotCol);
986	tableau.normalizeRow(row);
987	}
988	}
989
990	/// Perform pivots until the unknown has a non-negative sample value or until
991	/// no more upward pivots can be performed. Return success if we were able to
992	/// bring the row to a non-negative sample value, and failure otherwise.
993	LogicalResult Simplex::restoreRow(Unknown &u) {
994	assert(u.orientation == Orientation::Row &&
995	"unknown should be in row position");
996
997	while (tableau(u.pos, 1) < 0) {
998	std::optional<Pivot> maybePivot = findPivot(row: u.pos, direction: Direction::Up);
999	if (!maybePivot)
1000	break;
1001
1002	pivot(pair: *maybePivot);
1003	if (u.orientation == Orientation::Column)
1004	return success(); // the unknown is unbounded above.
1005	}
1006	return success(IsSuccess: tableau(u.pos, 1) >= 0);
1007	}
1008
1009	/// Find a row that can be used to pivot the column in the specified direction.
1010	/// This returns an empty optional if and only if the column is unbounded in the
1011	/// specified direction (ignoring skipRow, if skipRow is set).
1012	///
1013	/// If skipRow is set, this row is not considered, and (if it is restricted) its
1014	/// restriction may be violated by the returned pivot. Usually, skipRow is set
1015	/// because we don't want to move it to column position unless it is unbounded,
1016	/// and we are either trying to increase the value of skipRow or explicitly
1017	/// trying to make skipRow negative, so we are not concerned about this.
1018	///
1019	/// If the direction is up (resp. down) and a restricted row has a negative
1020	/// (positive) coefficient for the column, then this row imposes a bound on how
1021	/// much the sample value of the column can change. Such a row with constant
1022	/// term c and coefficient f for the column imposes a bound of c/|f| on the
1023	/// change in sample value (in the specified direction). (note that c is
1024	/// non-negative here since the row is restricted and the tableau is consistent)
1025	///
1026	/// We iterate through the rows and pick the row which imposes the most
1027	/// stringent bound, since pivoting with a row changes the row's sample value to
1028	/// 0 and hence saturates the bound it imposes. We break ties between rows that
1029	/// impose the same bound by considering a lexicographic ordering where we
1030	/// prefer unknowns with lower index value.
1031	std::optional<unsigned> Simplex::findPivotRow(std::optional<unsigned> skipRow,
1032	Direction direction,
1033	unsigned col) const {
1034	std::optional<unsigned> retRow;
1035	// Initialize these to zero in order to silence a warning about retElem and
1036	// retConst being used uninitialized in the initialization of `diff` below. In
1037	// reality, these are always initialized when that line is reached since these
1038	// are set whenever retRow is set.
1039	DynamicAPInt retElem, retConst;
1040	for (unsigned row = nRedundant, e = getNumRows(); row < e; ++row) {
1041	if (skipRow && row == *skipRow)
1042	continue;
1043	DynamicAPInt elem = tableau(row, col);
1044	if (elem == 0)
1045	continue;
1046	if (!unknownFromRow(row).restricted)
1047	continue;
1048	if (signMatchesDirection(elem, direction))
1049	continue;
1050	DynamicAPInt constTerm = tableau(row, 1);
1051
1052	if (!retRow) {
1053	retRow = row;
1054	retElem = elem;
1055	retConst = constTerm;
1056	continue;
1057	}
1058
1059	DynamicAPInt diff = retConst * elem - constTerm * retElem;
1060	if ((diff == 0 && rowUnknown[row] < rowUnknown[*retRow]) ||
1061	(diff != 0 && !signMatchesDirection(elem: diff, direction))) {
1062	retRow = row;
1063	retElem = elem;
1064	retConst = constTerm;
1065	}
1066	}
1067	return retRow;
1068	}
1069
1070	bool SimplexBase::isEmpty() const { return empty; }
1071
1072	void SimplexBase::swapRows(unsigned i, unsigned j) {
1073	if (i == j)
1074	return;
1075	tableau.swapRows(row: i, otherRow: j);
1076	std::swap(a&: rowUnknown[i], b&: rowUnknown[j]);
1077	unknownFromRow(row: i).pos = i;
1078	unknownFromRow(row: j).pos = j;
1079	}
1080
1081	void SimplexBase::swapColumns(unsigned i, unsigned j) {
1082	assert(i < getNumColumns() && j < getNumColumns() &&
1083	"Invalid columns provided!");
1084	if (i == j)
1085	return;
1086	tableau.swapColumns(column: i, otherColumn: j);
1087	std::swap(a&: colUnknown[i], b&: colUnknown[j]);
1088	unknownFromColumn(col: i).pos = i;
1089	unknownFromColumn(col: j).pos = j;
1090	}
1091
1092	/// Mark this tableau empty and push an entry to the undo stack.
1093	void SimplexBase::markEmpty() {
1094	// If the set is already empty, then we shouldn't add another UnmarkEmpty log
1095	// entry, since in that case the Simplex will be erroneously marked as
1096	// non-empty when rolling back past this point.
1097	if (empty)
1098	return;
1099	undoLog.emplace_back(Args: UndoLogEntry::UnmarkEmpty);
1100	empty = true;
1101	}
1102
1103	/// Add an inequality to the tableau. If coeffs is c_0, c_1, ... c_n, where n
1104	/// is the current number of variables, then the corresponding inequality is
1105	/// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} >= 0.
1106	///
1107	/// We add the inequality and mark it as restricted. We then try to make its
1108	/// sample value non-negative. If this is not possible, the tableau has become
1109	/// empty and we mark it as such.
1110	void Simplex::addInequality(ArrayRef<DynamicAPInt> coeffs) {
1111	unsigned conIndex = addRow(coeffs, /*makeRestricted=*/true);
1112	LogicalResult result = restoreRow(u&: con[conIndex]);
1113	if (result.failed())
1114	markEmpty();
1115	}
1116
1117	/// Add an equality to the tableau. If coeffs is c_0, c_1, ... c_n, where n
1118	/// is the current number of variables, then the corresponding equality is
1119	/// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} == 0.
1120	///
1121	/// We simply add two opposing inequalities, which force the expression to
1122	/// be zero.
1123	void SimplexBase::addEquality(ArrayRef<DynamicAPInt> coeffs) {
1124	addInequality(coeffs);
1125	SmallVector<DynamicAPInt, 8> negatedCoeffs;
1126	negatedCoeffs.reserve(N: coeffs.size());
1127	for (const DynamicAPInt &coeff : coeffs)
1128	negatedCoeffs.emplace_back(Args: -coeff);
1129	addInequality(coeffs: negatedCoeffs);
1130	}
1131
1132	unsigned SimplexBase::getNumVariables() const { return var.size(); }
1133	unsigned SimplexBase::getNumConstraints() const { return con.size(); }
1134
1135	/// Return a snapshot of the current state. This is just the current size of the
1136	/// undo log.
1137	unsigned SimplexBase::getSnapshot() const { return undoLog.size(); }
1138
1139	unsigned SimplexBase::getSnapshotBasis() {
1140	SmallVector<int, 8> basis;
1141	basis.reserve(N: colUnknown.size());
1142	for (int index : colUnknown) {
1143	if (index != nullIndex)
1144	basis.emplace_back(Args&: index);
1145	}
1146	savedBases.emplace_back(Args: std::move(basis));
1147
1148	undoLog.emplace_back(Args: UndoLogEntry::RestoreBasis);
1149	return undoLog.size() - 1;
1150	}
1151
1152	void SimplexBase::removeLastConstraintRowOrientation() {
1153	assert(con.back().orientation == Orientation::Row);
1154
1155	// Move this unknown to the last row and remove the last row from the
1156	// tableau.
1157	swapRows(i: con.back().pos, j: getNumRows() - 1);
1158	// It is not strictly necessary to shrink the tableau, but for now we
1159	// maintain the invariant that the tableau has exactly getNumRows()
1160	// rows.
1161	tableau.resizeVertically(newNRows: getNumRows() - 1);
1162	rowUnknown.pop_back();
1163	con.pop_back();
1164	}
1165
1166	// This doesn't find a pivot row only if the column has zero
1167	// coefficients for every row.
1168	//
1169	// If the unknown is a constraint, this can't happen, since it was added
1170	// initially as a row. Such a row could never have been pivoted to a column. So
1171	// a pivot row will always be found if we have a constraint.
1172	//
1173	// If we have a variable, then the column has zero coefficients for every row
1174	// iff no constraints have been added with a non-zero coefficient for this row.
1175	std::optional<unsigned> SimplexBase::findAnyPivotRow(unsigned col) {
1176	for (unsigned row = nRedundant, e = getNumRows(); row < e; ++row)
1177	if (tableau(row, col) != 0)
1178	return row;
1179	return {};
1180	}
1181
1182	// It's not valid to remove the constraint by deleting the column since this
1183	// would result in an invalid basis.
1184	void Simplex::undoLastConstraint() {
1185	if (con.back().orientation == Orientation::Column) {
1186	// We try to find any pivot row for this column that preserves tableau
1187	// consistency (except possibly the column itself, which is going to be
1188	// deallocated anyway).
1189	//
1190	// If no pivot row is found in either direction, then the unknown is
1191	// unbounded in both directions and we are free to perform any pivot at
1192	// all. To do this, we just need to find any row with a non-zero
1193	// coefficient for the column. findAnyPivotRow will always be able to
1194	// find such a row for a constraint.
1195	unsigned column = con.back().pos;
1196	if (std::optional<unsigned> maybeRow =
1197	findPivotRow(skipRow: {}, direction: Direction::Up, col: column)) {
1198	pivot(pivotRow: *maybeRow, pivotCol: column);
1199	} else if (std::optional<unsigned> maybeRow =
1200	findPivotRow(skipRow: {}, direction: Direction::Down, col: column)) {
1201	pivot(pivotRow: *maybeRow, pivotCol: column);
1202	} else {
1203	std::optional<unsigned> row = findAnyPivotRow(col: column);
1204	assert(row && "Pivot should always exist for a constraint!");
1205	pivot(pivotRow: *row, pivotCol: column);
1206	}
1207	}
1208	removeLastConstraintRowOrientation();
1209	}
1210
1211	// It's not valid to remove the constraint by deleting the column since this
1212	// would result in an invalid basis.
1213	void LexSimplexBase::undoLastConstraint() {
1214	if (con.back().orientation == Orientation::Column) {
1215	// When removing the last constraint during a rollback, we just need to find
1216	// any pivot at all, i.e., any row with non-zero coefficient for the
1217	// column, because when rolling back a lexicographic simplex, we always
1218	// end by restoring the exact basis that was present at the time of the
1219	// snapshot, so what pivots we perform while undoing doesn't matter as
1220	// long as we get the unknown to row orientation and remove it.
1221	unsigned column = con.back().pos;
1222	std::optional<unsigned> row = findAnyPivotRow(col: column);
1223	assert(row && "Pivot should always exist for a constraint!");
1224	pivot(pivotRow: *row, pivotCol: column);
1225	}
1226	removeLastConstraintRowOrientation();
1227	}
1228
1229	void SimplexBase::undo(UndoLogEntry entry) {
1230	if (entry == UndoLogEntry::RemoveLastConstraint) {
1231	// Simplex and LexSimplex handle this differently, so we call out to a
1232	// virtual function to handle this.
1233	undoLastConstraint();
1234	} else if (entry == UndoLogEntry::RemoveLastVariable) {
1235	// Whenever we are rolling back the addition of a variable, it is guaranteed
1236	// that the variable will be in column position.
1237	//
1238	// We can see this as follows: any constraint that depends on this variable
1239	// was added after this variable was added, so the addition of such
1240	// constraints should already have been rolled back by the time we get to
1241	// rolling back the addition of the variable. Therefore, no constraint
1242	// currently has a component along the variable, so the variable itself must
1243	// be part of the basis.
1244	assert(var.back().orientation == Orientation::Column &&
1245	"Variable to be removed must be in column orientation!");
1246
1247	if (var.back().isSymbol)
1248	nSymbol--;
1249
1250	// Move this variable to the last column and remove the column from the
1251	// tableau.
1252	swapColumns(i: var.back().pos, j: getNumColumns() - 1);
1253	tableau.resizeHorizontally(newNColumns: getNumColumns() - 1);
1254	var.pop_back();
1255	colUnknown.pop_back();
1256	} else if (entry == UndoLogEntry::UnmarkEmpty) {
1257	empty = false;
1258	} else if (entry == UndoLogEntry::UnmarkLastRedundant) {
1259	nRedundant--;
1260	} else if (entry == UndoLogEntry::RestoreBasis) {
1261	assert(!savedBases.empty() && "No bases saved!");
1262
1263	SmallVector<int, 8> basis = std::move(savedBases.back());
1264	savedBases.pop_back();
1265
1266	for (int index : basis) {
1267	Unknown &u = unknownFromIndex(index);
1268	if (u.orientation == Orientation::Column)
1269	continue;
1270	for (unsigned col = getNumFixedCols(), e = getNumColumns(); col < e;
1271	col++) {
1272	assert(colUnknown[col] != nullIndex &&
1273	"Column should not be a fixed column!");
1274	if (llvm::is_contained(Range&: basis, Element: colUnknown[col]))
1275	continue;
1276	if (tableau(u.pos, col) == 0)
1277	continue;
1278	pivot(pivotRow: u.pos, pivotCol: col);
1279	break;
1280	}
1281
1282	assert(u.orientation == Orientation::Column && "No pivot found!");
1283	}
1284	}
1285	}
1286
1287	/// Rollback to the specified snapshot.
1288	///
1289	/// We undo all the log entries until the log size when the snapshot was taken
1290	/// is reached.
1291	void SimplexBase::rollback(unsigned snapshot) {
1292	while (undoLog.size() > snapshot) {
1293	undo(entry: undoLog.back());
1294	undoLog.pop_back();
1295	}
1296	}
1297
1298	/// We add the usual floor division constraints:
1299	/// `0 <= coeffs - denom*q <= denom - 1`, where `q` is the new division
1300	/// variable.
1301	///
1302	/// This constrains the remainder `coeffs - denom*q` to be in the
1303	/// range `[0, denom - 1]`, which fixes the integer value of the quotient `q`.
1304	void SimplexBase::addDivisionVariable(ArrayRef<DynamicAPInt> coeffs,
1305	const DynamicAPInt &denom) {
1306	assert(denom > 0 && "Denominator must be positive!");
1307	appendVariable();
1308
1309	SmallVector<DynamicAPInt, 8> ineq(coeffs);
1310	DynamicAPInt constTerm = ineq.back();
1311	ineq.back() = -denom;
1312	ineq.emplace_back(Args&: constTerm);
1313	addInequality(coeffs: ineq);
1314
1315	for (DynamicAPInt &coeff : ineq)
1316	coeff = -coeff;
1317	ineq.back() += denom - 1;
1318	addInequality(coeffs: ineq);
1319	}
1320
1321	void SimplexBase::appendVariable(unsigned count) {
1322	if (count == 0)
1323	return;
1324	var.reserve(N: var.size() + count);
1325	colUnknown.reserve(N: colUnknown.size() + count);
1326	for (unsigned i = 0; i < count; ++i) {
1327	var.emplace_back(Args: Orientation::Column, /*restricted=*/Args: false,
1328	/*pos=*/Args: getNumColumns() + i);
1329	colUnknown.emplace_back(Args: var.size() - 1);
1330	}
1331	tableau.resizeHorizontally(newNColumns: getNumColumns() + count);
1332	undoLog.insert(I: undoLog.end(), NumToInsert: count, Elt: UndoLogEntry::RemoveLastVariable);
1333	}
1334
1335	/// Add all the constraints from the given IntegerRelation.
1336	void SimplexBase::intersectIntegerRelation(const IntegerRelation &rel) {
1337	assert(rel.getNumVars() == getNumVariables() &&
1338	"IntegerRelation must have same dimensionality as simplex");
1339	for (unsigned i = 0, e = rel.getNumInequalities(); i < e; ++i)
1340	addInequality(coeffs: rel.getInequality(idx: i));
1341	for (unsigned i = 0, e = rel.getNumEqualities(); i < e; ++i)
1342	addEquality(coeffs: rel.getEquality(idx: i));
1343	}
1344
1345	MaybeOptimum<Fraction> Simplex::computeRowOptimum(Direction direction,
1346	unsigned row) {
1347	// Keep trying to find a pivot for the row in the specified direction.
1348	while (std::optional<Pivot> maybePivot = findPivot(row, direction)) {
1349	// If findPivot returns a pivot involving the row itself, then the optimum
1350	// is unbounded, so we return std::nullopt.
1351	if (maybePivot->row == row)
1352	return OptimumKind::Unbounded;
1353	pivot(pair: *maybePivot);
1354	}
1355
1356	// The row has reached its optimal sample value, which we return.
1357	// The sample value is the entry in the constant column divided by the common
1358	// denominator for this row.
1359	return Fraction(tableau(row, 1), tableau(row, 0));
1360	}
1361
1362	/// Compute the optimum of the specified expression in the specified direction,
1363	/// or std::nullopt if it is unbounded.
1364	MaybeOptimum<Fraction> Simplex::computeOptimum(Direction direction,
1365	ArrayRef<DynamicAPInt> coeffs) {
1366	if (empty)
1367	return OptimumKind::Empty;
1368
1369	SimplexRollbackScopeExit scopeExit(*this);
1370	unsigned conIndex = addRow(coeffs);
1371	unsigned row = con[conIndex].pos;
1372	return computeRowOptimum(direction, row);
1373	}
1374
1375	MaybeOptimum<Fraction> Simplex::computeOptimum(Direction direction,
1376	Unknown &u) {
1377	if (empty)
1378	return OptimumKind::Empty;
1379	if (u.orientation == Orientation::Column) {
1380	unsigned column = u.pos;
1381	std::optional<unsigned> pivotRow = findPivotRow(skipRow: {}, direction, col: column);
1382	// If no pivot is returned, the constraint is unbounded in the specified
1383	// direction.
1384	if (!pivotRow)
1385	return OptimumKind::Unbounded;
1386	pivot(pivotRow: *pivotRow, pivotCol: column);
1387	}
1388
1389	unsigned row = u.pos;
1390	MaybeOptimum<Fraction> optimum = computeRowOptimum(direction, row);
1391	if (u.restricted && direction == Direction::Down &&
1392	(optimum.isUnbounded() || *optimum < Fraction(0, 1))) {
1393	if (restoreRow(u).failed())
1394	llvm_unreachable("Could not restore row!");
1395	}
1396	return optimum;
1397	}
1398
1399	bool Simplex::isBoundedAlongConstraint(unsigned constraintIndex) {
1400	assert(!empty && "It is not meaningful to ask whether a direction is bounded "
1401	"in an empty set.");
1402	// The constraint's perpendicular is already bounded below, since it is a
1403	// constraint. If it is also bounded above, we can return true.
1404	return computeOptimum(direction: Direction::Up, u&: con[constraintIndex]).isBounded();
1405	}
1406
1407	/// Redundant constraints are those that are in row orientation and lie in
1408	/// rows 0 to nRedundant - 1.
1409	bool Simplex::isMarkedRedundant(unsigned constraintIndex) const {
1410	const Unknown &u = con[constraintIndex];
1411	return u.orientation == Orientation::Row && u.pos < nRedundant;
1412	}
1413
1414	/// Mark the specified row redundant.
1415	///
1416	/// This is done by moving the unknown to the end of the block of redundant
1417	/// rows (namely, to row nRedundant) and incrementing nRedundant to
1418	/// accomodate the new redundant row.
1419	void Simplex::markRowRedundant(Unknown &u) {
1420	assert(u.orientation == Orientation::Row &&
1421	"Unknown should be in row position!");
1422	assert(u.pos >= nRedundant && "Unknown is already marked redundant!");
1423	swapRows(i: u.pos, j: nRedundant);
1424	++nRedundant;
1425	undoLog.emplace_back(Args: UndoLogEntry::UnmarkLastRedundant);
1426	}
1427
1428	/// Find a subset of constraints that is redundant and mark them redundant.
1429	void Simplex::detectRedundant(unsigned offset, unsigned count) {
1430	assert(offset + count <= con.size() && "invalid range!");
1431	// It is not meaningful to talk about redundancy for empty sets.
1432	if (empty)
1433	return;
1434
1435	// Iterate through the constraints and check for each one if it can attain
1436	// negative sample values. If it can, it's not redundant. Otherwise, it is.
1437	// We mark redundant constraints redundant.
1438	//
1439	// Constraints that get marked redundant in one iteration are not respected
1440	// when checking constraints in later iterations. This prevents, for example,
1441	// two identical constraints both being marked redundant since each is
1442	// redundant given the other one. In this example, only the first of the
1443	// constraints that is processed will get marked redundant, as it should be.
1444	for (unsigned i = 0; i < count; ++i) {
1445	Unknown &u = con[offset + i];
1446	if (u.orientation == Orientation::Column) {
1447	unsigned column = u.pos;
1448	std::optional<unsigned> pivotRow =
1449	findPivotRow(skipRow: {}, direction: Direction::Down, col: column);
1450	// If no downward pivot is returned, the constraint is unbounded below
1451	// and hence not redundant.
1452	if (!pivotRow)
1453	continue;
1454	pivot(pivotRow: *pivotRow, pivotCol: column);
1455	}
1456
1457	unsigned row = u.pos;
1458	MaybeOptimum<Fraction> minimum = computeRowOptimum(direction: Direction::Down, row);
1459	if (minimum.isUnbounded() || *minimum < Fraction(0, 1)) {
1460	// Constraint is unbounded below or can attain negative sample values and
1461	// hence is not redundant.
1462	if (restoreRow(u).failed())
1463	llvm_unreachable("Could not restore non-redundant row!");
1464	continue;
1465	}
1466
1467	markRowRedundant(u);
1468	}
1469	}
1470
1471	bool Simplex::isUnbounded() {
1472	if (empty)
1473	return false;
1474
1475	SmallVector<DynamicAPInt, 8> dir(var.size() + 1);
1476	for (unsigned i = 0; i < var.size(); ++i) {
1477	dir[i] = 1;
1478
1479	if (computeOptimum(direction: Direction::Up, coeffs: dir).isUnbounded())
1480	return true;
1481
1482	if (computeOptimum(direction: Direction::Down, coeffs: dir).isUnbounded())
1483	return true;
1484
1485	dir[i] = 0;
1486	}
1487	return false;
1488	}
1489
1490	/// Make a tableau to represent a pair of points in the original tableau.
1491	///
1492	/// The product constraints and variables are stored as: first A's, then B's.
1493	///
1494	/// The product tableau has row layout:
1495	/// A's redundant rows, B's redundant rows, A's other rows, B's other rows.
1496	///
1497	/// It has column layout:
1498	/// denominator, constant, A's columns, B's columns.
1499	Simplex Simplex::makeProduct(const Simplex &a, const Simplex &b) {
1500	unsigned numVar = a.getNumVariables() + b.getNumVariables();
1501	unsigned numCon = a.getNumConstraints() + b.getNumConstraints();
1502	Simplex result(numVar);
1503
1504	result.tableau.reserveRows(rows: numCon);
1505	result.empty = a.empty || b.empty;
1506
1507	auto concat = [](ArrayRef<Unknown> v, ArrayRef<Unknown> w) {
1508	SmallVector<Unknown, 8> result;
1509	result.reserve(N: v.size() + w.size());
1510	llvm::append_range(C&: result, R&: v);
1511	llvm::append_range(C&: result, R&: w);
1512	return result;
1513	};
1514	result.con = concat(a.con, b.con);
1515	result.var = concat(a.var, b.var);
1516
1517	auto indexFromBIndex = [&](int index) {
1518	return index >= 0 ? a.getNumVariables() + index
1519	: ~(a.getNumConstraints() + ~index);
1520	};
1521
1522	result.colUnknown.assign(NumElts: 2, Elt: nullIndex);
1523	for (unsigned i = 2, e = a.getNumColumns(); i < e; ++i) {
1524	result.colUnknown.emplace_back(Args: a.colUnknown[i]);
1525	result.unknownFromIndex(index: result.colUnknown.back()).pos =
1526	result.colUnknown.size() - 1;
1527	}
1528	for (unsigned i = 2, e = b.getNumColumns(); i < e; ++i) {
1529	result.colUnknown.emplace_back(Args: indexFromBIndex(b.colUnknown[i]));
1530	result.unknownFromIndex(index: result.colUnknown.back()).pos =
1531	result.colUnknown.size() - 1;
1532	}
1533
1534	auto appendRowFromA = [&](unsigned row) {
1535	unsigned resultRow = result.tableau.appendExtraRow();
1536	for (unsigned col = 0, e = a.getNumColumns(); col < e; ++col)
1537	result.tableau(resultRow, col) = a.tableau(row, col);
1538	result.rowUnknown.emplace_back(Args: a.rowUnknown[row]);
1539	result.unknownFromIndex(index: result.rowUnknown.back()).pos =
1540	result.rowUnknown.size() - 1;
1541	};
1542
1543	// Also fixes the corresponding entry in rowUnknown and var/con (as the case
1544	// may be).
1545	auto appendRowFromB = [&](unsigned row) {
1546	unsigned resultRow = result.tableau.appendExtraRow();
1547	result.tableau(resultRow, 0) = b.tableau(row, 0);
1548	result.tableau(resultRow, 1) = b.tableau(row, 1);
1549
1550	unsigned offset = a.getNumColumns() - 2;
1551	for (unsigned col = 2, e = b.getNumColumns(); col < e; ++col)
1552	result.tableau(resultRow, offset + col) = b.tableau(row, col);
1553	result.rowUnknown.emplace_back(Args: indexFromBIndex(b.rowUnknown[row]));
1554	result.unknownFromIndex(index: result.rowUnknown.back()).pos =
1555	result.rowUnknown.size() - 1;
1556	};
1557
1558	result.nRedundant = a.nRedundant + b.nRedundant;
1559	for (unsigned row = 0; row < a.nRedundant; ++row)
1560	appendRowFromA(row);
1561	for (unsigned row = 0; row < b.nRedundant; ++row)
1562	appendRowFromB(row);
1563	for (unsigned row = a.nRedundant, e = a.getNumRows(); row < e; ++row)
1564	appendRowFromA(row);
1565	for (unsigned row = b.nRedundant, e = b.getNumRows(); row < e; ++row)
1566	appendRowFromB(row);
1567
1568	return result;
1569	}
1570
1571	std::optional<SmallVector<Fraction, 8>> Simplex::getRationalSample() const {
1572	if (empty)
1573	return {};
1574
1575	SmallVector<Fraction, 8> sample;
1576	sample.reserve(N: var.size());
1577	// Push the sample value for each variable into the vector.
1578	for (const Unknown &u : var) {
1579	if (u.orientation == Orientation::Column) {
1580	// If the variable is in column position, its sample value is zero.
1581	sample.emplace_back(Args: 0, Args: 1);
1582	} else {
1583	// If the variable is in row position, its sample value is the
1584	// entry in the constant column divided by the denominator.
1585	DynamicAPInt denom = tableau(u.pos, 0);
1586	sample.emplace_back(Args: tableau(u.pos, 1), Args&: denom);
1587	}
1588	}
1589	return sample;
1590	}
1591
1592	void LexSimplexBase::addInequality(ArrayRef<DynamicAPInt> coeffs) {
1593	addRow(coeffs, /*makeRestricted=*/true);
1594	}
1595
1596	MaybeOptimum<SmallVector<Fraction, 8>> LexSimplex::getRationalSample() const {
1597	if (empty)
1598	return OptimumKind::Empty;
1599
1600	SmallVector<Fraction, 8> sample;
1601	sample.reserve(N: var.size());
1602	// Push the sample value for each variable into the vector.
1603	for (const Unknown &u : var) {
1604	// When the big M parameter is being used, each variable x is represented
1605	// as M + x, so its sample value is finite if and only if it is of the
1606	// form 1*M + c. If the coefficient of M is not one then the sample value
1607	// is infinite, and we return an empty optional.
1608
1609	if (u.orientation == Orientation::Column) {
1610	// If the variable is in column position, the sample value of M + x is
1611	// zero, so x = -M which is unbounded.
1612	return OptimumKind::Unbounded;
1613	}
1614
1615	// If the variable is in row position, its sample value is the
1616	// entry in the constant column divided by the denominator.
1617	DynamicAPInt denom = tableau(u.pos, 0);
1618	if (usingBigM)
1619	if (tableau(u.pos, 2) != denom)
1620	return OptimumKind::Unbounded;
1621	sample.emplace_back(Args: tableau(u.pos, 1), Args&: denom);
1622	}
1623	return sample;
1624	}
1625
1626	std::optional<SmallVector<DynamicAPInt, 8>>
1627	Simplex::getSamplePointIfIntegral() const {
1628	// If the tableau is empty, no sample point exists.
1629	if (empty)
1630	return {};
1631
1632	// The value will always exist since the Simplex is non-empty.
1633	SmallVector<Fraction, 8> rationalSample = *getRationalSample();
1634	SmallVector<DynamicAPInt, 8> integerSample;
1635	integerSample.reserve(N: var.size());
1636	for (const Fraction &coord : rationalSample) {
1637	// If the sample is non-integral, return std::nullopt.
1638	if (coord.num % coord.den != 0)
1639	return {};
1640	integerSample.emplace_back(Args: coord.num / coord.den);
1641	}
1642	return integerSample;
1643	}
1644
1645	/// Given a simplex for a polytope, construct a new simplex whose variables are
1646	/// identified with a pair of points (x, y) in the original polytope. Supports
1647	/// some operations needed for generalized basis reduction. In what follows,
1648	/// dotProduct(x, y) = x_1 * y_1 + x_2 * y_2 + ... x_n * y_n where n is the
1649	/// dimension of the original polytope.
1650	///
1651	/// This supports adding equality constraints dotProduct(dir, x - y) == 0. It
1652	/// also supports rolling back this addition, by maintaining a snapshot stack
1653	/// that contains a snapshot of the Simplex's state for each equality, just
1654	/// before that equality was added.
1655	class presburger::GBRSimplex {
1656	using Orientation = Simplex::Orientation;
1657
1658	public:
1659	GBRSimplex(const Simplex &originalSimplex)
1660	: simplex(Simplex::makeProduct(a: originalSimplex, b: originalSimplex)),
1661	simplexConstraintOffset(simplex.getNumConstraints()) {}
1662
1663	/// Add an equality dotProduct(dir, x - y) == 0.
1664	/// First pushes a snapshot for the current simplex state to the stack so
1665	/// that this can be rolled back later.
1666	void addEqualityForDirection(ArrayRef<DynamicAPInt> dir) {
1667	assert(llvm::any_of(dir, [](const DynamicAPInt &x) { return x != 0; }) &&
1668	"Direction passed is the zero vector!");
1669	snapshotStack.emplace_back(Args: simplex.getSnapshot());
1670	simplex.addEquality(coeffs: getCoeffsForDirection(dir));
1671	}
1672	/// Compute max(dotProduct(dir, x - y)).
1673	Fraction computeWidth(ArrayRef<DynamicAPInt> dir) {
1674	MaybeOptimum<Fraction> maybeWidth =
1675	simplex.computeOptimum(direction: Direction::Up, coeffs: getCoeffsForDirection(dir));
1676	assert(maybeWidth.isBounded() && "Width should be bounded!");
1677	return *maybeWidth;
1678	}
1679
1680	/// Compute max(dotProduct(dir, x - y)) and save the dual variables for only
1681	/// the direction equalities to `dual`.
1682	Fraction computeWidthAndDuals(ArrayRef<DynamicAPInt> dir,
1683	SmallVectorImpl<DynamicAPInt> &dual,
1684	DynamicAPInt &dualDenom) {
1685	// We can't just call into computeWidth or computeOptimum since we need to
1686	// access the state of the tableau after computing the optimum, and these
1687	// functions rollback the insertion of the objective function into the
1688	// tableau before returning. We instead add a row for the objective function
1689	// ourselves, call into computeOptimum, compute the duals from the tableau
1690	// state, and finally rollback the addition of the row before returning.
1691	SimplexRollbackScopeExit scopeExit(simplex);
1692	unsigned conIndex = simplex.addRow(coeffs: getCoeffsForDirection(dir));
1693	unsigned row = simplex.con[conIndex].pos;
1694	MaybeOptimum<Fraction> maybeWidth =
1695	simplex.computeRowOptimum(direction: Simplex::Direction::Up, row);
1696	assert(maybeWidth.isBounded() && "Width should be bounded!");
1697	dualDenom = simplex.tableau(row, 0);
1698	dual.clear();
1699	dual.reserve(N: (conIndex - simplexConstraintOffset) / 2);
1700
1701	// The increment is i += 2 because equalities are added as two inequalities,
1702	// one positive and one negative. Each iteration processes one equality.
1703	for (unsigned i = simplexConstraintOffset; i < conIndex; i += 2) {
1704	// The dual variable for an inequality in column orientation is the
1705	// negative of its coefficient at the objective row. If the inequality is
1706	// in row orientation, the corresponding dual variable is zero.
1707	//
1708	// We want the dual for the original equality, which corresponds to two
1709	// inequalities: a positive inequality, which has the same coefficients as
1710	// the equality, and a negative equality, which has negated coefficients.
1711	//
1712	// Note that at most one of these inequalities can be in column
1713	// orientation because the column unknowns should form a basis and hence
1714	// must be linearly independent. If the positive inequality is in column
1715	// position, its dual is the dual corresponding to the equality. If the
1716	// negative inequality is in column position, the negation of its dual is
1717	// the dual corresponding to the equality. If neither is in column
1718	// position, then that means that this equality is redundant, and its dual
1719	// is zero.
1720	//
1721	// Note that it is NOT valid to perform pivots during the computation of
1722	// the duals. This entire dual computation must be performed on the same
1723	// tableau configuration.
1724	assert((simplex.con[i].orientation != Orientation::Column ||
1725	simplex.con[i + 1].orientation != Orientation::Column) &&
1726	"Both inequalities for the equality cannot be in column "
1727	"orientation!");
1728	if (simplex.con[i].orientation == Orientation::Column)
1729	dual.emplace_back(Args: -simplex.tableau(row, simplex.con[i].pos));
1730	else if (simplex.con[i + 1].orientation == Orientation::Column)
1731	dual.emplace_back(Args&: simplex.tableau(row, simplex.con[i + 1].pos));
1732	else
1733	dual.emplace_back(Args: 0);
1734	}
1735	return *maybeWidth;
1736	}
1737
1738	/// Remove the last equality that was added through addEqualityForDirection.
1739	///
1740	/// We do this by rolling back to the snapshot at the top of the stack, which
1741	/// should be a snapshot taken just before the last equality was added.
1742	void removeLastEquality() {
1743	assert(!snapshotStack.empty() && "Snapshot stack is empty!");
1744	simplex.rollback(snapshot: snapshotStack.back());
1745	snapshotStack.pop_back();
1746	}
1747
1748	private:
1749	/// Returns coefficients of the expression 'dot_product(dir, x - y)',
1750	/// i.e., dir_1 * x_1 + dir_2 * x_2 + ... + dir_n * x_n
1751	/// - dir_1 * y_1 - dir_2 * y_2 - ... - dir_n * y_n,
1752	/// where n is the dimension of the original polytope.
1753	SmallVector<DynamicAPInt, 8>
1754	getCoeffsForDirection(ArrayRef<DynamicAPInt> dir) {
1755	assert(2 * dir.size() == simplex.getNumVariables() &&
1756	"Direction vector has wrong dimensionality");
1757	SmallVector<DynamicAPInt, 8> coeffs(dir);
1758	coeffs.reserve(N: dir.size() + 1);
1759	for (const DynamicAPInt &coeff : dir)
1760	coeffs.emplace_back(Args: -coeff);
1761	coeffs.emplace_back(Args: 0); // constant term
1762	return coeffs;
1763	}
1764
1765	Simplex simplex;
1766	/// The first index of the equality constraints, the index immediately after
1767	/// the last constraint in the initial product simplex.
1768	unsigned simplexConstraintOffset;
1769	/// A stack of snapshots, used for rolling back.
1770	SmallVector<unsigned, 8> snapshotStack;
1771	};
1772
1773	/// Reduce the basis to try and find a direction in which the polytope is
1774	/// "thin". This only works for bounded polytopes.
1775	///
1776	/// This is an implementation of the algorithm described in the paper
1777	/// "An Implementation of Generalized Basis Reduction for Integer Programming"
1778	/// by W. Cook, T. Rutherford, H. E. Scarf, D. Shallcross.
1779	///
1780	/// Let b_{level}, b_{level + 1}, ... b_n be the current basis.
1781	/// Let width_i(v) = max <v, x - y> where x and y are points in the original
1782	/// polytope such that <b_j, x - y> = 0 is satisfied for all level <= j < i.
1783	///
1784	/// In every iteration, we first replace b_{i+1} with b_{i+1} + u*b_i, where u
1785	/// is the integer such that width_i(b_{i+1} + u*b_i) is minimized. Let dual_i
1786	/// be the dual variable associated with the constraint <b_i, x - y> = 0 when
1787	/// computing width_{i+1}(b_{i+1}). It can be shown that dual_i is the
1788	/// minimizing value of u, if it were allowed to be fractional. Due to
1789	/// convexity, the minimizing integer value is either floor(dual_i) or
1790	/// ceil(dual_i), so we just need to check which of these gives a lower
1791	/// width_{i+1} value. If dual_i turned out to be an integer, then u = dual_i.
1792	///
1793	/// Now if width_i(b_{i+1}) < 0.75 * width_i(b_i), we swap b_i and (the new)
1794	/// b_{i + 1} and decrement i (unless i = level, in which case we stay at the
1795	/// same i). Otherwise, we increment i.
1796	///
1797	/// We keep f values and duals cached and invalidate them when necessary.
1798	/// Whenever possible, we use them instead of recomputing them. We implement the
1799	/// algorithm as follows.
1800	///
1801	/// In an iteration at i we need to compute:
1802	/// a) width_i(b_{i + 1})
1803	/// b) width_i(b_i)
1804	/// c) the integer u that minimizes width_i(b_{i + 1} + u*b_i)
1805	///
1806	/// If width_i(b_i) is not already cached, we compute it.
1807	///
1808	/// If the duals are not already cached, we compute width_{i+1}(b_{i+1}) and
1809	/// store the duals from this computation.
1810	///
1811	/// We call updateBasisWithUAndGetFCandidate, which finds the minimizing value
1812	/// of u as explained before, caches the duals from this computation, sets
1813	/// b_{i+1} to b_{i+1} + u*b_i, and returns the new value of width_i(b_{i+1}).
1814	///
1815	/// Now if width_i(b_{i+1}) < 0.75 * width_i(b_i), we swap b_i and b_{i+1} and
1816	/// decrement i, resulting in the basis
1817	/// ... b_{i - 1}, b_{i + 1} + u*b_i, b_i, b_{i+2}, ...
1818	/// with corresponding f values
1819	/// ... width_{i-1}(b_{i-1}), width_i(b_{i+1} + u*b_i), width_{i+1}(b_i), ...
1820	/// The values up to i - 1 remain unchanged. We have just gotten the middle
1821	/// value from updateBasisWithUAndGetFCandidate, so we can update that in the
1822	/// cache. The value at width_{i+1}(b_i) is unknown, so we evict this value from
1823	/// the cache. The iteration after decrementing needs exactly the duals from the
1824	/// computation of width_i(b_{i + 1} + u*b_i), so we keep these in the cache.
1825	///
1826	/// When incrementing i, no cached f values get invalidated. However, the cached
1827	/// duals do get invalidated as the duals for the higher levels are different.
1828	void Simplex::reduceBasis(IntMatrix &basis, unsigned level) {
1829	const Fraction epsilon(3, 4);
1830
1831	if (level == basis.getNumRows() - 1)
1832	return;
1833
1834	GBRSimplex gbrSimplex(*this);
1835	SmallVector<Fraction, 8> width;
1836	SmallVector<DynamicAPInt, 8> dual;
1837	DynamicAPInt dualDenom;
1838
1839	// Finds the value of u that minimizes width_i(b_{i+1} + u*b_i), caches the
1840	// duals from this computation, sets b_{i+1} to b_{i+1} + u*b_i, and returns
1841	// the new value of width_i(b_{i+1}).
1842	//
1843	// If dual_i is not an integer, the minimizing value must be either
1844	// floor(dual_i) or ceil(dual_i). We compute the expression for both and
1845	// choose the minimizing value.
1846	//
1847	// If dual_i is an integer, we don't need to perform these computations. We
1848	// know that in this case,
1849	// a) u = dual_i.
1850	// b) one can show that dual_j for j < i are the same duals we would have
1851	// gotten from computing width_i(b_{i + 1} + u*b_i), so the correct duals
1852	// are the ones already in the cache.
1853	// c) width_i(b_{i+1} + u*b_i) = min_{alpha} width_i(b_{i+1} + alpha * b_i),
1854	// which
1855	// one can show is equal to width_{i+1}(b_{i+1}). The latter value must
1856	// be in the cache, so we get it from there and return it.
1857	auto updateBasisWithUAndGetFCandidate = [&](unsigned i) -> Fraction {
1858	assert(i < level + dual.size() && "dual_i is not known!");
1859
1860	DynamicAPInt u = floorDiv(LHS: dual[i - level], RHS: dualDenom);
1861	basis.addToRow(sourceRow: i, targetRow: i + 1, scale: u);
1862	if (dual[i - level] % dualDenom != 0) {
1863	SmallVector<DynamicAPInt, 8> candidateDual[2];
1864	DynamicAPInt candidateDualDenom[2];
1865	Fraction widthI[2];
1866
1867	// Initially u is floor(dual) and basis reflects this.
1868	widthI[0] = gbrSimplex.computeWidthAndDuals(
1869	dir: basis.getRow(row: i + 1), dual&: candidateDual[0], dualDenom&: candidateDualDenom[0]);
1870
1871	// Now try ceil(dual), i.e. floor(dual) + 1.
1872	++u;
1873	basis.addToRow(sourceRow: i, targetRow: i + 1, scale: 1);
1874	widthI[1] = gbrSimplex.computeWidthAndDuals(
1875	dir: basis.getRow(row: i + 1), dual&: candidateDual[1], dualDenom&: candidateDualDenom[1]);
1876
1877	unsigned j = widthI[0] < widthI[1] ? 0 : 1;
1878	if (j == 0)
1879	// Subtract 1 to go from u = ceil(dual) back to floor(dual).
1880	basis.addToRow(sourceRow: i, targetRow: i + 1, scale: -1);
1881
1882	// width_i(b{i+1} + u*b_i) should be minimized at our value of u.
1883	// We assert that this holds by checking that the values of width_i at
1884	// u - 1 and u + 1 are greater than or equal to the value at u. If the
1885	// width is lesser at either of the adjacent values, then our computed
1886	// value of u is clearly not the minimizer. Otherwise by convexity the
1887	// computed value of u is really the minimizer.
1888
1889	// Check the value at u - 1.
1890	assert(gbrSimplex.computeWidth(scaleAndAddForAssert(
1891	basis.getRow(i + 1), DynamicAPInt(-1), basis.getRow(i))) >=
1892	widthI[j] &&
1893	"Computed u value does not minimize the width!");
1894	// Check the value at u + 1.
1895	assert(gbrSimplex.computeWidth(scaleAndAddForAssert(
1896	basis.getRow(i + 1), DynamicAPInt(+1), basis.getRow(i))) >=
1897	widthI[j] &&
1898	"Computed u value does not minimize the width!");
1899
1900	dual = std::move(candidateDual[j]);
1901	dualDenom = candidateDualDenom[j];
1902	return widthI[j];
1903	}
1904
1905	assert(i + 1 - level < width.size() && "width_{i+1} wasn't saved");
1906	// f_i(b_{i+1} + dual*b_i) == width_{i+1}(b_{i+1}) when `dual` minimizes the
1907	// LHS. (note: the basis has already been updated, so b_{i+1} + dual*b_i in
1908	// the above expression is equal to basis.getRow(i+1) below.)
1909	assert(gbrSimplex.computeWidth(basis.getRow(i + 1)) ==
1910	width[i + 1 - level]);
1911	return width[i + 1 - level];
1912	};
1913
1914	// In the ith iteration of the loop, gbrSimplex has constraints for directions
1915	// from `level` to i - 1.
1916	unsigned i = level;
1917	while (i < basis.getNumRows() - 1) {
1918	if (i >= level + width.size()) {
1919	// We don't even know the value of f_i(b_i), so let's find that first.
1920	// We have to do this first since later we assume that width already
1921	// contains values up to and including i.
1922
1923	assert((i == 0 || i - 1 < level + width.size()) &&
1924	"We are at level i but we don't know the value of width_{i-1}");
1925
1926	// We don't actually use these duals at all, but it doesn't matter
1927	// because this case should only occur when i is level, and there are no
1928	// duals in that case anyway.
1929	assert(i == level && "This case should only occur when i == level");
1930	width.emplace_back(
1931	Args: gbrSimplex.computeWidthAndDuals(dir: basis.getRow(row: i), dual, dualDenom));
1932	}
1933
1934	if (i >= level + dual.size()) {
1935	assert(i + 1 >= level + width.size() &&
1936	"We don't know dual_i but we know width_{i+1}");
1937	// We don't know dual for our level, so let's find it.
1938	gbrSimplex.addEqualityForDirection(dir: basis.getRow(row: i));
1939	width.emplace_back(Args: gbrSimplex.computeWidthAndDuals(dir: basis.getRow(row: i + 1),
1940	dual, dualDenom));
1941	gbrSimplex.removeLastEquality();
1942	}
1943
1944	// This variable stores width_i(b_{i+1} + u*b_i).
1945	Fraction widthICandidate = updateBasisWithUAndGetFCandidate(i);
1946	if (widthICandidate < epsilon * width[i - level]) {
1947	basis.swapRows(row: i, otherRow: i + 1);
1948	width[i - level] = widthICandidate;
1949	// The values of width_{i+1}(b_{i+1}) and higher may change after the
1950	// swap, so we remove the cached values here.
1951	width.resize(N: i - level + 1);
1952	if (i == level) {
1953	dual.clear();
1954	continue;
1955	}
1956
1957	gbrSimplex.removeLastEquality();
1958	i--;
1959	continue;
1960	}
1961
1962	// Invalidate duals since the higher level needs to recompute its own duals.
1963	dual.clear();
1964	gbrSimplex.addEqualityForDirection(dir: basis.getRow(row: i));
1965	i++;
1966	}
1967	}
1968
1969	/// Search for an integer sample point using a branch and bound algorithm.
1970	///
1971	/// Each row in the basis matrix is a vector, and the set of basis vectors
1972	/// should span the space. Initially this is the identity matrix,
1973	/// i.e., the basis vectors are just the variables.
1974	///
1975	/// In every level, a value is assigned to the level-th basis vector, as
1976	/// follows. Compute the minimum and maximum rational values of this direction.
1977	/// If only one integer point lies in this range, constrain the variable to
1978	/// have this value and recurse to the next variable.
1979	///
1980	/// If the range has multiple values, perform generalized basis reduction via
1981	/// reduceBasis and then compute the bounds again. Now we try constraining
1982	/// this direction in the first value in this range and "recurse" to the next
1983	/// level. If we fail to find a sample, we try assigning the direction the next
1984	/// value in this range, and so on.
1985	///
1986	/// If no integer sample is found from any of the assignments, or if the range
1987	/// contains no integer value, then of course the polytope is empty for the
1988	/// current assignment of the values in previous levels, so we return to
1989	/// the previous level.
1990	///
1991	/// If we reach the last level where all the variables have been assigned values
1992	/// already, then we simply return the current sample point if it is integral,
1993	/// and go back to the previous level otherwise.
1994	///
1995	/// To avoid potentially arbitrarily large recursion depths leading to stack
1996	/// overflows, this algorithm is implemented iteratively.
1997	std::optional<SmallVector<DynamicAPInt, 8>> Simplex::findIntegerSample() {
1998	if (empty)
1999	return {};
2000
2001	unsigned nDims = var.size();
2002	IntMatrix basis = IntMatrix::identity(dimension: nDims);
2003
2004	unsigned level = 0;
2005	// The snapshot just before constraining a direction to a value at each level.
2006	SmallVector<unsigned, 8> snapshotStack;
2007	// The maximum value in the range of the direction for each level.
2008	SmallVector<DynamicAPInt, 8> upperBoundStack;
2009	// The next value to try constraining the basis vector to at each level.
2010	SmallVector<DynamicAPInt, 8> nextValueStack;
2011
2012	snapshotStack.reserve(N: basis.getNumRows());
2013	upperBoundStack.reserve(N: basis.getNumRows());
2014	nextValueStack.reserve(N: basis.getNumRows());
2015	while (level != -1u) {
2016	if (level == basis.getNumRows()) {
2017	// We've assigned values to all variables. Return if we have a sample,
2018	// or go back up to the previous level otherwise.
2019	if (auto maybeSample = getSamplePointIfIntegral())
2020	return maybeSample;
2021	level--;
2022	continue;
2023	}
2024
2025	if (level >= upperBoundStack.size()) {
2026	// We haven't populated the stack values for this level yet, so we have
2027	// just come down a level ("recursed"). Find the lower and upper bounds.
2028	// If there is more than one integer point in the range, perform
2029	// generalized basis reduction.
2030	SmallVector<DynamicAPInt, 8> basisCoeffs =
2031	llvm::to_vector<8>(Range: basis.getRow(row: level));
2032	basisCoeffs.emplace_back(Args: 0);
2033
2034	auto [minRoundedUp, maxRoundedDown] = computeIntegerBounds(coeffs: basisCoeffs);
2035
2036	// We don't have any integer values in the range.
2037	// Pop the stack and return up a level.
2038	if (minRoundedUp.isEmpty() || maxRoundedDown.isEmpty()) {
2039	assert((minRoundedUp.isEmpty() && maxRoundedDown.isEmpty()) &&
2040	"If one bound is empty, both should be.");
2041	snapshotStack.pop_back();
2042	nextValueStack.pop_back();
2043	upperBoundStack.pop_back();
2044	level--;
2045	continue;
2046	}
2047
2048	// We already checked the empty case above.
2049	assert((minRoundedUp.isBounded() && maxRoundedDown.isBounded()) &&
2050	"Polyhedron should be bounded!");
2051
2052	// Heuristic: if the sample point is integral at this point, just return
2053	// it.
2054	if (auto maybeSample = getSamplePointIfIntegral())
2055	return *maybeSample;
2056
2057	if (*minRoundedUp < *maxRoundedDown) {
2058	reduceBasis(basis, level);
2059	basisCoeffs = llvm::to_vector<8>(Range: basis.getRow(row: level));
2060	basisCoeffs.emplace_back(Args: 0);
2061	std::tie(args&: minRoundedUp, args&: maxRoundedDown) =
2062	computeIntegerBounds(coeffs: basisCoeffs);
2063	}
2064
2065	snapshotStack.emplace_back(Args: getSnapshot());
2066	// The smallest value in the range is the next value to try.
2067	// The values in the optionals are guaranteed to exist since we know the
2068	// polytope is bounded.
2069	nextValueStack.emplace_back(Args&: *minRoundedUp);
2070	upperBoundStack.emplace_back(Args&: *maxRoundedDown);
2071	}
2072
2073	assert((snapshotStack.size() - 1 == level &&
2074	nextValueStack.size() - 1 == level &&
2075	upperBoundStack.size() - 1 == level) &&
2076	"Mismatched variable stack sizes!");
2077
2078	// Whether we "recursed" or "returned" from a lower level, we rollback
2079	// to the snapshot of the starting state at this level. (in the "recursed"
2080	// case this has no effect)
2081	rollback(snapshot: snapshotStack.back());
2082	DynamicAPInt nextValue = nextValueStack.back();
2083	++nextValueStack.back();
2084	if (nextValue > upperBoundStack.back()) {
2085	// We have exhausted the range and found no solution. Pop the stack and
2086	// return up a level.
2087	snapshotStack.pop_back();
2088	nextValueStack.pop_back();
2089	upperBoundStack.pop_back();
2090	level--;
2091	continue;
2092	}
2093
2094	// Try the next value in the range and "recurse" into the next level.
2095	SmallVector<DynamicAPInt, 8> basisCoeffs(basis.getRow(row: level).begin(),
2096	basis.getRow(row: level).end());
2097	basisCoeffs.emplace_back(Args: -nextValue);
2098	addEquality(coeffs: basisCoeffs);
2099	level++;
2100	}
2101
2102	return {};
2103	}
2104
2105	/// Compute the minimum and maximum integer values the expression can take. We
2106	/// compute each separately.
2107	std::pair<MaybeOptimum<DynamicAPInt>, MaybeOptimum<DynamicAPInt>>
2108	Simplex::computeIntegerBounds(ArrayRef<DynamicAPInt> coeffs) {
2109	MaybeOptimum<DynamicAPInt> minRoundedUp(
2110	computeOptimum(direction: Simplex::Direction::Down, coeffs).map(f&: ceil));
2111	MaybeOptimum<DynamicAPInt> maxRoundedDown(
2112	computeOptimum(direction: Simplex::Direction::Up, coeffs).map(f&: floor));
2113	return {minRoundedUp, maxRoundedDown};
2114	}
2115
2116	bool Simplex::isFlatAlong(ArrayRef<DynamicAPInt> coeffs) {
2117	assert(!isEmpty() && "cannot check for flatness of empty simplex!");
2118	auto upOpt = computeOptimum(direction: Simplex::Direction::Up, coeffs);
2119	auto downOpt = computeOptimum(direction: Simplex::Direction::Down, coeffs);
2120
2121	if (!upOpt.isBounded())
2122	return false;
2123	if (!downOpt.isBounded())
2124	return false;
2125
2126	return *upOpt == *downOpt;
2127	}
2128
2129	void SimplexBase::print(raw_ostream &os) const {
2130	os << "rows = " << getNumRows() << ", columns = " << getNumColumns() << "\n";
2131	if (empty)
2132	os << "Simplex marked empty!\n";
2133	os << "var: ";
2134	for (unsigned i = 0; i < var.size(); ++i) {
2135	if (i > 0)
2136	os << ", ";
2137	var[i].print(os);
2138	}
2139	os << "\ncon: ";
2140	for (unsigned i = 0; i < con.size(); ++i) {
2141	if (i > 0)
2142	os << ", ";
2143	con[i].print(os);
2144	}
2145	os << '\n';
2146	for (unsigned row = 0, e = getNumRows(); row < e; ++row) {
2147	if (row > 0)
2148	os << ", ";
2149	os << "r" << row << ": " << rowUnknown[row];
2150	}
2151	os << '\n';
2152	os << "c0: denom, c1: const";
2153	for (unsigned col = 2, e = getNumColumns(); col < e; ++col)
2154	os << ", c" << col << ": " << colUnknown[col];
2155	os << '\n';
2156	PrintTableMetrics ptm = {.maxPreIndent: 0, .maxPostIndent: 0, .preAlign: "-"};
2157	for (unsigned row = 0, numRows = getNumRows(); row < numRows; ++row)
2158	for (unsigned col = 0, numCols = getNumColumns(); col < numCols; ++col)
2159	updatePrintMetrics<DynamicAPInt>(val: tableau(row, col), m&: ptm);
2160	unsigned MIN_SPACING = 1;
2161	for (unsigned row = 0, numRows = getNumRows(); row < numRows; ++row) {
2162	for (unsigned col = 0, numCols = getNumColumns(); col < numCols; ++col) {
2163	printWithPrintMetrics<DynamicAPInt>(os, val: tableau(row, col), minSpacing: MIN_SPACING,
2164	m: ptm);
2165	}
2166	os << '\n';
2167	}
2168	os << '\n';
2169	}
2170
2171	void SimplexBase::dump() const { print(os&: llvm::errs()); }
2172
2173	bool Simplex::isRationalSubsetOf(const IntegerRelation &rel) {
2174	if (isEmpty())
2175	return true;
2176
2177	for (unsigned i = 0, e = rel.getNumInequalities(); i < e; ++i)
2178	if (findIneqType(coeffs: rel.getInequality(idx: i)) != IneqType::Redundant)
2179	return false;
2180
2181	for (unsigned i = 0, e = rel.getNumEqualities(); i < e; ++i)
2182	if (!isRedundantEquality(coeffs: rel.getEquality(idx: i)))
2183	return false;
2184
2185	return true;
2186	}
2187
2188	/// Returns the type of the inequality with coefficients `coeffs`.
2189	/// Possible types are:
2190	/// Redundant The inequality is satisfied by all points in the polytope
2191	/// Cut The inequality is satisfied by some points, but not by others
2192	/// Separate The inequality is not satisfied by any point
2193	///
2194	/// Internally, this computes the minimum and the maximum the inequality with
2195	/// coefficients `coeffs` can take. If the minimum is >= 0, the inequality holds
2196	/// for all points in the polytope, so it is redundant. If the minimum is <= 0
2197	/// and the maximum is >= 0, the points in between the minimum and the
2198	/// inequality do not satisfy it, the points in between the inequality and the
2199	/// maximum satisfy it. Hence, it is a cut inequality. If both are < 0, no
2200	/// points of the polytope satisfy the inequality, which means it is a separate
2201	/// inequality.
2202	Simplex::IneqType Simplex::findIneqType(ArrayRef<DynamicAPInt> coeffs) {
2203	MaybeOptimum<Fraction> minimum = computeOptimum(direction: Direction::Down, coeffs);
2204	if (minimum.isBounded() && *minimum >= Fraction(0, 1)) {
2205	return IneqType::Redundant;
2206	}
2207	MaybeOptimum<Fraction> maximum = computeOptimum(direction: Direction::Up, coeffs);
2208	if ((!minimum.isBounded() || *minimum <= Fraction(0, 1)) &&
2209	(!maximum.isBounded() || *maximum >= Fraction(0, 1))) {
2210	return IneqType::Cut;
2211	}
2212	return IneqType::Separate;
2213	}
2214
2215	/// Checks whether the type of the inequality with coefficients `coeffs`
2216	/// is Redundant.
2217	bool Simplex::isRedundantInequality(ArrayRef<DynamicAPInt> coeffs) {
2218	assert(!empty &&
2219	"It is not meaningful to ask about redundancy in an empty set!");
2220	return findIneqType(coeffs) == IneqType::Redundant;
2221	}
2222
2223	/// Check whether the equality given by `coeffs == 0` is redundant given
2224	/// the existing constraints. This is redundant when `coeffs` is already
2225	/// always zero under the existing constraints. `coeffs` is always zero
2226	/// when the minimum and maximum value that `coeffs` can take are both zero.
2227	bool Simplex::isRedundantEquality(ArrayRef<DynamicAPInt> coeffs) {
2228	assert(!empty &&
2229	"It is not meaningful to ask about redundancy in an empty set!");
2230	MaybeOptimum<Fraction> minimum = computeOptimum(direction: Direction::Down, coeffs);
2231	MaybeOptimum<Fraction> maximum = computeOptimum(direction: Direction::Up, coeffs);
2232	assert((!minimum.isEmpty() && !maximum.isEmpty()) &&
2233	"Optima should be non-empty for a non-empty set");
2234	return minimum.isBounded() && maximum.isBounded() &&
2235	*maximum == Fraction(0, 1) && *minimum == Fraction(0, 1);
2236	}
2237