brintos / llvm-project-archived public Read only

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