1	//===- IntegerRelation.cpp - MLIR IntegerRelation 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	// A class to represent an relation over integer tuples. A relation is
10	// represented as a constraint system over a space of tuples of integer valued
11	// variables supporting symbolic variables and existential quantification.
12	//
13	//===----------------------------------------------------------------------===//
14
15	#include "mlir/Analysis/Presburger/IntegerRelation.h"
16	#include "mlir/Analysis/Presburger/Fraction.h"
17	#include "mlir/Analysis/Presburger/LinearTransform.h"
18	#include "mlir/Analysis/Presburger/PWMAFunction.h"
19	#include "mlir/Analysis/Presburger/PresburgerRelation.h"
20	#include "mlir/Analysis/Presburger/PresburgerSpace.h"
21	#include "mlir/Analysis/Presburger/Simplex.h"
22	#include "mlir/Analysis/Presburger/Utils.h"
23	#include "llvm/ADT/DenseMap.h"
24	#include "llvm/ADT/STLExtras.h"
25	#include "llvm/ADT/Sequence.h"
26	#include "llvm/ADT/SmallBitVector.h"
27	#include "llvm/Support/Debug.h"
28	#include "llvm/Support/LogicalResult.h"
29	#include "llvm/Support/raw_ostream.h"
30	#include <algorithm>
31	#include <cassert>
32	#include <functional>
33	#include <memory>
34	#include <numeric>
35	#include <optional>
36	#include <sstream>
37	#include <string>
38	#include <utility>
39	#include <vector>
40
41	#define DEBUG_TYPE "presburger"
42
43	using namespace mlir;
44	using namespace presburger;
45
46	using llvm::SmallDenseMap;
47
48	std::unique_ptr<IntegerRelation> IntegerRelation::clone() const {
49	return std::make_unique<IntegerRelation>(args: *this);
50	}
51
52	std::unique_ptr<IntegerPolyhedron> IntegerPolyhedron::clone() const {
53	return std::make_unique<IntegerPolyhedron>(args: *this);
54	}
55
56	void IntegerRelation::setSpace(const PresburgerSpace &oSpace) {
57	assert(space.getNumVars() == oSpace.getNumVars() && "invalid space!");
58	space = oSpace;
59	}
60
61	void IntegerRelation::setSpaceExceptLocals(const PresburgerSpace &oSpace) {
62	assert(oSpace.getNumLocalVars() == 0 && "no locals should be present!");
63	assert(oSpace.getNumVars() <= getNumVars() && "invalid space!");
64	unsigned newNumLocals = getNumVars() - oSpace.getNumVars();
65	space = oSpace;
66	space.insertVar(kind: VarKind::Local, pos: 0, num: newNumLocals);
67	}
68
69	void IntegerRelation::setId(VarKind kind, unsigned i, Identifier id) {
70	assert(space.isUsingIds() &&
71	"space must be using identifiers to set an identifier");
72	assert(kind != VarKind::Local && "local variables cannot have identifiers");
73	assert(i < space.getNumVarKind(kind) && "invalid variable index");
74	space.setId(kind, pos: i, id);
75	}
76
77	ArrayRef<Identifier> IntegerRelation::getIds(VarKind kind) {
78	if (!space.isUsingIds())
79	space.resetIds();
80	return space.getIds(kind);
81	}
82
83	void IntegerRelation::append(const IntegerRelation &other) {
84	assert(space.isEqual(other.getSpace()) && "Spaces must be equal.");
85
86	inequalities.reserveRows(rows: inequalities.getNumRows() +
87	other.getNumInequalities());
88	equalities.reserveRows(rows: equalities.getNumRows() + other.getNumEqualities());
89
90	for (unsigned r = 0, e = other.getNumInequalities(); r < e; r++) {
91	addInequality(inEq: other.getInequality(idx: r));
92	}
93	for (unsigned r = 0, e = other.getNumEqualities(); r < e; r++) {
94	addEquality(eq: other.getEquality(idx: r));
95	}
96	}
97
98	IntegerRelation IntegerRelation::intersect(IntegerRelation other) const {
99	IntegerRelation result = *this;
100	result.mergeLocalVars(other);
101	result.append(other);
102	return result;
103	}
104
105	bool IntegerRelation::isEqual(const IntegerRelation &other) const {
106	assert(space.isCompatible(other.getSpace()) && "Spaces must be compatible.");
107	return PresburgerRelation(*this).isEqual(set: PresburgerRelation(other));
108	}
109
110	bool IntegerRelation::isObviouslyEqual(const IntegerRelation &other) const {
111	if (!space.isEqual(other: other.getSpace()))
112	return false;
113	if (getNumEqualities() != other.getNumEqualities())
114	return false;
115	if (getNumInequalities() != other.getNumInequalities())
116	return false;
117
118	unsigned cols = getNumCols();
119	for (unsigned i = 0, eqs = getNumEqualities(); i < eqs; ++i) {
120	for (unsigned j = 0; j < cols; ++j) {
121	if (atEq(i, j) != other.atEq(i, j))
122	return false;
123	}
124	}
125	for (unsigned i = 0, ineqs = getNumInequalities(); i < ineqs; ++i) {
126	for (unsigned j = 0; j < cols; ++j) {
127	if (atIneq(i, j) != other.atIneq(i, j))
128	return false;
129	}
130	}
131	return true;
132	}
133
134	bool IntegerRelation::isSubsetOf(const IntegerRelation &other) const {
135	assert(space.isCompatible(other.getSpace()) && "Spaces must be compatible.");
136	return PresburgerRelation(*this).isSubsetOf(set: PresburgerRelation(other));
137	}
138
139	MaybeOptimum<SmallVector<Fraction, 8>>
140	IntegerRelation::findRationalLexMin() const {
141	assert(getNumSymbolVars() == 0 && "Symbols are not supported!");
142	MaybeOptimum<SmallVector<Fraction, 8>> maybeLexMin =
143	LexSimplex(*this).findRationalLexMin();
144
145	if (!maybeLexMin.isBounded())
146	return maybeLexMin;
147
148	// The Simplex returns the lexmin over all the variables including locals. But
149	// locals are not actually part of the space and should not be returned in the
150	// result. Since the locals are placed last in the list of variables, they
151	// will be minimized last in the lexmin. So simply truncating out the locals
152	// from the end of the answer gives the desired lexmin over the dimensions.
153	assert(maybeLexMin->size() == getNumVars() &&
154	"Incorrect number of vars in lexMin!");
155	maybeLexMin->resize(N: getNumDimAndSymbolVars());
156	return maybeLexMin;
157	}
158
159	MaybeOptimum<SmallVector<DynamicAPInt, 8>>
160	IntegerRelation::findIntegerLexMin() const {
161	assert(getNumSymbolVars() == 0 && "Symbols are not supported!");
162	MaybeOptimum<SmallVector<DynamicAPInt, 8>> maybeLexMin =
163	LexSimplex(*this).findIntegerLexMin();
164
165	if (!maybeLexMin.isBounded())
166	return maybeLexMin.getKind();
167
168	// The Simplex returns the lexmin over all the variables including locals. But
169	// locals are not actually part of the space and should not be returned in the
170	// result. Since the locals are placed last in the list of variables, they
171	// will be minimized last in the lexmin. So simply truncating out the locals
172	// from the end of the answer gives the desired lexmin over the dimensions.
173	assert(maybeLexMin->size() == getNumVars() &&
174	"Incorrect number of vars in lexMin!");
175	maybeLexMin->resize(N: getNumDimAndSymbolVars());
176	return maybeLexMin;
177	}
178
179	static bool rangeIsZero(ArrayRef<DynamicAPInt> range) {
180	return llvm::all_of(Range&: range, P: [](const DynamicAPInt &x) { return x == 0; });
181	}
182
183	static void removeConstraintsInvolvingVarRange(IntegerRelation &poly,
184	unsigned begin, unsigned count) {
185	// We loop until i > 0 and index into i - 1 to avoid sign issues.
186	//
187	// We iterate backwards so that whether we remove constraint i - 1 or not, the
188	// next constraint to be tested is always i - 2.
189	for (unsigned i = poly.getNumEqualities(); i > 0; i--)
190	if (!rangeIsZero(range: poly.getEquality(idx: i - 1).slice(N: begin, M: count)))
191	poly.removeEquality(pos: i - 1);
192	for (unsigned i = poly.getNumInequalities(); i > 0; i--)
193	if (!rangeIsZero(range: poly.getInequality(idx: i - 1).slice(N: begin, M: count)))
194	poly.removeInequality(pos: i - 1);
195	}
196
197	IntegerRelation::CountsSnapshot IntegerRelation::getCounts() const {
198	return {getSpace(), getNumInequalities(), getNumEqualities()};
199	}
200
201	void IntegerRelation::truncateVarKind(VarKind kind, unsigned num) {
202	unsigned curNum = getNumVarKind(kind);
203	assert(num <= curNum && "Can't truncate to more vars!");
204	removeVarRange(kind, varStart: num, varLimit: curNum);
205	}
206
207	void IntegerRelation::truncateVarKind(VarKind kind,
208	const CountsSnapshot &counts) {
209	truncateVarKind(kind, num: counts.getSpace().getNumVarKind(kind));
210	}
211
212	void IntegerRelation::truncate(const CountsSnapshot &counts) {
213	truncateVarKind(kind: VarKind::Domain, counts);
214	truncateVarKind(kind: VarKind::Range, counts);
215	truncateVarKind(kind: VarKind::Symbol, counts);
216	truncateVarKind(kind: VarKind::Local, counts);
217	removeInequalityRange(start: counts.getNumIneqs(), end: getNumInequalities());
218	removeEqualityRange(start: counts.getNumEqs(), end: getNumEqualities());
219	}
220
221	PresburgerRelation IntegerRelation::computeReprWithOnlyDivLocals() const {
222	// If there are no locals, we're done.
223	if (getNumLocalVars() == 0)
224	return PresburgerRelation(*this);
225
226	// Move all the non-div locals to the end, as the current API to
227	// SymbolicLexOpt requires these to form a contiguous range.
228	//
229	// Take a copy so we can perform mutations.
230	IntegerRelation copy = *this;
231	std::vector<MaybeLocalRepr> reprs(getNumLocalVars());
232	copy.getLocalReprs(repr: &reprs);
233
234	// Iterate through all the locals. The last `numNonDivLocals` are the locals
235	// that have been scanned already and do not have division representations.
236	unsigned numNonDivLocals = 0;
237	unsigned offset = copy.getVarKindOffset(kind: VarKind::Local);
238	for (unsigned i = 0, e = copy.getNumLocalVars(); i < e - numNonDivLocals;) {
239	if (!reprs[i]) {
240	// Whenever we come across a local that does not have a division
241	// representation, we swap it to the `numNonDivLocals`-th last position
242	// and increment `numNonDivLocal`s. `reprs` also needs to be swapped.
243	copy.swapVar(posA: offset + i, posB: offset + e - numNonDivLocals - 1);
244	std::swap(a&: reprs[i], b&: reprs[e - numNonDivLocals - 1]);
245	++numNonDivLocals;
246	continue;
247	}
248	++i;
249	}
250
251	// If there are no non-div locals, we're done.
252	if (numNonDivLocals == 0)
253	return PresburgerRelation(*this);
254
255	// We computeSymbolicIntegerLexMin by considering the non-div locals as
256	// "non-symbols" and considering everything else as "symbols". This will
257	// compute a function mapping assignments to "symbols" to the
258	// lexicographically minimal valid assignment of "non-symbols", when a
259	// satisfying assignment exists. It separately returns the set of assignments
260	// to the "symbols" such that a satisfying assignment to the "non-symbols"
261	// exists but the lexmin is unbounded. We basically want to find the set of
262	// values of the "symbols" such that an assignment to the "non-symbols"
263	// exists, which is the union of the domain of the returned lexmin function
264	// and the returned set of assignments to the "symbols" that makes the lexmin
265	// unbounded.
266	SymbolicLexOpt lexminResult =
267	SymbolicLexSimplex(copy, /*symbolOffset*/ 0,
268	IntegerPolyhedron(PresburgerSpace::getSetSpace(
269	/*numDims=*/copy.getNumVars() - numNonDivLocals)))
270	.computeSymbolicIntegerLexMin();
271	PresburgerRelation result =
272	lexminResult.lexopt.getDomain().unionSet(set: lexminResult.unboundedDomain);
273
274	// The result set might lie in the wrong space -- all its ids are dims.
275	// Set it to the desired space and return.
276	PresburgerSpace space = getSpace();
277	space.removeVarRange(kind: VarKind::Local, varStart: 0, varLimit: getNumLocalVars());
278	result.setSpace(space);
279	return result;
280	}
281
282	SymbolicLexOpt IntegerRelation::findSymbolicIntegerLexMin() const {
283	// Symbol and Domain vars will be used as symbols for symbolic lexmin.
284	// In other words, for every value of the symbols and domain, return the
285	// lexmin value of the (range, locals).
286	llvm::SmallBitVector isSymbol(getNumVars(), false);
287	isSymbol.set(I: getVarKindOffset(kind: VarKind::Symbol),
288	E: getVarKindEnd(kind: VarKind::Symbol));
289	isSymbol.set(I: getVarKindOffset(kind: VarKind::Domain),
290	E: getVarKindEnd(kind: VarKind::Domain));
291	// Compute the symbolic lexmin of the dims and locals, with the symbols being
292	// the actual symbols of this set.
293	// The resultant space of lexmin is the space of the relation itself.
294	SymbolicLexOpt result =
295	SymbolicLexSimplex(*this,
296	IntegerPolyhedron(PresburgerSpace::getSetSpace(
297	/*numDims=*/getNumDomainVars(),
298	/*numSymbols=*/getNumSymbolVars())),
299	isSymbol)
300	.computeSymbolicIntegerLexMin();
301
302	// We want to return only the lexmin over the dims, so strip the locals from
303	// the computed lexmin.
304	result.lexopt.removeOutputs(start: result.lexopt.getNumOutputs() - getNumLocalVars(),
305	end: result.lexopt.getNumOutputs());
306	return result;
307	}
308
309	/// findSymbolicIntegerLexMax is implemented using findSymbolicIntegerLexMin as
310	/// follows:
311	/// 1. A new relation is created which is `this` relation with the sign of
312	/// each dimension variable in range flipped;
313	/// 2. findSymbolicIntegerLexMin is called on the range negated relation to
314	/// compute the negated lexmax of `this` relation;
315	/// 3. The sign of the negated lexmax is flipped and returned.
316	SymbolicLexOpt IntegerRelation::findSymbolicIntegerLexMax() const {
317	IntegerRelation flippedRel = *this;
318	// Flip range sign by flipping the sign of range variables in all constraints.
319	for (unsigned j = getNumDomainVars(),
320	b = getNumDomainVars() + getNumRangeVars();
321	j < b; j++) {
322	for (unsigned i = 0, a = getNumEqualities(); i < a; i++)
323	flippedRel.atEq(i, j) = -1 * atEq(i, j);
324	for (unsigned i = 0, a = getNumInequalities(); i < a; i++)
325	flippedRel.atIneq(i, j) = -1 * atIneq(i, j);
326	}
327	// Compute negated lexmax by computing lexmin.
328	SymbolicLexOpt flippedSymbolicIntegerLexMax =
329	flippedRel.findSymbolicIntegerLexMin(),
330	symbolicIntegerLexMax(
331	flippedSymbolicIntegerLexMax.lexopt.getSpace());
332	// Get lexmax by flipping range sign in the PWMA constraints.
333	for (auto &flippedPiece :
334	flippedSymbolicIntegerLexMax.lexopt.getAllPieces()) {
335	IntMatrix mat = flippedPiece.output.getOutputMatrix();
336	for (unsigned i = 0, e = mat.getNumRows(); i < e; i++)
337	mat.negateRow(row: i);
338	MultiAffineFunction maf(flippedPiece.output.getSpace(), mat);
339	PWMAFunction::Piece piece = {.domain: flippedPiece.domain, .output: maf};
340	symbolicIntegerLexMax.lexopt.addPiece(piece);
341	}
342	symbolicIntegerLexMax.unboundedDomain =
343	flippedSymbolicIntegerLexMax.unboundedDomain;
344	return symbolicIntegerLexMax;
345	}
346
347	PresburgerRelation
348	IntegerRelation::subtract(const PresburgerRelation &set) const {
349	return PresburgerRelation(*this).subtract(set);
350	}
351
352	unsigned IntegerRelation::insertVar(VarKind kind, unsigned pos, unsigned num) {
353	assert(pos <= getNumVarKind(kind));
354
355	unsigned insertPos = space.insertVar(kind, pos, num);
356	inequalities.insertColumns(pos: insertPos, count: num);
357	equalities.insertColumns(pos: insertPos, count: num);
358	return insertPos;
359	}
360
361	unsigned IntegerRelation::appendVar(VarKind kind, unsigned num) {
362	unsigned pos = getNumVarKind(kind);
363	return insertVar(kind, pos, num);
364	}
365
366	void IntegerRelation::addEquality(ArrayRef<DynamicAPInt> eq) {
367	assert(eq.size() == getNumCols());
368	unsigned row = equalities.appendExtraRow();
369	for (unsigned i = 0, e = eq.size(); i < e; ++i)
370	equalities(row, i) = eq[i];
371	}
372
373	void IntegerRelation::addInequality(ArrayRef<DynamicAPInt> inEq) {
374	assert(inEq.size() == getNumCols());
375	unsigned row = inequalities.appendExtraRow();
376	for (unsigned i = 0, e = inEq.size(); i < e; ++i)
377	inequalities(row, i) = inEq[i];
378	}
379
380	void IntegerRelation::removeVar(VarKind kind, unsigned pos) {
381	removeVarRange(kind, varStart: pos, varLimit: pos + 1);
382	}
383
384	void IntegerRelation::removeVar(unsigned pos) { removeVarRange(varStart: pos, varLimit: pos + 1); }
385
386	void IntegerRelation::removeVarRange(VarKind kind, unsigned varStart,
387	unsigned varLimit) {
388	assert(varLimit <= getNumVarKind(kind));
389
390	if (varStart >= varLimit)
391	return;
392
393	// Remove eliminated variables from the constraints.
394	unsigned offset = getVarKindOffset(kind);
395	equalities.removeColumns(pos: offset + varStart, count: varLimit - varStart);
396	inequalities.removeColumns(pos: offset + varStart, count: varLimit - varStart);
397
398	// Remove eliminated variables from the space.
399	space.removeVarRange(kind, varStart, varLimit);
400	}
401
402	void IntegerRelation::removeVarRange(unsigned varStart, unsigned varLimit) {
403	assert(varLimit <= getNumVars());
404
405	if (varStart >= varLimit)
406	return;
407
408	// Helper function to remove vars of the specified kind in the given range
409	// [start, limit), The range is absolute (i.e. it is not relative to the kind
410	// of variable). Also updates `limit` to reflect the deleted variables.
411	auto removeVarKindInRange = [this](VarKind kind, unsigned &start,
412	unsigned &limit) {
413	if (start >= limit)
414	return;
415
416	unsigned offset = getVarKindOffset(kind);
417	unsigned num = getNumVarKind(kind);
418
419	// Get `start`, `limit` relative to the specified kind.
420	unsigned relativeStart =
421	start <= offset ? 0 : std::min(a: num, b: start - offset);
422	unsigned relativeLimit =
423	limit <= offset ? 0 : std::min(a: num, b: limit - offset);
424
425	// Remove vars of the specified kind in the relative range.
426	removeVarRange(kind, varStart: relativeStart, varLimit: relativeLimit);
427
428	// Update `limit` to reflect deleted variables.
429	// `start` does not need to be updated because any variables that are
430	// deleted are after position `start`.
431	limit -= relativeLimit - relativeStart;
432	};
433
434	removeVarKindInRange(VarKind::Domain, varStart, varLimit);
435	removeVarKindInRange(VarKind::Range, varStart, varLimit);
436	removeVarKindInRange(VarKind::Symbol, varStart, varLimit);
437	removeVarKindInRange(VarKind::Local, varStart, varLimit);
438	}
439
440	void IntegerRelation::removeEquality(unsigned pos) {
441	equalities.removeRow(pos);
442	}
443
444	void IntegerRelation::removeInequality(unsigned pos) {
445	inequalities.removeRow(pos);
446	}
447
448	void IntegerRelation::removeEqualityRange(unsigned start, unsigned end) {
449	if (start >= end)
450	return;
451	equalities.removeRows(pos: start, count: end - start);
452	}
453
454	void IntegerRelation::removeInequalityRange(unsigned start, unsigned end) {
455	if (start >= end)
456	return;
457	inequalities.removeRows(pos: start, count: end - start);
458	}
459
460	void IntegerRelation::swapVar(unsigned posA, unsigned posB) {
461	assert(posA < getNumVars() && "invalid position A");
462	assert(posB < getNumVars() && "invalid position B");
463
464	if (posA == posB)
465	return;
466
467	VarKind kindA = space.getVarKindAt(pos: posA);
468	VarKind kindB = space.getVarKindAt(pos: posB);
469	unsigned relativePosA = posA - getVarKindOffset(kind: kindA);
470	unsigned relativePosB = posB - getVarKindOffset(kind: kindB);
471	space.swapVar(kindA, kindB, posA: relativePosA, posB: relativePosB);
472
473	inequalities.swapColumns(column: posA, otherColumn: posB);
474	equalities.swapColumns(column: posA, otherColumn: posB);
475	}
476
477	void IntegerRelation::clearConstraints() {
478	equalities.resizeVertically(newNRows: 0);
479	inequalities.resizeVertically(newNRows: 0);
480	}
481
482	/// Gather all lower and upper bounds of the variable at `pos`, and
483	/// optionally any equalities on it. In addition, the bounds are to be
484	/// independent of variables in position range [`offset`, `offset` + `num`).
485	void IntegerRelation::getLowerAndUpperBoundIndices(
486	unsigned pos, SmallVectorImpl<unsigned> *lbIndices,
487	SmallVectorImpl<unsigned> *ubIndices, SmallVectorImpl<unsigned> *eqIndices,
488	unsigned offset, unsigned num) const {
489	assert(pos < getNumVars() && "invalid position");
490	assert(offset + num < getNumCols() && "invalid range");
491
492	// Checks for a constraint that has a non-zero coeff for the variables in
493	// the position range [offset, offset + num) while ignoring `pos`.
494	auto containsConstraintDependentOnRange = [&](unsigned r, bool isEq) {
495	unsigned c, f;
496	auto cst = isEq ? getEquality(idx: r) : getInequality(idx: r);
497	for (c = offset, f = offset + num; c < f; ++c) {
498	if (c == pos)
499	continue;
500	if (cst[c] != 0)
501	break;
502	}
503	return c < f;
504	};
505
506	// Gather all lower bounds and upper bounds of the variable. Since the
507	// canonical form c_1*x_1 + c_2*x_2 + ... + c_0 >= 0, a constraint is a lower
508	// bound for x_i if c_i >= 1, and an upper bound if c_i <= -1.
509	for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
510	// The bounds are to be independent of [offset, offset + num) columns.
511	if (containsConstraintDependentOnRange(r, /*isEq=*/false))
512	continue;
513	if (atIneq(i: r, j: pos) >= 1) {
514	// Lower bound.
515	lbIndices->emplace_back(Args&: r);
516	} else if (atIneq(i: r, j: pos) <= -1) {
517	// Upper bound.
518	ubIndices->emplace_back(Args&: r);
519	}
520	}
521
522	// An equality is both a lower and upper bound. Record any equalities
523	// involving the pos^th variable.
524	if (!eqIndices)
525	return;
526
527	for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
528	if (atEq(i: r, j: pos) == 0)
529	continue;
530	if (containsConstraintDependentOnRange(r, /*isEq=*/true))
531	continue;
532	eqIndices->emplace_back(Args&: r);
533	}
534	}
535
536	bool IntegerRelation::hasConsistentState() const {
537	if (!inequalities.hasConsistentState())
538	return false;
539	if (!equalities.hasConsistentState())
540	return false;
541	return true;
542	}
543
544	void IntegerRelation::setAndEliminate(unsigned pos,
545	ArrayRef<DynamicAPInt> values) {
546	if (values.empty())
547	return;
548	assert(pos + values.size() <= getNumVars() &&
549	"invalid position or too many values");
550	// Setting x_j = p in sum_i a_i x_i + c is equivalent to adding p*a_j to the
551	// constant term and removing the var x_j. We do this for all the vars
552	// pos, pos + 1, ... pos + values.size() - 1.
553	unsigned constantColPos = getNumCols() - 1;
554	for (unsigned i = 0, numVals = values.size(); i < numVals; ++i)
555	inequalities.addToColumn(sourceColumn: i + pos, targetColumn: constantColPos, scale: values[i]);
556	for (unsigned i = 0, numVals = values.size(); i < numVals; ++i)
557	equalities.addToColumn(sourceColumn: i + pos, targetColumn: constantColPos, scale: values[i]);
558	removeVarRange(varStart: pos, varLimit: pos + values.size());
559	}
560
561	void IntegerRelation::clearAndCopyFrom(const IntegerRelation &other) {
562	*this = other;
563	}
564
565	std::optional<unsigned>
566	IntegerRelation::findConstraintWithNonZeroAt(unsigned colIdx, bool isEq) const {
567	assert(colIdx < getNumCols() && "position out of bounds");
568	auto at = [&](unsigned rowIdx) -> DynamicAPInt {
569	return isEq ? atEq(i: rowIdx, j: colIdx) : atIneq(i: rowIdx, j: colIdx);
570	};
571	unsigned e = isEq ? getNumEqualities() : getNumInequalities();
572	for (unsigned rowIdx = 0; rowIdx < e; ++rowIdx) {
573	if (at(rowIdx) != 0)
574	return rowIdx;
575	}
576	return std::nullopt;
577	}
578
579	void IntegerRelation::normalizeConstraintsByGCD() {
580	for (unsigned i = 0, e = getNumEqualities(); i < e; ++i)
581	equalities.normalizeRow(row: i);
582	for (unsigned i = 0, e = getNumInequalities(); i < e; ++i)
583	inequalities.normalizeRow(row: i);
584	}
585
586	bool IntegerRelation::hasInvalidConstraint() const {
587	assert(hasConsistentState());
588	auto check = [&](bool isEq) -> bool {
589	unsigned numCols = getNumCols();
590	unsigned numRows = isEq ? getNumEqualities() : getNumInequalities();
591	for (unsigned i = 0, e = numRows; i < e; ++i) {
592	unsigned j;
593	for (j = 0; j < numCols - 1; ++j) {
594	DynamicAPInt v = isEq ? atEq(i, j) : atIneq(i, j);
595	// Skip rows with non-zero variable coefficients.
596	if (v != 0)
597	break;
598	}
599	if (j < numCols - 1) {
600	continue;
601	}
602	// Check validity of constant term at 'numCols - 1' w.r.t 'isEq'.
603	// Example invalid constraints include: '1 == 0' or '-1 >= 0'
604	DynamicAPInt v = isEq ? atEq(i, j: numCols - 1) : atIneq(i, j: numCols - 1);
605	if ((isEq && v != 0) || (!isEq && v < 0)) {
606	return true;
607	}
608	}
609	return false;
610	};
611	if (check(/*isEq=*/true))
612	return true;
613	return check(/*isEq=*/false);
614	}
615
616	/// Eliminate variable from constraint at `rowIdx` based on coefficient at
617	/// pivotRow, pivotCol. Columns in range [elimColStart, pivotCol) will not be
618	/// updated as they have already been eliminated.
619	static void eliminateFromConstraint(IntegerRelation *constraints,
620	unsigned rowIdx, unsigned pivotRow,
621	unsigned pivotCol, unsigned elimColStart,
622	bool isEq) {
623	// Skip if equality 'rowIdx' if same as 'pivotRow'.
624	if (isEq && rowIdx == pivotRow)
625	return;
626	auto at = [&](unsigned i, unsigned j) -> DynamicAPInt {
627	return isEq ? constraints->atEq(i, j) : constraints->atIneq(i, j);
628	};
629	DynamicAPInt leadCoeff = at(rowIdx, pivotCol);
630	// Skip if leading coefficient at 'rowIdx' is already zero.
631	if (leadCoeff == 0)
632	return;
633	DynamicAPInt pivotCoeff = constraints->atEq(i: pivotRow, j: pivotCol);
634	int sign = (leadCoeff * pivotCoeff > 0) ? -1 : 1;
635	DynamicAPInt lcm = llvm::lcm(A: pivotCoeff, B: leadCoeff);
636	DynamicAPInt pivotMultiplier = sign * (lcm / abs(X: pivotCoeff));
637	DynamicAPInt rowMultiplier = lcm / abs(X: leadCoeff);
638
639	unsigned numCols = constraints->getNumCols();
640	for (unsigned j = 0; j < numCols; ++j) {
641	// Skip updating column 'j' if it was just eliminated.
642	if (j >= elimColStart && j < pivotCol)
643	continue;
644	DynamicAPInt v = pivotMultiplier * constraints->atEq(i: pivotRow, j) +
645	rowMultiplier * at(rowIdx, j);
646	isEq ? constraints->atEq(i: rowIdx, j) = v
647	: constraints->atIneq(i: rowIdx, j) = v;
648	}
649	}
650
651	/// Returns the position of the variable that has the minimum <number of lower
652	/// bounds> times <number of upper bounds> from the specified range of
653	/// variables [start, end). It is often best to eliminate in the increasing
654	/// order of these counts when doing Fourier-Motzkin elimination since FM adds
655	/// that many new constraints.
656	static unsigned getBestVarToEliminate(const IntegerRelation &cst,
657	unsigned start, unsigned end) {
658	assert(start < cst.getNumVars() && end < cst.getNumVars() + 1);
659
660	auto getProductOfNumLowerUpperBounds = [&](unsigned pos) {
661	unsigned numLb = 0;
662	unsigned numUb = 0;
663	for (unsigned r = 0, e = cst.getNumInequalities(); r < e; r++) {
664	if (cst.atIneq(i: r, j: pos) > 0) {
665	++numLb;
666	} else if (cst.atIneq(i: r, j: pos) < 0) {
667	++numUb;
668	}
669	}
670	return numLb * numUb;
671	};
672
673	unsigned minLoc = start;
674	unsigned min = getProductOfNumLowerUpperBounds(start);
675	for (unsigned c = start + 1; c < end; c++) {
676	unsigned numLbUbProduct = getProductOfNumLowerUpperBounds(c);
677	if (numLbUbProduct < min) {
678	min = numLbUbProduct;
679	minLoc = c;
680	}
681	}
682	return minLoc;
683	}
684
685	// Checks for emptiness of the set by eliminating variables successively and
686	// using the GCD test (on all equality constraints) and checking for trivially
687	// invalid constraints. Returns 'true' if the constraint system is found to be
688	// empty; false otherwise.
689	bool IntegerRelation::isEmpty() const {
690	if (isEmptyByGCDTest() || hasInvalidConstraint())
691	return true;
692
693	IntegerRelation tmpCst(*this);
694
695	// First, eliminate as many local variables as possible using equalities.
696	tmpCst.removeRedundantLocalVars();
697	if (tmpCst.isEmptyByGCDTest() || tmpCst.hasInvalidConstraint())
698	return true;
699
700	// Eliminate as many variables as possible using Gaussian elimination.
701	unsigned currentPos = 0;
702	while (currentPos < tmpCst.getNumVars()) {
703	tmpCst.gaussianEliminateVars(posStart: currentPos, posLimit: tmpCst.getNumVars());
704	++currentPos;
705	// We check emptiness through trivial checks after eliminating each ID to
706	// detect emptiness early. Since the checks isEmptyByGCDTest() and
707	// hasInvalidConstraint() are linear time and single sweep on the constraint
708	// buffer, this appears reasonable - but can optimize in the future.
709	if (tmpCst.hasInvalidConstraint() || tmpCst.isEmptyByGCDTest())
710	return true;
711	}
712
713	// Eliminate the remaining using FM.
714	for (unsigned i = 0, e = tmpCst.getNumVars(); i < e; i++) {
715	tmpCst.fourierMotzkinEliminate(
716	pos: getBestVarToEliminate(cst: tmpCst, start: 0, end: tmpCst.getNumVars()));
717	// Check for a constraint explosion. This rarely happens in practice, but
718	// this check exists as a safeguard against improperly constructed
719	// constraint systems or artificially created arbitrarily complex systems
720	// that aren't the intended use case for IntegerRelation. This is
721	// needed since FM has a worst case exponential complexity in theory.
722	if (tmpCst.getNumConstraints() >= kExplosionFactor * getNumVars()) {
723	LLVM_DEBUG(llvm::dbgs() << "FM constraint explosion detected\n");
724	return false;
725	}
726
727	// FM wouldn't have modified the equalities in any way. So no need to again
728	// run GCD test. Check for trivial invalid constraints.
729	if (tmpCst.hasInvalidConstraint())
730	return true;
731	}
732	return false;
733	}
734
735	bool IntegerRelation::isObviouslyEmpty() const {
736	return isEmptyByGCDTest() || hasInvalidConstraint();
737	}
738
739	// Runs the GCD test on all equality constraints. Returns 'true' if this test
740	// fails on any equality. Returns 'false' otherwise.
741	// This test can be used to disprove the existence of a solution. If it returns
742	// true, no integer solution to the equality constraints can exist.
743	//
744	// GCD test definition:
745	//
746	// The equality constraint:
747	//
748	// c_1*x_1 + c_2*x_2 + ... + c_n*x_n = c_0
749	//
750	// has an integer solution iff:
751	//
752	// GCD of c_1, c_2, ..., c_n divides c_0.
753	bool IntegerRelation::isEmptyByGCDTest() const {
754	assert(hasConsistentState());
755	unsigned numCols = getNumCols();
756	for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
757	DynamicAPInt gcd = abs(X: atEq(i, j: 0));
758	for (unsigned j = 1; j < numCols - 1; ++j) {
759	gcd = llvm::gcd(A: gcd, B: abs(X: atEq(i, j)));
760	}
761	DynamicAPInt v = abs(X: atEq(i, j: numCols - 1));
762	if (gcd > 0 && (v % gcd != 0)) {
763	return true;
764	}
765	}
766	return false;
767	}
768
769	// Returns a matrix where each row is a vector along which the polytope is
770	// bounded. The span of the returned vectors is guaranteed to contain all
771	// such vectors. The returned vectors are NOT guaranteed to be linearly
772	// independent. This function should not be called on empty sets.
773	//
774	// It is sufficient to check the perpendiculars of the constraints, as the set
775	// of perpendiculars which are bounded must span all bounded directions.
776	IntMatrix IntegerRelation::getBoundedDirections() const {
777	// Note that it is necessary to add the equalities too (which the constructor
778	// does) even though we don't need to check if they are bounded; whether an
779	// inequality is bounded or not depends on what other constraints, including
780	// equalities, are present.
781	Simplex simplex(*this);
782
783	assert(!simplex.isEmpty() && "It is not meaningful to ask whether a "
784	"direction is bounded in an empty set.");
785
786	SmallVector<unsigned, 8> boundedIneqs;
787	// The constructor adds the inequalities to the simplex first, so this
788	// processes all the inequalities.
789	for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
790	if (simplex.isBoundedAlongConstraint(constraintIndex: i))
791	boundedIneqs.emplace_back(Args&: i);
792	}
793
794	// The direction vector is given by the coefficients and does not include the
795	// constant term, so the matrix has one fewer column.
796	unsigned dirsNumCols = getNumCols() - 1;
797	IntMatrix dirs(boundedIneqs.size() + getNumEqualities(), dirsNumCols);
798
799	// Copy the bounded inequalities.
800	unsigned row = 0;
801	for (unsigned i : boundedIneqs) {
802	for (unsigned col = 0; col < dirsNumCols; ++col)
803	dirs(row, col) = atIneq(i, j: col);
804	++row;
805	}
806
807	// Copy the equalities. All the equalities' perpendiculars are bounded.
808	for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
809	for (unsigned col = 0; col < dirsNumCols; ++col)
810	dirs(row, col) = atEq(i, j: col);
811	++row;
812	}
813
814	return dirs;
815	}
816
817	bool IntegerRelation::isIntegerEmpty() const { return !findIntegerSample(); }
818
819	/// Let this set be S. If S is bounded then we directly call into the GBR
820	/// sampling algorithm. Otherwise, there are some unbounded directions, i.e.,
821	/// vectors v such that S extends to infinity along v or -v. In this case we
822	/// use an algorithm described in the integer set library (isl) manual and used
823	/// by the isl_set_sample function in that library. The algorithm is:
824	///
825	/// 1) Apply a unimodular transform T to S to obtain S*T, such that all
826	/// dimensions in which S*T is bounded lie in the linear span of a prefix of the
827	/// dimensions.
828	///
829	/// 2) Construct a set B by removing all constraints that involve
830	/// the unbounded dimensions and then deleting the unbounded dimensions. Note
831	/// that B is a Bounded set.
832	///
833	/// 3) Try to obtain a sample from B using the GBR sampling
834	/// algorithm. If no sample is found, return that S is empty.
835	///
836	/// 4) Otherwise, substitute the obtained sample into S*T to obtain a set
837	/// C. C is a full-dimensional Cone and always contains a sample.
838	///
839	/// 5) Obtain an integer sample from C.
840	///
841	/// 6) Return T*v, where v is the concatenation of the samples from B and C.
842	///
843	/// The following is a sketch of a proof that
844	/// a) If the algorithm returns empty, then S is empty.
845	/// b) If the algorithm returns a sample, it is a valid sample in S.
846	///
847	/// The algorithm returns empty only if B is empty, in which case S*T is
848	/// certainly empty since B was obtained by removing constraints and then
849	/// deleting unconstrained dimensions from S*T. Since T is unimodular, a vector
850	/// v is in S*T iff T*v is in S. So in this case, since
851	/// S*T is empty, S is empty too.
852	///
853	/// Otherwise, the algorithm substitutes the sample from B into S*T. All the
854	/// constraints of S*T that did not involve unbounded dimensions are satisfied
855	/// by this substitution. All dimensions in the linear span of the dimensions
856	/// outside the prefix are unbounded in S*T (step 1). Substituting values for
857	/// the bounded dimensions cannot make these dimensions bounded, and these are
858	/// the only remaining dimensions in C, so C is unbounded along every vector (in
859	/// the positive or negative direction, or both). C is hence a full-dimensional
860	/// cone and therefore always contains an integer point.
861	///
862	/// Concatenating the samples from B and C gives a sample v in S*T, so the
863	/// returned sample T*v is a sample in S.
864	std::optional<SmallVector<DynamicAPInt, 8>>
865	IntegerRelation::findIntegerSample() const {
866	// First, try the GCD test heuristic.
867	if (isEmptyByGCDTest())
868	return {};
869
870	Simplex simplex(*this);
871	if (simplex.isEmpty())
872	return {};
873
874	// For a bounded set, we directly call into the GBR sampling algorithm.
875	if (!simplex.isUnbounded())
876	return simplex.findIntegerSample();
877
878	// The set is unbounded. We cannot directly use the GBR algorithm.
879	//
880	// m is a matrix containing, in each row, a vector in which S is
881	// bounded, such that the linear span of all these dimensions contains all
882	// bounded dimensions in S.
883	IntMatrix m = getBoundedDirections();
884	// In column echelon form, each row of m occupies only the first rank(m)
885	// columns and has zeros on the other columns. The transform T that brings S
886	// to column echelon form is unimodular as well, so this is a suitable
887	// transform to use in step 1 of the algorithm.
888	std::pair<unsigned, LinearTransform> result =
889	LinearTransform::makeTransformToColumnEchelon(m);
890	const LinearTransform &transform = result.second;
891	// 1) Apply T to S to obtain S*T.
892	IntegerRelation transformedSet = transform.applyTo(rel: *this);
893
894	// 2) Remove the unbounded dimensions and constraints involving them to
895	// obtain a bounded set.
896	IntegerRelation boundedSet(transformedSet);
897	unsigned numBoundedDims = result.first;
898	unsigned numUnboundedDims = getNumVars() - numBoundedDims;
899	removeConstraintsInvolvingVarRange(poly&: boundedSet, begin: numBoundedDims,
900	count: numUnboundedDims);
901	boundedSet.removeVarRange(varStart: numBoundedDims, varLimit: boundedSet.getNumVars());
902
903	// 3) Try to obtain a sample from the bounded set.
904	std::optional<SmallVector<DynamicAPInt, 8>> boundedSample =
905	Simplex(boundedSet).findIntegerSample();
906	if (!boundedSample)
907	return {};
908	assert(boundedSet.containsPoint(*boundedSample) &&
909	"Simplex returned an invalid sample!");
910
911	// 4) Substitute the values of the bounded dimensions into S*T to obtain a
912	// full-dimensional cone, which necessarily contains an integer sample.
913	transformedSet.setAndEliminate(pos: 0, values: *boundedSample);
914	IntegerRelation &cone = transformedSet;
915
916	// 5) Obtain an integer sample from the cone.
917	//
918	// We shrink the cone such that for any rational point in the shrunken cone,
919	// rounding up each of the point's coordinates produces a point that still
920	// lies in the original cone.
921	//
922	// Rounding up a point x adds a number e_i in [0, 1) to each coordinate x_i.
923	// For each inequality sum_i a_i x_i + c >= 0 in the original cone, the
924	// shrunken cone will have the inequality tightened by some amount s, such
925	// that if x satisfies the shrunken cone's tightened inequality, then x + e
926	// satisfies the original inequality, i.e.,
927	//
928	// sum_i a_i x_i + c + s >= 0 implies sum_i a_i (x_i + e_i) + c >= 0
929	//
930	// for any e_i values in [0, 1). In fact, we will handle the slightly more
931	// general case where e_i can be in [0, 1]. For example, consider the
932	// inequality 2x_1 - 3x_2 - 7x_3 - 6 >= 0, and let x = (3, 0, 0). How low
933	// could the LHS go if we added a number in [0, 1] to each coordinate? The LHS
934	// is minimized when we add 1 to the x_i with negative coefficient a_i and
935	// keep the other x_i the same. In the example, we would get x = (3, 1, 1),
936	// changing the value of the LHS by -3 + -7 = -10.
937	//
938	// In general, the value of the LHS can change by at most the sum of the
939	// negative a_i, so we accomodate this by shifting the inequality by this
940	// amount for the shrunken cone.
941	for (unsigned i = 0, e = cone.getNumInequalities(); i < e; ++i) {
942	for (unsigned j = 0; j < cone.getNumVars(); ++j) {
943	DynamicAPInt coeff = cone.atIneq(i, j);
944	if (coeff < 0)
945	cone.atIneq(i, j: cone.getNumVars()) += coeff;
946	}
947	}
948
949	// Obtain an integer sample in the cone by rounding up a rational point from
950	// the shrunken cone. Shrinking the cone amounts to shifting its apex
951	// "inwards" without changing its "shape"; the shrunken cone is still a
952	// full-dimensional cone and is hence non-empty.
953	Simplex shrunkenConeSimplex(cone);
954	assert(!shrunkenConeSimplex.isEmpty() && "Shrunken cone cannot be empty!");
955
956	// The sample will always exist since the shrunken cone is non-empty.
957	SmallVector<Fraction, 8> shrunkenConeSample =
958	*shrunkenConeSimplex.getRationalSample();
959
960	SmallVector<DynamicAPInt, 8> coneSample(
961	llvm::map_range(C&: shrunkenConeSample, F: ceil));
962
963	// 6) Return transform * concat(boundedSample, coneSample).
964	SmallVector<DynamicAPInt, 8> &sample = *boundedSample;
965	sample.append(in_start: coneSample.begin(), in_end: coneSample.end());
966	return transform.postMultiplyWithColumn(colVec: sample);
967	}
968
969	/// Helper to evaluate an affine expression at a point.
970	/// The expression is a list of coefficients for the dimensions followed by the
971	/// constant term.
972	static DynamicAPInt valueAt(ArrayRef<DynamicAPInt> expr,
973	ArrayRef<DynamicAPInt> point) {
974	assert(expr.size() == 1 + point.size() &&
975	"Dimensionalities of point and expression don't match!");
976	DynamicAPInt value = expr.back();
977	for (unsigned i = 0; i < point.size(); ++i)
978	value += expr[i] * point[i];
979	return value;
980	}
981
982	/// A point satisfies an equality iff the value of the equality at the
983	/// expression is zero, and it satisfies an inequality iff the value of the
984	/// inequality at that point is non-negative.
985	bool IntegerRelation::containsPoint(ArrayRef<DynamicAPInt> point) const {
986	for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
987	if (valueAt(expr: getEquality(idx: i), point) != 0)
988	return false;
989	}
990	for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
991	if (valueAt(expr: getInequality(idx: i), point) < 0)
992	return false;
993	}
994	return true;
995	}
996
997	/// Just substitute the values given and check if an integer sample exists for
998	/// the local vars.
999	///
1000	/// TODO: this could be made more efficient by handling divisions separately.
1001	/// Instead of finding an integer sample over all the locals, we can first
1002	/// compute the values of the locals that have division representations and
1003	/// only use the integer emptiness check for the locals that don't have this.
1004	/// Handling this correctly requires ordering the divs, though.
1005	std::optional<SmallVector<DynamicAPInt, 8>>
1006	IntegerRelation::containsPointNoLocal(ArrayRef<DynamicAPInt> point) const {
1007	assert(point.size() == getNumVars() - getNumLocalVars() &&
1008	"Point should contain all vars except locals!");
1009	assert(getVarKindOffset(VarKind::Local) == getNumVars() - getNumLocalVars() &&
1010	"This function depends on locals being stored last!");
1011	IntegerRelation copy = *this;
1012	copy.setAndEliminate(pos: 0, values: point);
1013	return copy.findIntegerSample();
1014	}
1015
1016	DivisionRepr
1017	IntegerRelation::getLocalReprs(std::vector<MaybeLocalRepr> *repr) const {
1018	SmallVector<bool, 8> foundRepr(getNumVars(), false);
1019	for (unsigned i = 0, e = getNumDimAndSymbolVars(); i < e; ++i)
1020	foundRepr[i] = true;
1021
1022	unsigned localOffset = getVarKindOffset(kind: VarKind::Local);
1023	DivisionRepr divs(getNumVars(), getNumLocalVars());
1024	bool changed;
1025	do {
1026	// Each time changed is true, at end of this iteration, one or more local
1027	// vars have been detected as floor divs.
1028	changed = false;
1029	for (unsigned i = 0, e = getNumLocalVars(); i < e; ++i) {
1030	if (!foundRepr[i + localOffset]) {
1031	MaybeLocalRepr res =
1032	computeSingleVarRepr(cst: *this, foundRepr, pos: localOffset + i,
1033	dividend: divs.getDividend(i), divisor&: divs.getDenom(i));
1034	if (!res) {
1035	// No representation was found, so clear the representation and
1036	// continue.
1037	divs.clearRepr(i);
1038	continue;
1039	}
1040	foundRepr[localOffset + i] = true;
1041	if (repr)
1042	(*repr)[i] = res;
1043	changed = true;
1044	}
1045	}
1046	} while (changed);
1047
1048	return divs;
1049	}
1050
1051	/// Tightens inequalities given that we are dealing with integer spaces. This is
1052	/// analogous to the GCD test but applied to inequalities. The constant term can
1053	/// be reduced to the preceding multiple of the GCD of the coefficients, i.e.,
1054	/// 64*i - 100 >= 0 => 64*i - 128 >= 0 (since 'i' is an integer). This is a
1055	/// fast method - linear in the number of coefficients.
1056	// Example on how this affects practical cases: consider the scenario:
1057	// 64*i >= 100, j = 64*i; without a tightening, elimination of i would yield
1058	// j >= 100 instead of the tighter (exact) j >= 128.
1059	void IntegerRelation::gcdTightenInequalities() {
1060	unsigned numCols = getNumCols();
1061	for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
1062	// Normalize the constraint and tighten the constant term by the GCD.
1063	DynamicAPInt gcd = inequalities.normalizeRow(row: i, nCols: getNumCols() - 1);
1064	if (gcd > 1)
1065	atIneq(i, j: numCols - 1) = floorDiv(LHS: atIneq(i, j: numCols - 1), RHS: gcd);
1066	}
1067	}
1068
1069	// Eliminates all variable variables in column range [posStart, posLimit).
1070	// Returns the number of variables eliminated.
1071	unsigned IntegerRelation::gaussianEliminateVars(unsigned posStart,
1072	unsigned posLimit) {
1073	// Return if variable positions to eliminate are out of range.
1074	assert(posLimit <= getNumVars());
1075	assert(hasConsistentState());
1076
1077	if (posStart >= posLimit)
1078	return 0;
1079
1080	gcdTightenInequalities();
1081
1082	unsigned pivotCol = 0;
1083	for (pivotCol = posStart; pivotCol < posLimit; ++pivotCol) {
1084	// Find a row which has a non-zero coefficient in column 'j'.
1085	std::optional<unsigned> pivotRow =
1086	findConstraintWithNonZeroAt(colIdx: pivotCol, /*isEq=*/true);
1087	// No pivot row in equalities with non-zero at 'pivotCol'.
1088	if (!pivotRow) {
1089	// If inequalities are also non-zero in 'pivotCol', it can be eliminated.
1090	if ((pivotRow = findConstraintWithNonZeroAt(colIdx: pivotCol, /*isEq=*/false)))
1091	break;
1092	continue;
1093	}
1094
1095	// Eliminate variable at 'pivotCol' from each equality row.
1096	for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
1097	eliminateFromConstraint(constraints: this, rowIdx: i, pivotRow: *pivotRow, pivotCol, elimColStart: posStart,
1098	/*isEq=*/true);
1099	equalities.normalizeRow(row: i);
1100	}
1101
1102	// Eliminate variable at 'pivotCol' from each inequality row.
1103	for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
1104	eliminateFromConstraint(constraints: this, rowIdx: i, pivotRow: *pivotRow, pivotCol, elimColStart: posStart,
1105	/*isEq=*/false);
1106	inequalities.normalizeRow(row: i);
1107	}
1108	removeEquality(pos: *pivotRow);
1109	gcdTightenInequalities();
1110	}
1111	// Update position limit based on number eliminated.
1112	posLimit = pivotCol;
1113	// Remove eliminated columns from all constraints.
1114	removeVarRange(varStart: posStart, varLimit: posLimit);
1115	return posLimit - posStart;
1116	}
1117
1118	bool IntegerRelation::gaussianEliminate() {
1119	gcdTightenInequalities();
1120	unsigned firstVar = 0, vars = getNumVars();
1121	unsigned nowDone, eqs;
1122	std::optional<unsigned> pivotRow;
1123	for (nowDone = 0, eqs = getNumEqualities(); nowDone < eqs; ++nowDone) {
1124	// Finds the first non-empty column.
1125	for (; firstVar < vars; ++firstVar) {
1126	if ((pivotRow = findConstraintWithNonZeroAt(colIdx: firstVar, /*isEq=*/true)))
1127	break;
1128	}
1129	// The matrix has been normalized to row echelon form.
1130	if (firstVar >= vars)
1131	break;
1132
1133	// The first pivot row found is below where it should currently be placed.
1134	if (*pivotRow > nowDone) {
1135	equalities.swapRows(row: *pivotRow, otherRow: nowDone);
1136	*pivotRow = nowDone;
1137	}
1138
1139	// Normalize all lower equations and all inequalities.
1140	for (unsigned i = nowDone + 1; i < eqs; ++i) {
1141	eliminateFromConstraint(constraints: this, rowIdx: i, pivotRow: *pivotRow, pivotCol: firstVar, elimColStart: 0, isEq: true);
1142	equalities.normalizeRow(row: i);
1143	}
1144	for (unsigned i = 0, ineqs = getNumInequalities(); i < ineqs; ++i) {
1145	eliminateFromConstraint(constraints: this, rowIdx: i, pivotRow: *pivotRow, pivotCol: firstVar, elimColStart: 0, isEq: false);
1146	inequalities.normalizeRow(row: i);
1147	}
1148	gcdTightenInequalities();
1149	}
1150
1151	// No redundant rows.
1152	if (nowDone == eqs)
1153	return false;
1154
1155	// Check to see if the redundant rows constant is zero, a non-zero value means
1156	// the set is empty.
1157	for (unsigned i = nowDone; i < eqs; ++i) {
1158	if (atEq(i, j: vars) == 0)
1159	continue;
1160
1161	*this = getEmpty(space: getSpace());
1162	return true;
1163	}
1164	// Eliminate rows that are confined to be all zeros.
1165	removeEqualityRange(start: nowDone, end: eqs);
1166	return true;
1167	}
1168
1169	// A more complex check to eliminate redundant inequalities. Uses FourierMotzkin
1170	// to check if a constraint is redundant.
1171	void IntegerRelation::removeRedundantInequalities() {
1172	SmallVector<bool, 32> redun(getNumInequalities(), false);
1173	// To check if an inequality is redundant, we replace the inequality by its
1174	// complement (for eg., i - 1 >= 0 by i <= 0), and check if the resulting
1175	// system is empty. If it is, the inequality is redundant.
1176	IntegerRelation tmpCst(*this);
1177	for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
1178	// Change the inequality to its complement.
1179	tmpCst.inequalities.negateRow(row: r);
1180	--tmpCst.atIneq(i: r, j: tmpCst.getNumCols() - 1);
1181	if (tmpCst.isEmpty()) {
1182	redun[r] = true;
1183	// Zero fill the redundant inequality.
1184	inequalities.fillRow(row: r, /*value=*/0);
1185	tmpCst.inequalities.fillRow(row: r, /*value=*/0);
1186	} else {
1187	// Reverse the change (to avoid recreating tmpCst each time).
1188	++tmpCst.atIneq(i: r, j: tmpCst.getNumCols() - 1);
1189	tmpCst.inequalities.negateRow(row: r);
1190	}
1191	}
1192
1193	unsigned pos = 0;
1194	for (unsigned r = 0, e = getNumInequalities(); r < e; ++r) {
1195	if (!redun[r])
1196	inequalities.copyRow(sourceRow: r, targetRow: pos++);
1197	}
1198	inequalities.resizeVertically(newNRows: pos);
1199	}
1200
1201	// A more complex check to eliminate redundant inequalities and equalities. Uses
1202	// Simplex to check if a constraint is redundant.
1203	void IntegerRelation::removeRedundantConstraints() {
1204	// First, we run gcdTightenInequalities. This allows us to catch some
1205	// constraints which are not redundant when considering rational solutions
1206	// but are redundant in terms of integer solutions.
1207	gcdTightenInequalities();
1208	Simplex simplex(*this);
1209	simplex.detectRedundant();
1210
1211	unsigned pos = 0;
1212	unsigned numIneqs = getNumInequalities();
1213	// Scan to get rid of all inequalities marked redundant, in-place. In Simplex,
1214	// the first constraints added are the inequalities.
1215	for (unsigned r = 0; r < numIneqs; r++) {
1216	if (!simplex.isMarkedRedundant(constraintIndex: r))
1217	inequalities.copyRow(sourceRow: r, targetRow: pos++);
1218	}
1219	inequalities.resizeVertically(newNRows: pos);
1220
1221	// Scan to get rid of all equalities marked redundant, in-place. In Simplex,
1222	// after the inequalities, a pair of constraints for each equality is added.
1223	// An equality is redundant if both the inequalities in its pair are
1224	// redundant.
1225	pos = 0;
1226	for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
1227	if (!(simplex.isMarkedRedundant(constraintIndex: numIneqs + 2 * r) &&
1228	simplex.isMarkedRedundant(constraintIndex: numIneqs + 2 * r + 1)))
1229	equalities.copyRow(sourceRow: r, targetRow: pos++);
1230	}
1231	equalities.resizeVertically(newNRows: pos);
1232	}
1233
1234	std::optional<DynamicAPInt> IntegerRelation::computeVolume() const {
1235	assert(getNumSymbolVars() == 0 && "Symbols are not yet supported!");
1236
1237	Simplex simplex(*this);
1238	// If the polytope is rationally empty, there are certainly no integer
1239	// points.
1240	if (simplex.isEmpty())
1241	return DynamicAPInt(0);
1242
1243	// Just find the maximum and minimum integer value of each non-local var
1244	// separately, thus finding the number of integer values each such var can
1245	// take. Multiplying these together gives a valid overapproximation of the
1246	// number of integer points in the relation. The result this gives is
1247	// equivalent to projecting (rationally) the relation onto its non-local vars
1248	// and returning the number of integer points in a minimal axis-parallel
1249	// hyperrectangular overapproximation of that.
1250	//
1251	// We also handle the special case where one dimension is unbounded and
1252	// another dimension can take no integer values. In this case, the volume is
1253	// zero.
1254	//
1255	// If there is no such empty dimension, if any dimension is unbounded we
1256	// just return the result as unbounded.
1257	DynamicAPInt count(1);
1258	SmallVector<DynamicAPInt, 8> dim(getNumVars() + 1);
1259	bool hasUnboundedVar = false;
1260	for (unsigned i = 0, e = getNumDimAndSymbolVars(); i < e; ++i) {
1261	dim[i] = 1;
1262	auto [min, max] = simplex.computeIntegerBounds(coeffs: dim);
1263	dim[i] = 0;
1264
1265	assert((!min.isEmpty() && !max.isEmpty()) &&
1266	"Polytope should be rationally non-empty!");
1267
1268	// One of the dimensions is unbounded. Note this fact. We will return
1269	// unbounded if none of the other dimensions makes the volume zero.
1270	if (min.isUnbounded() || max.isUnbounded()) {
1271	hasUnboundedVar = true;
1272	continue;
1273	}
1274
1275	// In this case there are no valid integer points and the volume is
1276	// definitely zero.
1277	if (min.getBoundedOptimum() > max.getBoundedOptimum())
1278	return DynamicAPInt(0);
1279
1280	count *= (*max - *min + 1);
1281	}
1282
1283	if (count == 0)
1284	return DynamicAPInt(0);
1285	if (hasUnboundedVar)
1286	return {};
1287	return count;
1288	}
1289
1290	void IntegerRelation::eliminateRedundantLocalVar(unsigned posA, unsigned posB) {
1291	assert(posA < getNumLocalVars() && "Invalid local var position");
1292	assert(posB < getNumLocalVars() && "Invalid local var position");
1293
1294	unsigned localOffset = getVarKindOffset(kind: VarKind::Local);
1295	posA += localOffset;
1296	posB += localOffset;
1297	inequalities.addToColumn(sourceColumn: posB, targetColumn: posA, scale: 1);
1298	equalities.addToColumn(sourceColumn: posB, targetColumn: posA, scale: 1);
1299	removeVar(pos: posB);
1300	}
1301
1302	/// mergeAndAlignSymbols's implementation can be broken down into two steps:
1303	/// 1. Merge and align identifiers into `other` from `this. If an identifier
1304	/// from `this` exists in `other` then we align it. Otherwise, we assume it is a
1305	/// new identifier and insert it into `other` in the same position as `this`.
1306	/// 2. Add identifiers that are in `other` but not `this to `this`.
1307	void IntegerRelation::mergeAndAlignSymbols(IntegerRelation &other) {
1308	assert(space.isUsingIds() && other.space.isUsingIds() &&
1309	"both relations need to have identifers to merge and align");
1310
1311	unsigned i = 0;
1312	for (const Identifier identifier : space.getIds(kind: VarKind::Symbol)) {
1313	// Search in `other` starting at position `i` since the left of `i` is
1314	// aligned.
1315	const Identifier *findBegin =
1316	other.space.getIds(kind: VarKind::Symbol).begin() + i;
1317	const Identifier *findEnd = other.space.getIds(kind: VarKind::Symbol).end();
1318	const Identifier *itr = std::find(first: findBegin, last: findEnd, val: identifier);
1319	if (itr != findEnd) {
1320	other.swapVar(posA: other.getVarKindOffset(kind: VarKind::Symbol) + i,
1321	posB: other.getVarKindOffset(kind: VarKind::Symbol) + i +
1322	std::distance(first: findBegin, last: itr));
1323	} else {
1324	other.insertVar(kind: VarKind::Symbol, pos: i);
1325	other.space.setId(kind: VarKind::Symbol, pos: i, id: identifier);
1326	}
1327	++i;
1328	}
1329
1330	for (unsigned e = other.getNumVarKind(kind: VarKind::Symbol); i < e; ++i) {
1331	insertVar(kind: VarKind::Symbol, pos: i);
1332	space.setId(kind: VarKind::Symbol, pos: i, id: other.space.getId(kind: VarKind::Symbol, pos: i));
1333	}
1334	}
1335
1336	/// Adds additional local ids to the sets such that they both have the union
1337	/// of the local ids in each set, without changing the set of points that
1338	/// lie in `this` and `other`.
1339	///
1340	/// To detect local ids that always take the same value, each local id is
1341	/// represented as a floordiv with constant denominator in terms of other ids.
1342	/// After extracting these divisions, local ids in `other` with the same
1343	/// division representation as some other local id in any set are considered
1344	/// duplicate and are merged.
1345	///
1346	/// It is possible that division representation for some local id cannot be
1347	/// obtained, and thus these local ids are not considered for detecting
1348	/// duplicates.
1349	unsigned IntegerRelation::mergeLocalVars(IntegerRelation &other) {
1350	IntegerRelation &relA = *this;
1351	IntegerRelation &relB = other;
1352
1353	unsigned oldALocals = relA.getNumLocalVars();
1354
1355	// Merge function that merges the local variables in both sets by treating
1356	// them as the same variable.
1357	auto merge = [&relA, &relB, oldALocals](unsigned i, unsigned j) -> bool {
1358	// We only merge from local at pos j to local at pos i, where j > i.
1359	if (i >= j)
1360	return false;
1361
1362	// If i < oldALocals, we are trying to merge duplicate divs. Since we do not
1363	// want to merge duplicates in A, we ignore this call.
1364	if (j < oldALocals)
1365	return false;
1366
1367	// Merge local at pos j into local at position i.
1368	relA.eliminateRedundantLocalVar(posA: i, posB: j);
1369	relB.eliminateRedundantLocalVar(posA: i, posB: j);
1370	return true;
1371	};
1372
1373	presburger::mergeLocalVars(relA&: *this, relB&: other, merge);
1374
1375	// Since we do not remove duplicate divisions in relA, this is guranteed to be
1376	// non-negative.
1377	return relA.getNumLocalVars() - oldALocals;
1378	}
1379
1380	bool IntegerRelation::hasOnlyDivLocals() const {
1381	return getLocalReprs().hasAllReprs();
1382	}
1383
1384	void IntegerRelation::removeDuplicateDivs() {
1385	DivisionRepr divs = getLocalReprs();
1386	auto merge = [this](unsigned i, unsigned j) -> bool {
1387	eliminateRedundantLocalVar(posA: i, posB: j);
1388	return true;
1389	};
1390	divs.removeDuplicateDivs(merge);
1391	}
1392
1393	void IntegerRelation::simplify() {
1394	bool changed = true;
1395	// Repeat until we reach a fixed point.
1396	while (changed) {
1397	if (isObviouslyEmpty())
1398	return;
1399	changed = false;
1400	normalizeConstraintsByGCD();
1401	changed |= gaussianEliminate();
1402	changed |= removeDuplicateConstraints();
1403	}
1404	// Current set is not empty.
1405	}
1406
1407	/// Removes local variables using equalities. Each equality is checked if it
1408	/// can be reduced to the form: `e = affine-expr`, where `e` is a local
1409	/// variable and `affine-expr` is an affine expression not containing `e`.
1410	/// If an equality satisfies this form, the local variable is replaced in
1411	/// each constraint and then removed. The equality used to replace this local
1412	/// variable is also removed.
1413	void IntegerRelation::removeRedundantLocalVars() {
1414	// Normalize the equality constraints to reduce coefficients of local
1415	// variables to 1 wherever possible.
1416	for (unsigned i = 0, e = getNumEqualities(); i < e; ++i)
1417	equalities.normalizeRow(row: i);
1418
1419	while (true) {
1420	unsigned i, e, j, f;
1421	for (i = 0, e = getNumEqualities(); i < e; ++i) {
1422	// Find a local variable to eliminate using ith equality.
1423	for (j = getNumDimAndSymbolVars(), f = getNumVars(); j < f; ++j)
1424	if (abs(X: atEq(i, j)) == 1)
1425	break;
1426
1427	// Local variable can be eliminated using ith equality.
1428	if (j < f)
1429	break;
1430	}
1431
1432	// No equality can be used to eliminate a local variable.
1433	if (i == e)
1434	break;
1435
1436	// Use the ith equality to simplify other equalities. If any changes
1437	// are made to an equality constraint, it is normalized by GCD.
1438	for (unsigned k = 0, t = getNumEqualities(); k < t; ++k) {
1439	if (atEq(i: k, j) != 0) {
1440	eliminateFromConstraint(constraints: this, rowIdx: k, pivotRow: i, pivotCol: j, elimColStart: j, /*isEq=*/true);
1441	equalities.normalizeRow(row: k);
1442	}
1443	}
1444
1445	// Use the ith equality to simplify inequalities.
1446	for (unsigned k = 0, t = getNumInequalities(); k < t; ++k)
1447	eliminateFromConstraint(constraints: this, rowIdx: k, pivotRow: i, pivotCol: j, elimColStart: j, /*isEq=*/false);
1448
1449	// Remove the ith equality and the found local variable.
1450	removeVar(pos: j);
1451	removeEquality(pos: i);
1452	}
1453	}
1454
1455	void IntegerRelation::convertVarKind(VarKind srcKind, unsigned varStart,
1456	unsigned varLimit, VarKind dstKind,
1457	unsigned pos) {
1458	assert(varLimit <= getNumVarKind(srcKind) && "invalid id range");
1459
1460	if (varStart >= varLimit)
1461	return;
1462
1463	unsigned srcOffset = getVarKindOffset(kind: srcKind);
1464	unsigned dstOffset = getVarKindOffset(kind: dstKind);
1465	unsigned convertCount = varLimit - varStart;
1466	int forwardMoveOffset = dstOffset > srcOffset ? -convertCount : 0;
1467
1468	equalities.moveColumns(srcPos: srcOffset + varStart, num: convertCount,
1469	dstPos: dstOffset + pos + forwardMoveOffset);
1470	inequalities.moveColumns(srcPos: srcOffset + varStart, num: convertCount,
1471	dstPos: dstOffset + pos + forwardMoveOffset);
1472
1473	space.convertVarKind(srcKind, srcPos: varStart, num: varLimit - varStart, dstKind, dstPos: pos);
1474	}
1475
1476	void IntegerRelation::addBound(BoundType type, unsigned pos,
1477	const DynamicAPInt &value) {
1478	assert(pos < getNumCols());
1479	if (type == BoundType::EQ) {
1480	unsigned row = equalities.appendExtraRow();
1481	equalities(row, pos) = 1;
1482	equalities(row, getNumCols() - 1) = -value;
1483	} else {
1484	unsigned row = inequalities.appendExtraRow();
1485	inequalities(row, pos) = type == BoundType::LB ? 1 : -1;
1486	inequalities(row, getNumCols() - 1) =
1487	type == BoundType::LB ? -value : value;
1488	}
1489	}
1490
1491	void IntegerRelation::addBound(BoundType type, ArrayRef<DynamicAPInt> expr,
1492	const DynamicAPInt &value) {
1493	assert(type != BoundType::EQ && "EQ not implemented");
1494	assert(expr.size() == getNumCols());
1495	unsigned row = inequalities.appendExtraRow();
1496	for (unsigned i = 0, e = expr.size(); i < e; ++i)
1497	inequalities(row, i) = type == BoundType::LB ? expr[i] : -expr[i];
1498	inequalities(inequalities.getNumRows() - 1, getNumCols() - 1) +=
1499	type == BoundType::LB ? -value : value;
1500	}
1501
1502	/// Adds a new local variable as the floordiv of an affine function of other
1503	/// variables, the coefficients of which are provided in 'dividend' and with
1504	/// respect to a positive constant 'divisor'. Two constraints are added to the
1505	/// system to capture equivalence with the floordiv.
1506	/// q = expr floordiv c <=> c*q <= expr <= c*q + c - 1.
1507	void IntegerRelation::addLocalFloorDiv(ArrayRef<DynamicAPInt> dividend,
1508	const DynamicAPInt &divisor) {
1509	assert(dividend.size() == getNumCols() && "incorrect dividend size");
1510	assert(divisor > 0 && "positive divisor expected");
1511
1512	appendVar(kind: VarKind::Local);
1513
1514	SmallVector<DynamicAPInt, 8> dividendCopy(dividend);
1515	dividendCopy.insert(I: dividendCopy.end() - 1, Elt: DynamicAPInt(0));
1516	addInequality(
1517	inEq: getDivLowerBound(dividend: dividendCopy, divisor, localVarIdx: dividendCopy.size() - 2));
1518	addInequality(
1519	inEq: getDivUpperBound(dividend: dividendCopy, divisor, localVarIdx: dividendCopy.size() - 2));
1520	}
1521
1522	int IntegerRelation::findEqualityToConstant(unsigned pos, bool symbolic) const {
1523	assert(pos < getNumVars() && "invalid position");
1524	for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
1525	DynamicAPInt v = atEq(i: r, j: pos);
1526	if (v * v != 1)
1527	continue;
1528	unsigned c;
1529	unsigned f = symbolic ? getNumDimVars() : getNumVars();
1530	// This checks for zeros in all positions other than 'pos' in [0, f)
1531	for (c = 0; c < f; c++) {
1532	if (c == pos)
1533	continue;
1534	if (atEq(i: r, j: c) != 0) {
1535	// Dependent on another variable.
1536	break;
1537	}
1538	}
1539	if (c == f)
1540	// Equality is free of other variables.
1541	return r;
1542	}
1543	return -1;
1544	}
1545
1546	LogicalResult IntegerRelation::constantFoldVar(unsigned pos) {
1547	assert(pos < getNumVars() && "invalid position");
1548	int rowIdx;
1549	if ((rowIdx = findEqualityToConstant(pos)) == -1)
1550	return failure();
1551
1552	// atEq(rowIdx, pos) is either -1 or 1.
1553	assert(atEq(rowIdx, pos) * atEq(rowIdx, pos) == 1);
1554	DynamicAPInt constVal = -atEq(i: rowIdx, j: getNumCols() - 1) / atEq(i: rowIdx, j: pos);
1555	setAndEliminate(pos, values: constVal);
1556	return success();
1557	}
1558
1559	void IntegerRelation::constantFoldVarRange(unsigned pos, unsigned num) {
1560	for (unsigned s = pos, t = pos, e = pos + num; s < e; s++) {
1561	if (constantFoldVar(pos: t).failed())
1562	t++;
1563	}
1564	}
1565
1566	/// Returns a non-negative constant bound on the extent (upper bound - lower
1567	/// bound) of the specified variable if it is found to be a constant; returns
1568	/// std::nullopt if it's not a constant. This methods treats symbolic variables
1569	/// specially, i.e., it looks for constant differences between affine
1570	/// expressions involving only the symbolic variables. See comments at function
1571	/// definition for example. 'lb', if provided, is set to the lower bound
1572	/// associated with the constant difference. Note that 'lb' is purely symbolic
1573	/// and thus will contain the coefficients of the symbolic variables and the
1574	/// constant coefficient.
1575	// Egs: 0 <= i <= 15, return 16.
1576	// s0 + 2 <= i <= s0 + 17, returns 16. (s0 has to be a symbol)
1577	// s0 + s1 + 16 <= d0 <= s0 + s1 + 31, returns 16.
1578	// s0 - 7 <= 8*j <= s0 returns 1 with lb = s0, lbDivisor = 8 (since lb =
1579	// ceil(s0 - 7 / 8) = floor(s0 / 8)).
1580	std::optional<DynamicAPInt> IntegerRelation::getConstantBoundOnDimSize(
1581	unsigned pos, SmallVectorImpl<DynamicAPInt> *lb,
1582	DynamicAPInt *boundFloorDivisor, SmallVectorImpl<DynamicAPInt> *ub,
1583	unsigned *minLbPos, unsigned *minUbPos) const {
1584	assert(pos < getNumDimVars() && "Invalid variable position");
1585
1586	// Find an equality for 'pos'^th variable that equates it to some function
1587	// of the symbolic variables (+ constant).
1588	int eqPos = findEqualityToConstant(pos, /*symbolic=*/true);
1589	if (eqPos != -1) {
1590	auto eq = getEquality(idx: eqPos);
1591	// If the equality involves a local var, we do not handle it.
1592	// FlatLinearConstraints can instead be used to detect the local variable as
1593	// an affine function (potentially div/mod) of other variables and use
1594	// affine expressions/maps to represent output.
1595	if (!std::all_of(first: eq.begin() + getNumDimAndSymbolVars(), last: eq.end() - 1,
1596	pred: [](const DynamicAPInt &coeff) { return coeff == 0; }))
1597	return std::nullopt;
1598
1599	// This variable can only take a single value.
1600	if (lb) {
1601	// Set lb to that symbolic value.
1602	lb->resize(N: getNumSymbolVars() + 1);
1603	if (ub)
1604	ub->resize(N: getNumSymbolVars() + 1);
1605	for (unsigned c = 0, f = getNumSymbolVars() + 1; c < f; c++) {
1606	DynamicAPInt v = atEq(i: eqPos, j: pos);
1607	// atEq(eqRow, pos) is either -1 or 1.
1608	assert(v * v == 1);
1609	(*lb)[c] = v < 0 ? atEq(i: eqPos, j: getNumDimVars() + c) / -v
1610	: -atEq(i: eqPos, j: getNumDimVars() + c) / v;
1611	// Since this is an equality, ub = lb.
1612	if (ub)
1613	(*ub)[c] = (*lb)[c];
1614	}
1615	assert(boundFloorDivisor &&
1616	"both lb and divisor or none should be provided");
1617	*boundFloorDivisor = 1;
1618	}
1619	if (minLbPos)
1620	*minLbPos = eqPos;
1621	if (minUbPos)
1622	*minUbPos = eqPos;
1623	return DynamicAPInt(1);
1624	}
1625
1626	// Check if the variable appears at all in any of the inequalities.
1627	unsigned r, e;
1628	for (r = 0, e = getNumInequalities(); r < e; r++) {
1629	if (atIneq(i: r, j: pos) != 0)
1630	break;
1631	}
1632	if (r == e)
1633	// If it doesn't, there isn't a bound on it.
1634	return std::nullopt;
1635
1636	// Positions of constraints that are lower/upper bounds on the variable.
1637	SmallVector<unsigned, 4> lbIndices, ubIndices;
1638
1639	// Gather all symbolic lower bounds and upper bounds of the variable, i.e.,
1640	// the bounds can only involve symbolic (and local) variables. Since the
1641	// canonical form c_1*x_1 + c_2*x_2 + ... + c_0 >= 0, a constraint is a lower
1642	// bound for x_i if c_i >= 1, and an upper bound if c_i <= -1.
1643	getLowerAndUpperBoundIndices(pos, lbIndices: &lbIndices, ubIndices: &ubIndices,
1644	/*eqIndices=*/nullptr, /*offset=*/0,
1645	/*num=*/getNumDimVars());
1646
1647	std::optional<DynamicAPInt> minDiff;
1648	unsigned minLbPosition = 0, minUbPosition = 0;
1649	for (auto ubPos : ubIndices) {
1650	for (auto lbPos : lbIndices) {
1651	// Look for a lower bound and an upper bound that only differ by a
1652	// constant, i.e., pairs of the form 0 <= c_pos - f(c_i's) <= diffConst.
1653	// For example, if ii is the pos^th variable, we are looking for
1654	// constraints like ii >= i, ii <= ii + 50, 50 being the difference. The
1655	// minimum among all such constant differences is kept since that's the
1656	// constant bounding the extent of the pos^th variable.
1657	unsigned j, e;
1658	for (j = 0, e = getNumCols() - 1; j < e; j++)
1659	if (atIneq(i: ubPos, j) != -atIneq(i: lbPos, j)) {
1660	break;
1661	}
1662	if (j < getNumCols() - 1)
1663	continue;
1664	DynamicAPInt diff = ceilDiv(LHS: atIneq(i: ubPos, j: getNumCols() - 1) +
1665	atIneq(i: lbPos, j: getNumCols() - 1) + 1,
1666	RHS: atIneq(i: lbPos, j: pos));
1667	// This bound is non-negative by definition.
1668	diff = std::max<DynamicAPInt>(a: diff, b: DynamicAPInt(0));
1669	if (minDiff == std::nullopt || diff < minDiff) {
1670	minDiff = diff;
1671	minLbPosition = lbPos;
1672	minUbPosition = ubPos;
1673	}
1674	}
1675	}
1676	if (lb && minDiff) {
1677	// Set lb to the symbolic lower bound.
1678	lb->resize(N: getNumSymbolVars() + 1);
1679	if (ub)
1680	ub->resize(N: getNumSymbolVars() + 1);
1681	// The lower bound is the ceildiv of the lb constraint over the coefficient
1682	// of the variable at 'pos'. We express the ceildiv equivalently as a floor
1683	// for uniformity. For eg., if the lower bound constraint was: 32*d0 - N +
1684	// 31 >= 0, the lower bound for d0 is ceil(N - 31, 32), i.e., floor(N, 32).
1685	*boundFloorDivisor = atIneq(i: minLbPosition, j: pos);
1686	assert(*boundFloorDivisor == -atIneq(minUbPosition, pos));
1687	for (unsigned c = 0, e = getNumSymbolVars() + 1; c < e; c++) {
1688	(*lb)[c] = -atIneq(i: minLbPosition, j: getNumDimVars() + c);
1689	}
1690	if (ub) {
1691	for (unsigned c = 0, e = getNumSymbolVars() + 1; c < e; c++)
1692	(*ub)[c] = atIneq(i: minUbPosition, j: getNumDimVars() + c);
1693	}
1694	// The lower bound leads to a ceildiv while the upper bound is a floordiv
1695	// whenever the coefficient at pos != 1. ceildiv (val / d) = floordiv (val +
1696	// d - 1 / d); hence, the addition of 'atIneq(minLbPosition, pos) - 1' to
1697	// the constant term for the lower bound.
1698	(*lb)[getNumSymbolVars()] += atIneq(i: minLbPosition, j: pos) - 1;
1699	}
1700	if (minLbPos)
1701	*minLbPos = minLbPosition;
1702	if (minUbPos)
1703	*minUbPos = minUbPosition;
1704	return minDiff;
1705	}
1706
1707	template <bool isLower>
1708	std::optional<DynamicAPInt>
1709	IntegerRelation::computeConstantLowerOrUpperBound(unsigned pos) {
1710	assert(pos < getNumVars() && "invalid position");
1711	// Project to 'pos'.
1712	projectOut(pos: 0, num: pos);
1713	projectOut(pos: 1, num: getNumVars() - 1);
1714	// Check if there's an equality equating the '0'^th variable to a constant.
1715	int eqRowIdx = findEqualityToConstant(/*pos=*/0, /*symbolic=*/false);
1716	if (eqRowIdx != -1)
1717	// atEq(rowIdx, 0) is either -1 or 1.
1718	return -atEq(i: eqRowIdx, j: getNumCols() - 1) / atEq(i: eqRowIdx, j: 0);
1719
1720	// Check if the variable appears at all in any of the inequalities.
1721	unsigned r, e;
1722	for (r = 0, e = getNumInequalities(); r < e; r++) {
1723	if (atIneq(i: r, j: 0) != 0)
1724	break;
1725	}
1726	if (r == e)
1727	// If it doesn't, there isn't a bound on it.
1728	return std::nullopt;
1729
1730	std::optional<DynamicAPInt> minOrMaxConst;
1731
1732	// Take the max across all const lower bounds (or min across all constant
1733	// upper bounds).
1734	for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
1735	if (isLower) {
1736	if (atIneq(i: r, j: 0) <= 0)
1737	// Not a lower bound.
1738	continue;
1739	} else if (atIneq(i: r, j: 0) >= 0) {
1740	// Not an upper bound.
1741	continue;
1742	}
1743	unsigned c, f;
1744	for (c = 0, f = getNumCols() - 1; c < f; c++)
1745	if (c != 0 && atIneq(i: r, j: c) != 0)
1746	break;
1747	if (c < getNumCols() - 1)
1748	// Not a constant bound.
1749	continue;
1750
1751	DynamicAPInt boundConst =
1752	isLower ? ceilDiv(LHS: -atIneq(i: r, j: getNumCols() - 1), RHS: atIneq(i: r, j: 0))
1753	: floorDiv(LHS: atIneq(i: r, j: getNumCols() - 1), RHS: -atIneq(i: r, j: 0));
1754	if (isLower) {
1755	if (minOrMaxConst == std::nullopt || boundConst > minOrMaxConst)
1756	minOrMaxConst = boundConst;
1757	} else {
1758	if (minOrMaxConst == std::nullopt || boundConst < minOrMaxConst)
1759	minOrMaxConst = boundConst;
1760	}
1761	}
1762	return minOrMaxConst;
1763	}
1764
1765	std::optional<DynamicAPInt>
1766	IntegerRelation::getConstantBound(BoundType type, unsigned pos) const {
1767	if (type == BoundType::LB)
1768	return IntegerRelation(*this)
1769	.computeConstantLowerOrUpperBound</*isLower=*/true>(pos);
1770	if (type == BoundType::UB)
1771	return IntegerRelation(*this)
1772	.computeConstantLowerOrUpperBound</*isLower=*/false>(pos);
1773
1774	assert(type == BoundType::EQ && "expected EQ");
1775	std::optional<DynamicAPInt> lb =
1776	IntegerRelation(*this).computeConstantLowerOrUpperBound</*isLower=*/true>(
1777	pos);
1778	std::optional<DynamicAPInt> ub =
1779	IntegerRelation(*this)
1780	.computeConstantLowerOrUpperBound</*isLower=*/false>(pos);
1781	return (lb && ub && *lb == *ub) ? std::optional<DynamicAPInt>(*ub)
1782	: std::nullopt;
1783	}
1784
1785	// A simple (naive and conservative) check for hyper-rectangularity.
1786	bool IntegerRelation::isHyperRectangular(unsigned pos, unsigned num) const {
1787	assert(pos < getNumCols() - 1);
1788	// Check for two non-zero coefficients in the range [pos, pos + sum).
1789	for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
1790	unsigned sum = 0;
1791	for (unsigned c = pos; c < pos + num; c++) {
1792	if (atIneq(i: r, j: c) != 0)
1793	sum++;
1794	}
1795	if (sum > 1)
1796	return false;
1797	}
1798	for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
1799	unsigned sum = 0;
1800	for (unsigned c = pos; c < pos + num; c++) {
1801	if (atEq(i: r, j: c) != 0)
1802	sum++;
1803	}
1804	if (sum > 1)
1805	return false;
1806	}
1807	return true;
1808	}
1809
1810	/// Removes duplicate constraints, trivially true constraints, and constraints
1811	/// that can be detected as redundant as a result of differing only in their
1812	/// constant term part. A constraint of the form <non-negative constant> >= 0 is
1813	/// considered trivially true.
1814	// Uses a DenseSet to hash and detect duplicates followed by a linear scan to
1815	// remove duplicates in place.
1816	void IntegerRelation::removeTrivialRedundancy() {
1817	gcdTightenInequalities();
1818	normalizeConstraintsByGCD();
1819
1820	// A map used to detect redundancy stemming from constraints that only differ
1821	// in their constant term. The value stored is <row position, const term>
1822	// for a given row.
1823	SmallDenseMap<ArrayRef<DynamicAPInt>, std::pair<unsigned, DynamicAPInt>>
1824	rowsWithoutConstTerm;
1825
1826	// Check if constraint is of the form <non-negative-constant> >= 0.
1827	auto isTriviallyValid = [&](unsigned r) -> bool {
1828	for (unsigned c = 0, e = getNumCols() - 1; c < e; c++) {
1829	if (atIneq(i: r, j: c) != 0)
1830	return false;
1831	}
1832	return atIneq(i: r, j: getNumCols() - 1) >= 0;
1833	};
1834
1835	// Detect and mark redundant constraints.
1836	SmallVector<bool, 256> redunIneq(getNumInequalities(), false);
1837	for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
1838	DynamicAPInt *rowStart = &inequalities(r, 0);
1839	if (isTriviallyValid(r)) {
1840	redunIneq[r] = true;
1841	continue;
1842	}
1843
1844	// Among constraints that only differ in the constant term part, mark
1845	// everything other than the one with the smallest constant term redundant.
1846	// (eg: among i - 16j - 5 >= 0, i - 16j - 1 >=0, i - 16j - 7 >= 0, the
1847	// former two are redundant).
1848	DynamicAPInt constTerm = atIneq(i: r, j: getNumCols() - 1);
1849	auto rowWithoutConstTerm =
1850	ArrayRef<DynamicAPInt>(rowStart, getNumCols() - 1);
1851	const auto &ret =
1852	rowsWithoutConstTerm.insert(KV: {rowWithoutConstTerm, {r, constTerm}});
1853	if (!ret.second) {
1854	// Check if the other constraint has a higher constant term.
1855	auto &val = ret.first->second;
1856	if (val.second > constTerm) {
1857	// The stored row is redundant. Mark it so, and update with this one.
1858	redunIneq[val.first] = true;
1859	val = {r, constTerm};
1860	} else {
1861	// The one stored makes this one redundant.
1862	redunIneq[r] = true;
1863	}
1864	}
1865	}
1866
1867	// Scan to get rid of all rows marked redundant, in-place.
1868	unsigned pos = 0;
1869	for (unsigned r = 0, e = getNumInequalities(); r < e; r++)
1870	if (!redunIneq[r])
1871	inequalities.copyRow(sourceRow: r, targetRow: pos++);
1872
1873	inequalities.resizeVertically(newNRows: pos);
1874
1875	// TODO: consider doing this for equalities as well, but probably not worth
1876	// the savings.
1877	}
1878
1879	#undef DEBUG_TYPE
1880	#define DEBUG_TYPE "fm"
1881
1882	/// Eliminates variable at the specified position using Fourier-Motzkin
1883	/// variable elimination. This technique is exact for rational spaces but
1884	/// conservative (in "rare" cases) for integer spaces. The operation corresponds
1885	/// to a projection operation yielding the (convex) set of integer points
1886	/// contained in the rational shadow of the set. An emptiness test that relies
1887	/// on this method will guarantee emptiness, i.e., it disproves the existence of
1888	/// a solution if it says it's empty.
1889	/// If a non-null isResultIntegerExact is passed, it is set to true if the
1890	/// result is also integer exact. If it's set to false, the obtained solution
1891	/// *may* not be exact, i.e., it may contain integer points that do not have an
1892	/// integer pre-image in the original set.
1893	///
1894	/// Eg:
1895	/// j >= 0, j <= i + 1
1896	/// i >= 0, i <= N + 1
1897	/// Eliminating i yields,
1898	/// j >= 0, 0 <= N + 1, j - 1 <= N + 1
1899	///
1900	/// If darkShadow = true, this method computes the dark shadow on elimination;
1901	/// the dark shadow is a convex integer subset of the exact integer shadow. A
1902	/// non-empty dark shadow proves the existence of an integer solution. The
1903	/// elimination in such a case could however be an under-approximation, and thus
1904	/// should not be used for scanning sets or used by itself for dependence
1905	/// checking.
1906	///
1907	/// Eg: 2-d set, * represents grid points, 'o' represents a point in the set.
1908	/// ^
1909	/// |
1910	/// | * * * * o o
1911	/// i | * * o o o o
1912	/// | o * * * * *
1913	/// --------------->
1914	/// j ->
1915	///
1916	/// Eliminating i from this system (projecting on the j dimension):
1917	/// rational shadow / integer light shadow: 1 <= j <= 6
1918	/// dark shadow: 3 <= j <= 6
1919	/// exact integer shadow: j = 1 \union 3 <= j <= 6
1920	/// holes/splinters: j = 2
1921	///
1922	/// darkShadow = false, isResultIntegerExact = nullptr are default values.
1923	// TODO: a slight modification to yield dark shadow version of FM (tightened),
1924	// which can prove the existence of a solution if there is one.
1925	void IntegerRelation::fourierMotzkinEliminate(unsigned pos, bool darkShadow,
1926	bool *isResultIntegerExact) {
1927	LLVM_DEBUG(llvm::dbgs() << "FM input (eliminate pos " << pos << "):\n");
1928	LLVM_DEBUG(dump());
1929	assert(pos < getNumVars() && "invalid position");
1930	assert(hasConsistentState());
1931
1932	// Check if this variable can be eliminated through a substitution.
1933	for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
1934	if (atEq(i: r, j: pos) != 0) {
1935	// Use Gaussian elimination here (since we have an equality).
1936	LogicalResult ret = gaussianEliminateVar(position: pos);
1937	(void)ret;
1938	assert(ret.succeeded() && "Gaussian elimination guaranteed to succeed");
1939	LLVM_DEBUG(llvm::dbgs() << "FM output (through Gaussian elimination):\n");
1940	LLVM_DEBUG(dump());
1941	return;
1942	}
1943	}
1944
1945	// A fast linear time tightening.
1946	gcdTightenInequalities();
1947
1948	// Check if the variable appears at all in any of the inequalities.
1949	if (isColZero(pos)) {
1950	// If it doesn't appear, just remove the column and return.
1951	// TODO: refactor removeColumns to use it from here.
1952	removeVar(pos);
1953	LLVM_DEBUG(llvm::dbgs() << "FM output:\n");
1954	LLVM_DEBUG(dump());
1955	return;
1956	}
1957
1958	// Positions of constraints that are lower bounds on the variable.
1959	SmallVector<unsigned, 4> lbIndices;
1960	// Positions of constraints that are lower bounds on the variable.
1961	SmallVector<unsigned, 4> ubIndices;
1962	// Positions of constraints that do not involve the variable.
1963	std::vector<unsigned> nbIndices;
1964	nbIndices.reserve(n: getNumInequalities());
1965
1966	// Gather all lower bounds and upper bounds of the variable. Since the
1967	// canonical form c_1*x_1 + c_2*x_2 + ... + c_0 >= 0, a constraint is a lower
1968	// bound for x_i if c_i >= 1, and an upper bound if c_i <= -1.
1969	for (unsigned r = 0, e = getNumInequalities(); r < e; r++) {
1970	if (atIneq(i: r, j: pos) == 0) {
1971	// Var does not appear in bound.
1972	nbIndices.emplace_back(args&: r);
1973	} else if (atIneq(i: r, j: pos) >= 1) {
1974	// Lower bound.
1975	lbIndices.emplace_back(Args&: r);
1976	} else {
1977	// Upper bound.
1978	ubIndices.emplace_back(Args&: r);
1979	}
1980	}
1981
1982	PresburgerSpace newSpace = getSpace();
1983	VarKind idKindRemove = newSpace.getVarKindAt(pos);
1984	unsigned relativePos = pos - newSpace.getVarKindOffset(kind: idKindRemove);
1985	newSpace.removeVarRange(kind: idKindRemove, varStart: relativePos, varLimit: relativePos + 1);
1986
1987	/// Create the new system which has one variable less.
1988	IntegerRelation newRel(lbIndices.size() * ubIndices.size() + nbIndices.size(),
1989	getNumEqualities(), getNumCols() - 1, newSpace);
1990
1991	// This will be used to check if the elimination was integer exact.
1992	bool allLCMsAreOne = true;
1993
1994	// Let x be the variable we are eliminating.
1995	// For each lower bound, lb <= c_l*x, and each upper bound c_u*x <= ub, (note
1996	// that c_l, c_u >= 1) we have:
1997	// lb*lcm(c_l, c_u)/c_l <= lcm(c_l, c_u)*x <= ub*lcm(c_l, c_u)/c_u
1998	// We thus generate a constraint:
1999	// lcm(c_l, c_u)/c_l*lb <= lcm(c_l, c_u)/c_u*ub.
2000	// Note if c_l = c_u = 1, all integer points captured by the resulting
2001	// constraint correspond to integer points in the original system (i.e., they
2002	// have integer pre-images). Hence, if the lcm's are all 1, the elimination is
2003	// integer exact.
2004	for (auto ubPos : ubIndices) {
2005	for (auto lbPos : lbIndices) {
2006	SmallVector<DynamicAPInt, 4> ineq;
2007	ineq.reserve(N: newRel.getNumCols());
2008	DynamicAPInt lbCoeff = atIneq(i: lbPos, j: pos);
2009	// Note that in the comments above, ubCoeff is the negation of the
2010	// coefficient in the canonical form as the view taken here is that of the
2011	// term being moved to the other size of '>='.
2012	DynamicAPInt ubCoeff = -atIneq(i: ubPos, j: pos);
2013	// TODO: refactor this loop to avoid all branches inside.
2014	for (unsigned l = 0, e = getNumCols(); l < e; l++) {
2015	if (l == pos)
2016	continue;
2017	assert(lbCoeff >= 1 && ubCoeff >= 1 && "bounds wrongly identified");
2018	DynamicAPInt lcm = llvm::lcm(A: lbCoeff, B: ubCoeff);
2019	ineq.emplace_back(Args: atIneq(i: ubPos, j: l) * (lcm / ubCoeff) +
2020	atIneq(i: lbPos, j: l) * (lcm / lbCoeff));
2021	assert(lcm > 0 && "lcm should be positive!");
2022	if (lcm != 1)
2023	allLCMsAreOne = false;
2024	}
2025	if (darkShadow) {
2026	// The dark shadow is a convex subset of the exact integer shadow. If
2027	// there is a point here, it proves the existence of a solution.
2028	ineq[ineq.size() - 1] += lbCoeff * ubCoeff - lbCoeff - ubCoeff + 1;
2029	}
2030	// TODO: we need to have a way to add inequalities in-place in
2031	// IntegerRelation instead of creating and copying over.
2032	newRel.addInequality(inEq: ineq);
2033	}
2034	}
2035
2036	LLVM_DEBUG(llvm::dbgs() << "FM isResultIntegerExact: " << allLCMsAreOne
2037	<< "\n");
2038	if (allLCMsAreOne && isResultIntegerExact)
2039	*isResultIntegerExact = true;
2040
2041	// Copy over the constraints not involving this variable.
2042	for (auto nbPos : nbIndices) {
2043	SmallVector<DynamicAPInt, 4> ineq;
2044	ineq.reserve(N: getNumCols() - 1);
2045	for (unsigned l = 0, e = getNumCols(); l < e; l++) {
2046	if (l == pos)
2047	continue;
2048	ineq.emplace_back(Args&: atIneq(i: nbPos, j: l));
2049	}
2050	newRel.addInequality(inEq: ineq);
2051	}
2052
2053	assert(newRel.getNumConstraints() ==
2054	lbIndices.size() * ubIndices.size() + nbIndices.size());
2055
2056	// Copy over the equalities.
2057	for (unsigned r = 0, e = getNumEqualities(); r < e; r++) {
2058	SmallVector<DynamicAPInt, 4> eq;
2059	eq.reserve(N: newRel.getNumCols());
2060	for (unsigned l = 0, e = getNumCols(); l < e; l++) {
2061	if (l == pos)
2062	continue;
2063	eq.emplace_back(Args&: atEq(i: r, j: l));
2064	}
2065	newRel.addEquality(eq);
2066	}
2067
2068	// GCD tightening and normalization allows detection of more trivially
2069	// redundant constraints.
2070	newRel.gcdTightenInequalities();
2071	newRel.normalizeConstraintsByGCD();
2072	newRel.removeTrivialRedundancy();
2073	clearAndCopyFrom(other: newRel);
2074	LLVM_DEBUG(llvm::dbgs() << "FM output:\n");
2075	LLVM_DEBUG(dump());
2076	}
2077
2078	#undef DEBUG_TYPE
2079	#define DEBUG_TYPE "presburger"
2080
2081	void IntegerRelation::projectOut(unsigned pos, unsigned num) {
2082	if (num == 0)
2083	return;
2084
2085	// 'pos' can be at most getNumCols() - 2 if num > 0.
2086	assert((getNumCols() < 2 || pos <= getNumCols() - 2) && "invalid position");
2087	assert(pos + num < getNumCols() && "invalid range");
2088
2089	// Eliminate as many variables as possible using Gaussian elimination.
2090	unsigned currentPos = pos;
2091	unsigned numToEliminate = num;
2092	unsigned numGaussianEliminated = 0;
2093
2094	while (currentPos < getNumVars()) {
2095	unsigned curNumEliminated =
2096	gaussianEliminateVars(posStart: currentPos, posLimit: currentPos + numToEliminate);
2097	++currentPos;
2098	numToEliminate -= curNumEliminated + 1;
2099	numGaussianEliminated += curNumEliminated;
2100	}
2101
2102	// Eliminate the remaining using Fourier-Motzkin.
2103	for (unsigned i = 0; i < num - numGaussianEliminated; i++) {
2104	unsigned numToEliminate = num - numGaussianEliminated - i;
2105	fourierMotzkinEliminate(
2106	pos: getBestVarToEliminate(cst: *this, start: pos, end: pos + numToEliminate));
2107	}
2108
2109	// Fast/trivial simplifications.
2110	gcdTightenInequalities();
2111	// Normalize constraints after tightening since the latter impacts this, but
2112	// not the other way round.
2113	normalizeConstraintsByGCD();
2114	}
2115
2116	namespace {
2117
2118	enum BoundCmpResult { Greater, Less, Equal, Unknown };
2119
2120	/// Compares two affine bounds whose coefficients are provided in 'first' and
2121	/// 'second'. The last coefficient is the constant term.
2122	static BoundCmpResult compareBounds(ArrayRef<DynamicAPInt> a,
2123	ArrayRef<DynamicAPInt> b) {
2124	assert(a.size() == b.size());
2125
2126	// For the bounds to be comparable, their corresponding variable
2127	// coefficients should be equal; the constant terms are then compared to
2128	// determine less/greater/equal.
2129
2130	if (!std::equal(first1: a.begin(), last1: a.end() - 1, first2: b.begin()))
2131	return Unknown;
2132
2133	if (a.back() == b.back())
2134	return Equal;
2135
2136	return a.back() < b.back() ? Less : Greater;
2137	}
2138	} // namespace
2139
2140	// Returns constraints that are common to both A & B.
2141	static void getCommonConstraints(const IntegerRelation &a,
2142	const IntegerRelation &b, IntegerRelation &c) {
2143	c = IntegerRelation(a.getSpace());
2144	// a naive O(n^2) check should be enough here given the input sizes.
2145	for (unsigned r = 0, e = a.getNumInequalities(); r < e; ++r) {
2146	for (unsigned s = 0, f = b.getNumInequalities(); s < f; ++s) {
2147	if (a.getInequality(idx: r) == b.getInequality(idx: s)) {
2148	c.addInequality(inEq: a.getInequality(idx: r));
2149	break;
2150	}
2151	}
2152	}
2153	for (unsigned r = 0, e = a.getNumEqualities(); r < e; ++r) {
2154	for (unsigned s = 0, f = b.getNumEqualities(); s < f; ++s) {
2155	if (a.getEquality(idx: r) == b.getEquality(idx: s)) {
2156	c.addEquality(eq: a.getEquality(idx: r));
2157	break;
2158	}
2159	}
2160	}
2161	}
2162
2163	// Computes the bounding box with respect to 'other' by finding the min of the
2164	// lower bounds and the max of the upper bounds along each of the dimensions.
2165	LogicalResult
2166	IntegerRelation::unionBoundingBox(const IntegerRelation &otherCst) {
2167	assert(space.isEqual(otherCst.getSpace()) && "Spaces should match.");
2168	assert(getNumLocalVars() == 0 && "local ids not supported yet here");
2169
2170	// Get the constraints common to both systems; these will be added as is to
2171	// the union.
2172	IntegerRelation commonCst(PresburgerSpace::getRelationSpace());
2173	getCommonConstraints(a: *this, b: otherCst, c&: commonCst);
2174
2175	std::vector<SmallVector<DynamicAPInt, 8>> boundingLbs;
2176	std::vector<SmallVector<DynamicAPInt, 8>> boundingUbs;
2177	boundingLbs.reserve(n: 2 * getNumDimVars());
2178	boundingUbs.reserve(n: 2 * getNumDimVars());
2179
2180	// To hold lower and upper bounds for each dimension.
2181	SmallVector<DynamicAPInt, 4> lb, otherLb, ub, otherUb;
2182	// To compute min of lower bounds and max of upper bounds for each dimension.
2183	SmallVector<DynamicAPInt, 4> minLb(getNumSymbolVars() + 1);
2184	SmallVector<DynamicAPInt, 4> maxUb(getNumSymbolVars() + 1);
2185	// To compute final new lower and upper bounds for the union.
2186	SmallVector<DynamicAPInt, 8> newLb(getNumCols()), newUb(getNumCols());
2187
2188	DynamicAPInt lbFloorDivisor, otherLbFloorDivisor;
2189	for (unsigned d = 0, e = getNumDimVars(); d < e; ++d) {
2190	auto extent = getConstantBoundOnDimSize(pos: d, lb: &lb, boundFloorDivisor: &lbFloorDivisor, ub: &ub);
2191	if (!extent.has_value())
2192	// TODO: symbolic extents when necessary.
2193	// TODO: handle union if a dimension is unbounded.
2194	return failure();
2195
2196	auto otherExtent = otherCst.getConstantBoundOnDimSize(
2197	pos: d, lb: &otherLb, boundFloorDivisor: &otherLbFloorDivisor, ub: &otherUb);
2198	if (!otherExtent.has_value() || lbFloorDivisor != otherLbFloorDivisor)
2199	// TODO: symbolic extents when necessary.
2200	return failure();
2201
2202	assert(lbFloorDivisor > 0 && "divisor always expected to be positive");
2203
2204	auto res = compareBounds(a: lb, b: otherLb);
2205	// Identify min.
2206	if (res == BoundCmpResult::Less || res == BoundCmpResult::Equal) {
2207	minLb = lb;
2208	// Since the divisor is for a floordiv, we need to convert to ceildiv,
2209	// i.e., i >= expr floordiv div <=> i >= (expr - div + 1) ceildiv div <=>
2210	// div * i >= expr - div + 1.
2211	minLb.back() -= lbFloorDivisor - 1;
2212	} else if (res == BoundCmpResult::Greater) {
2213	minLb = otherLb;
2214	minLb.back() -= otherLbFloorDivisor - 1;
2215	} else {
2216	// Uncomparable - check for constant lower/upper bounds.
2217	auto constLb = getConstantBound(type: BoundType::LB, pos: d);
2218	auto constOtherLb = otherCst.getConstantBound(type: BoundType::LB, pos: d);
2219	if (!constLb.has_value() || !constOtherLb.has_value())
2220	return failure();
2221	std::fill(first: minLb.begin(), last: minLb.end(), value: 0);
2222	minLb.back() = std::min(a: *constLb, b: *constOtherLb);
2223	}
2224
2225	// Do the same for ub's but max of upper bounds. Identify max.
2226	auto uRes = compareBounds(a: ub, b: otherUb);
2227	if (uRes == BoundCmpResult::Greater || uRes == BoundCmpResult::Equal) {
2228	maxUb = ub;
2229	} else if (uRes == BoundCmpResult::Less) {
2230	maxUb = otherUb;
2231	} else {
2232	// Uncomparable - check for constant lower/upper bounds.
2233	auto constUb = getConstantBound(type: BoundType::UB, pos: d);
2234	auto constOtherUb = otherCst.getConstantBound(type: BoundType::UB, pos: d);
2235	if (!constUb.has_value() || !constOtherUb.has_value())
2236	return failure();
2237	std::fill(first: maxUb.begin(), last: maxUb.end(), value: 0);
2238	maxUb.back() = std::max(a: *constUb, b: *constOtherUb);
2239	}
2240
2241	std::fill(first: newLb.begin(), last: newLb.end(), value: 0);
2242	std::fill(first: newUb.begin(), last: newUb.end(), value: 0);
2243
2244	// The divisor for lb, ub, otherLb, otherUb at this point is lbDivisor,
2245	// and so it's the divisor for newLb and newUb as well.
2246	newLb[d] = lbFloorDivisor;
2247	newUb[d] = -lbFloorDivisor;
2248	// Copy over the symbolic part + constant term.
2249	std::copy(first: minLb.begin(), last: minLb.end(), result: newLb.begin() + getNumDimVars());
2250	std::transform(first: newLb.begin() + getNumDimVars(), last: newLb.end(),
2251	result: newLb.begin() + getNumDimVars(),
2252	unary_op: std::negate<DynamicAPInt>());
2253	std::copy(first: maxUb.begin(), last: maxUb.end(), result: newUb.begin() + getNumDimVars());
2254
2255	boundingLbs.emplace_back(args&: newLb);
2256	boundingUbs.emplace_back(args&: newUb);
2257	}
2258
2259	// Clear all constraints and add the lower/upper bounds for the bounding box.
2260	clearConstraints();
2261	for (unsigned d = 0, e = getNumDimVars(); d < e; ++d) {
2262	addInequality(inEq: boundingLbs[d]);
2263	addInequality(inEq: boundingUbs[d]);
2264	}
2265
2266	// Add the constraints that were common to both systems.
2267	append(other: commonCst);
2268	removeTrivialRedundancy();
2269
2270	// TODO: copy over pure symbolic constraints from this and 'other' over to the
2271	// union (since the above are just the union along dimensions); we shouldn't
2272	// be discarding any other constraints on the symbols.
2273
2274	return success();
2275	}
2276
2277	bool IntegerRelation::isColZero(unsigned pos) const {
2278	return !findConstraintWithNonZeroAt(colIdx: pos, /*isEq=*/false) &&
2279	!findConstraintWithNonZeroAt(colIdx: pos, /*isEq=*/true);
2280	}
2281
2282	/// Find positions of inequalities and equalities that do not have a coefficient
2283	/// for [pos, pos + num) variables.
2284	static void getIndependentConstraints(const IntegerRelation &cst, unsigned pos,
2285	unsigned num,
2286	SmallVectorImpl<unsigned> &nbIneqIndices,
2287	SmallVectorImpl<unsigned> &nbEqIndices) {
2288	assert(pos < cst.getNumVars() && "invalid start position");
2289	assert(pos + num <= cst.getNumVars() && "invalid limit");
2290
2291	for (unsigned r = 0, e = cst.getNumInequalities(); r < e; r++) {
2292	// The bounds are to be independent of [offset, offset + num) columns.
2293	unsigned c;
2294	for (c = pos; c < pos + num; ++c) {
2295	if (cst.atIneq(i: r, j: c) != 0)
2296	break;
2297	}
2298	if (c == pos + num)
2299	nbIneqIndices.emplace_back(Args&: r);
2300	}
2301
2302	for (unsigned r = 0, e = cst.getNumEqualities(); r < e; r++) {
2303	// The bounds are to be independent of [offset, offset + num) columns.
2304	unsigned c;
2305	for (c = pos; c < pos + num; ++c) {
2306	if (cst.atEq(i: r, j: c) != 0)
2307	break;
2308	}
2309	if (c == pos + num)
2310	nbEqIndices.emplace_back(Args&: r);
2311	}
2312	}
2313
2314	void IntegerRelation::removeIndependentConstraints(unsigned pos, unsigned num) {
2315	assert(pos + num <= getNumVars() && "invalid range");
2316
2317	// Remove constraints that are independent of these variables.
2318	SmallVector<unsigned, 4> nbIneqIndices, nbEqIndices;
2319	getIndependentConstraints(cst: *this, /*pos=*/0, num, nbIneqIndices, nbEqIndices);
2320
2321	// Iterate in reverse so that indices don't have to be updated.
2322	// TODO: This method can be made more efficient (because removal of each
2323	// inequality leads to much shifting/copying in the underlying buffer).
2324	for (auto nbIndex : llvm::reverse(C&: nbIneqIndices))
2325	removeInequality(pos: nbIndex);
2326	for (auto nbIndex : llvm::reverse(C&: nbEqIndices))
2327	removeEquality(pos: nbIndex);
2328	}
2329
2330	IntegerPolyhedron IntegerRelation::getDomainSet() const {
2331	IntegerRelation copyRel = *this;
2332
2333	// Convert Range variables to Local variables.
2334	copyRel.convertVarKind(srcKind: VarKind::Range, varStart: 0, varLimit: getNumVarKind(kind: VarKind::Range),
2335	dstKind: VarKind::Local);
2336
2337	// Convert Domain variables to SetDim(Range) variables.
2338	copyRel.convertVarKind(srcKind: VarKind::Domain, varStart: 0, varLimit: getNumVarKind(kind: VarKind::Domain),
2339	dstKind: VarKind::SetDim);
2340
2341	return IntegerPolyhedron(std::move(copyRel));
2342	}
2343
2344	bool IntegerRelation::removeDuplicateConstraints() {
2345	bool changed = false;
2346	SmallDenseMap<ArrayRef<DynamicAPInt>, unsigned> hashTable;
2347	unsigned ineqs = getNumInequalities(), cols = getNumCols();
2348
2349	if (ineqs <= 1)
2350	return changed;
2351
2352	// Check if the non-constant part of the constraint is the same.
2353	ArrayRef<DynamicAPInt> row = getInequality(idx: 0).drop_back();
2354	hashTable.insert(KV: {row, 0});
2355	for (unsigned k = 1; k < ineqs; ++k) {
2356	row = getInequality(idx: k).drop_back();
2357	if (hashTable.try_emplace(Key: row, Args&: k).second)
2358	continue;
2359
2360	// For identical cases, keep only the smaller part of the constant term.
2361	unsigned l = hashTable[row];
2362	changed = true;
2363	if (atIneq(i: k, j: cols - 1) <= atIneq(i: l, j: cols - 1))
2364	inequalities.swapRows(row: k, otherRow: l);
2365	removeInequality(pos: k);
2366	--k;
2367	--ineqs;
2368	}
2369
2370	// Check the neg form of each inequality, need an extra vector to store it.
2371	SmallVector<DynamicAPInt> negIneq(cols - 1);
2372	for (unsigned k = 0; k < ineqs; ++k) {
2373	row = getInequality(idx: k).drop_back();
2374	negIneq.assign(in_start: row.begin(), in_end: row.end());
2375	for (DynamicAPInt &ele : negIneq)
2376	ele = -ele;
2377	if (!hashTable.contains(Val: negIneq))
2378	continue;
2379
2380	// For cases where the neg is the same as other inequalities, check that the
2381	// sum of their constant terms is positive.
2382	unsigned l = hashTable[row];
2383	auto sum = atIneq(i: l, j: cols - 1) + atIneq(i: k, j: cols - 1);
2384	if (sum > 0 || l == k)
2385	continue;
2386
2387	// A sum of constant terms equal to zero combines two inequalities into one
2388	// equation, less than zero means the set is empty.
2389	changed = true;
2390	if (k < l)
2391	std::swap(a&: l, b&: k);
2392	if (sum == 0) {
2393	addEquality(eq: getInequality(idx: k));
2394	removeInequality(pos: k);
2395	removeInequality(pos: l);
2396	} else {
2397	*this = getEmpty(space: getSpace());
2398	}
2399	break;
2400	}
2401
2402	return changed;
2403	}
2404
2405	IntegerPolyhedron IntegerRelation::getRangeSet() const {
2406	IntegerRelation copyRel = *this;
2407
2408	// Convert Domain variables to Local variables.
2409	copyRel.convertVarKind(srcKind: VarKind::Domain, varStart: 0, varLimit: getNumVarKind(kind: VarKind::Domain),
2410	dstKind: VarKind::Local);
2411
2412	// We do not need to do anything to Range variables since they are already in
2413	// SetDim position.
2414
2415	return IntegerPolyhedron(std::move(copyRel));
2416	}
2417
2418	void IntegerRelation::intersectDomain(const IntegerPolyhedron &poly) {
2419	assert(getDomainSet().getSpace().isCompatible(poly.getSpace()) &&
2420	"Domain set is not compatible with poly");
2421
2422	// Treating the poly as a relation, convert it from `0 -> R` to `R -> 0`.
2423	IntegerRelation rel = poly;
2424	rel.inverse();
2425
2426	// Append dummy range variables to make the spaces compatible.
2427	rel.appendVar(kind: VarKind::Range, num: getNumRangeVars());
2428
2429	// Intersect in place.
2430	mergeLocalVars(other&: rel);
2431	append(other: rel);
2432	}
2433
2434	void IntegerRelation::intersectRange(const IntegerPolyhedron &poly) {
2435	assert(getRangeSet().getSpace().isCompatible(poly.getSpace()) &&
2436	"Range set is not compatible with poly");
2437
2438	IntegerRelation rel = poly;
2439
2440	// Append dummy domain variables to make the spaces compatible.
2441	rel.appendVar(kind: VarKind::Domain, num: getNumDomainVars());
2442
2443	mergeLocalVars(other&: rel);
2444	append(other: rel);
2445	}
2446
2447	void IntegerRelation::inverse() {
2448	unsigned numRangeVars = getNumVarKind(kind: VarKind::Range);
2449	convertVarKind(srcKind: VarKind::Domain, varStart: 0, varLimit: getVarKindEnd(kind: VarKind::Domain),
2450	dstKind: VarKind::Range);
2451	convertVarKind(srcKind: VarKind::Range, varStart: 0, varLimit: numRangeVars, dstKind: VarKind::Domain);
2452	}
2453
2454	void IntegerRelation::compose(const IntegerRelation &rel) {
2455	assert(getRangeSet().getSpace().isCompatible(rel.getDomainSet().getSpace()) &&
2456	"Range of `this` should be compatible with Domain of `rel`");
2457
2458	IntegerRelation copyRel = rel;
2459
2460	// Let relation `this` be R1: A -> B, and `rel` be R2: B -> C.
2461	// We convert R1 to A -> (B X C), and R2 to B X C then intersect the range of
2462	// R1 with R2. After this, we get R1: A -> C, by projecting out B.
2463	// TODO: Using nested spaces here would help, since we could directly
2464	// intersect the range with another relation.
2465	unsigned numBVars = getNumRangeVars();
2466
2467	// Convert R1 from A -> B to A -> (B X C).
2468	appendVar(kind: VarKind::Range, num: copyRel.getNumRangeVars());
2469
2470	// Convert R2 to B X C.
2471	copyRel.convertVarKind(srcKind: VarKind::Domain, varStart: 0, varLimit: numBVars, dstKind: VarKind::Range, pos: 0);
2472
2473	// Intersect R2 to range of R1.
2474	intersectRange(poly: IntegerPolyhedron(copyRel));
2475
2476	// Project out B in R1.
2477	convertVarKind(srcKind: VarKind::Range, varStart: 0, varLimit: numBVars, dstKind: VarKind::Local);
2478	}
2479
2480	void IntegerRelation::applyDomain(const IntegerRelation &rel) {
2481	inverse();
2482	compose(rel);
2483	inverse();
2484	}
2485
2486	void IntegerRelation::applyRange(const IntegerRelation &rel) { compose(rel); }
2487
2488	void IntegerRelation::printSpace(raw_ostream &os) const {
2489	space.print(os);
2490	os << getNumConstraints() << " constraints\n";
2491	}
2492
2493	void IntegerRelation::removeTrivialEqualities() {
2494	for (int i = getNumEqualities() - 1; i >= 0; --i)
2495	if (rangeIsZero(range: getEquality(idx: i)))
2496	removeEquality(pos: i);
2497	}
2498
2499	bool IntegerRelation::isFullDim() {
2500	if (getNumVars() == 0)
2501	return true;
2502	if (isEmpty())
2503	return false;
2504
2505	// If there is a non-trivial equality, the space cannot be full-dimensional.
2506	removeTrivialEqualities();
2507	if (getNumEqualities() > 0)
2508	return false;
2509
2510	// The polytope is full-dimensional iff it is not flat along any of the
2511	// inequality directions.
2512	Simplex simplex(*this);
2513	return llvm::none_of(Range: llvm::seq<int>(Size: getNumInequalities()), P: [&](int i) {
2514	return simplex.isFlatAlong(coeffs: getInequality(idx: i));
2515	});
2516	}
2517
2518	void IntegerRelation::mergeAndCompose(const IntegerRelation &other) {
2519	assert(getNumDomainVars() == other.getNumRangeVars() &&
2520	"Domain of this and range of other do not match");
2521	// assert(std::equal(values.begin(), values.begin() +
2522	// other.getNumDomainVars(),
2523	// otherValues.begin() + other.getNumDomainVars()) &&
2524	// "Domain of this and range of other do not match");
2525
2526	IntegerRelation result = other;
2527
2528	const unsigned thisDomain = getNumDomainVars();
2529	const unsigned thisRange = getNumRangeVars();
2530	const unsigned otherDomain = other.getNumDomainVars();
2531	const unsigned otherRange = other.getNumRangeVars();
2532
2533	// Add dimension variables temporarily to merge symbol and local vars.
2534	// Convert `this` from
2535	// [thisDomain] -> [thisRange]
2536	// to
2537	// [otherDomain thisDomain] -> [otherRange thisRange].
2538	// and `result` from
2539	// [otherDomain] -> [otherRange]
2540	// to
2541	// [otherDomain thisDomain] -> [otherRange thisRange]
2542	insertVar(kind: VarKind::Domain, pos: 0, num: otherDomain);
2543	insertVar(kind: VarKind::Range, pos: 0, num: otherRange);
2544	result.insertVar(kind: VarKind::Domain, pos: otherDomain, num: thisDomain);
2545	result.insertVar(kind: VarKind::Range, pos: otherRange, num: thisRange);
2546
2547	// Merge symbol and local variables.
2548	mergeAndAlignSymbols(other&: result);
2549	mergeLocalVars(other&: result);
2550
2551	// Convert `result` from [otherDomain thisDomain] -> [otherRange thisRange] to
2552	// [otherDomain] -> [thisRange]
2553	result.removeVarRange(kind: VarKind::Domain, varStart: otherDomain, varLimit: otherDomain + thisDomain);
2554	result.convertToLocal(kind: VarKind::Range, varStart: 0, varLimit: otherRange);
2555	// Convert `this` from [otherDomain thisDomain] -> [otherRange thisRange] to
2556	// [otherDomain] -> [thisRange]
2557	convertToLocal(kind: VarKind::Domain, varStart: otherDomain, varLimit: otherDomain + thisDomain);
2558	removeVarRange(kind: VarKind::Range, varStart: 0, varLimit: otherRange);
2559
2560	// Add and match domain of `result` to domain of `this`.
2561	for (unsigned i = 0, e = result.getNumDomainVars(); i < e; ++i)
2562	if (result.getSpace().getId(kind: VarKind::Domain, pos: i).hasValue())
2563	space.setId(kind: VarKind::Domain, pos: i,
2564	id: result.getSpace().getId(kind: VarKind::Domain, pos: i));
2565	// Add and match range of `this` to range of `result`.
2566	for (unsigned i = 0, e = getNumRangeVars(); i < e; ++i)
2567	if (space.getId(kind: VarKind::Range, pos: i).hasValue())
2568	result.space.setId(kind: VarKind::Range, pos: i, id: space.getId(kind: VarKind::Range, pos: i));
2569
2570	// Append `this` to `result` and simplify constraints.
2571	result.append(other: *this);
2572	result.removeRedundantLocalVars();
2573
2574	*this = result;
2575	}
2576
2577	void IntegerRelation::print(raw_ostream &os) const {
2578	assert(hasConsistentState());
2579	printSpace(os);
2580	PrintTableMetrics ptm = {.maxPreIndent: 0, .maxPostIndent: 0, .preAlign: "-"};
2581	for (unsigned i = 0, e = getNumEqualities(); i < e; ++i)
2582	for (unsigned j = 0, f = getNumCols(); j < f; ++j)
2583	updatePrintMetrics<DynamicAPInt>(val: atEq(i, j), m&: ptm);
2584	for (unsigned i = 0, e = getNumInequalities(); i < e; ++i)
2585	for (unsigned j = 0, f = getNumCols(); j < f; ++j)
2586	updatePrintMetrics<DynamicAPInt>(val: atIneq(i, j), m&: ptm);
2587	// Print using PrintMetrics.
2588	constexpr unsigned kMinSpacing = 1;
2589	for (unsigned i = 0, e = getNumEqualities(); i < e; ++i) {
2590	for (unsigned j = 0, f = getNumCols(); j < f; ++j) {
2591	printWithPrintMetrics<DynamicAPInt>(os, val: atEq(i, j), minSpacing: kMinSpacing, m: ptm);
2592	}
2593	os << " = 0\n";
2594	}
2595	for (unsigned i = 0, e = getNumInequalities(); i < e; ++i) {
2596	for (unsigned j = 0, f = getNumCols(); j < f; ++j) {
2597	printWithPrintMetrics<DynamicAPInt>(os, val: atIneq(i, j), minSpacing: kMinSpacing, m: ptm);
2598	}
2599	os << " >= 0\n";
2600	}
2601	os << '\n';
2602	}
2603
2604	void IntegerRelation::dump() const { print(os&: llvm::errs()); }
2605
2606	unsigned IntegerPolyhedron::insertVar(VarKind kind, unsigned pos,
2607	unsigned num) {
2608	assert((kind != VarKind::Domain || num == 0) &&
2609	"Domain has to be zero in a set");
2610	return IntegerRelation::insertVar(kind, pos, num);
2611	}
2612	IntegerPolyhedron
2613	IntegerPolyhedron::intersect(const IntegerPolyhedron &other) const {
2614	return IntegerPolyhedron(IntegerRelation::intersect(other));
2615	}
2616
2617	PresburgerSet IntegerPolyhedron::subtract(const PresburgerSet &other) const {
2618	return PresburgerSet(IntegerRelation::subtract(set: other));
2619	}
2620