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