brintos / llvm-project-archived public Read only

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