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