572 lines · cpp
1//===- Utils.cpp - General utilities for Presburger library ---------------===//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// Utility functions required by the Presburger Library.10//11//===----------------------------------------------------------------------===//12 13#include "mlir/Analysis/Presburger/Utils.h"14#include "mlir/Analysis/Presburger/IntegerRelation.h"15#include "mlir/Analysis/Presburger/PresburgerSpace.h"16#include "llvm/ADT/STLFunctionalExtras.h"17#include "llvm/ADT/SmallBitVector.h"18#include "llvm/Support/raw_ostream.h"19#include <cassert>20#include <cstdint>21#include <numeric>22#include <optional>23 24using namespace mlir;25using namespace presburger;26using llvm::dynamicAPIntFromInt64;27 28/// Normalize a division's `dividend` and the `divisor` by their GCD. For29/// example: if the dividend and divisor are [2,0,4] and 4 respectively,30/// they get normalized to [1,0,2] and 2. The divisor must be non-negative;31/// it is allowed for the divisor to be zero, but nothing is done in this case.32static void normalizeDivisionByGCD(MutableArrayRef<DynamicAPInt> dividend,33 DynamicAPInt &divisor) {34 assert(divisor > 0 && "divisor must be non-negative!");35 if (dividend.empty())36 return;37 // We take the absolute value of dividend's coefficients to make sure that38 // `gcd` is positive.39 DynamicAPInt gcd = llvm::gcd(abs(dividend.front()), divisor);40 41 // The reason for ignoring the constant term is as follows.42 // For a division:43 // floor((a + m.f(x))/(m.d))44 // It can be replaced by:45 // floor((floor(a/m) + f(x))/d)46 // Since `{a/m}/d` in the dividend satisfies 0 <= {a/m}/d < 1/d, it will not47 // influence the result of the floor division and thus, can be ignored.48 for (size_t i = 1, m = dividend.size() - 1; i < m; i++) {49 gcd = llvm::gcd(abs(dividend[i]), gcd);50 if (gcd == 1)51 return;52 }53 54 // Normalize the dividend and the denominator.55 llvm::transform(dividend, dividend.begin(),56 [gcd](DynamicAPInt &n) { return floorDiv(n, gcd); });57 divisor /= gcd;58}59 60/// Check if the pos^th variable can be represented as a division using upper61/// bound inequality at position `ubIneq` and lower bound inequality at position62/// `lbIneq`.63///64/// Let `var` be the pos^th variable, then `var` is equivalent to65/// `expr floordiv divisor` if there are constraints of the form:66/// 0 <= expr - divisor * var <= divisor - 167/// Rearranging, we have:68/// divisor * var - expr + (divisor - 1) >= 0 <-- Lower bound for 'var'69/// -divisor * var + expr >= 0 <-- Upper bound for 'var'70///71/// For example:72/// 32*k >= 16*i + j - 31 <-- Lower bound for 'k'73/// 32*k <= 16*i + j <-- Upper bound for 'k'74/// expr = 16*i + j, divisor = 3275/// k = ( 16*i + j ) floordiv 3276///77/// 4q >= i + j - 2 <-- Lower bound for 'q'78/// 4q <= i + j + 1 <-- Upper bound for 'q'79/// expr = i + j + 1, divisor = 480/// q = (i + j + 1) floordiv 481//82/// This function also supports detecting divisions from bounds that are83/// strictly tighter than the division bounds described above, since tighter84/// bounds imply the division bounds. For example:85/// 4q - i - j + 2 >= 0 <-- Lower bound for 'q'86/// -4q + i + j >= 0 <-- Tight upper bound for 'q'87///88/// To extract floor divisions with tighter bounds, we assume that the89/// constraints are of the form:90/// c <= expr - divisior * var <= divisor - 1, where 0 <= c <= divisor - 191/// Rearranging, we have:92/// divisor * var - expr + (divisor - 1) >= 0 <-- Lower bound for 'var'93/// -divisor * var + expr - c >= 0 <-- Upper bound for 'var'94///95/// If successful, `expr` is set to dividend of the division and `divisor` is96/// set to the denominator of the division, which will be positive.97/// The final division expression is normalized by GCD.98static LogicalResult getDivRepr(const IntegerRelation &cst, unsigned pos,99 unsigned ubIneq, unsigned lbIneq,100 MutableArrayRef<DynamicAPInt> expr,101 DynamicAPInt &divisor) {102 103 assert(pos <= cst.getNumVars() && "Invalid variable position");104 assert(ubIneq <= cst.getNumInequalities() &&105 "Invalid upper bound inequality position");106 assert(lbIneq <= cst.getNumInequalities() &&107 "Invalid upper bound inequality position");108 assert(expr.size() == cst.getNumCols() && "Invalid expression size");109 assert(cst.atIneq(lbIneq, pos) > 0 && "lbIneq is not a lower bound!");110 assert(cst.atIneq(ubIneq, pos) < 0 && "ubIneq is not an upper bound!");111 112 // Extract divisor from the lower bound.113 divisor = cst.atIneq(lbIneq, pos);114 115 // First, check if the constraints are opposite of each other except the116 // constant term.117 unsigned i = 0, e = 0;118 for (i = 0, e = cst.getNumVars(); i < e; ++i)119 if (cst.atIneq(ubIneq, i) != -cst.atIneq(lbIneq, i))120 break;121 122 if (i < e)123 return failure();124 125 // Then, check if the constant term is of the proper form.126 // Due to the form of the upper/lower bound inequalities, the sum of their127 // constants is `divisor - 1 - c`. From this, we can extract c:128 DynamicAPInt constantSum = cst.atIneq(lbIneq, cst.getNumCols() - 1) +129 cst.atIneq(ubIneq, cst.getNumCols() - 1);130 DynamicAPInt c = divisor - 1 - constantSum;131 132 // Check if `c` satisfies the condition `0 <= c <= divisor - 1`.133 // This also implictly checks that `divisor` is positive.134 if (!(0 <= c && c <= divisor - 1)) // NOLINT135 return failure();136 137 // The inequality pair can be used to extract the division.138 // Set `expr` to the dividend of the division except the constant term, which139 // is set below.140 for (i = 0, e = cst.getNumVars(); i < e; ++i)141 if (i != pos)142 expr[i] = cst.atIneq(ubIneq, i);143 144 // From the upper bound inequality's form, its constant term is equal to the145 // constant term of `expr`, minus `c`. From this,146 // constant term of `expr` = constant term of upper bound + `c`.147 expr.back() = cst.atIneq(ubIneq, cst.getNumCols() - 1) + c;148 normalizeDivisionByGCD(expr, divisor);149 150 return success();151}152 153/// Check if the pos^th variable can be represented as a division using154/// equality at position `eqInd`.155///156/// For example:157/// 32*k == 16*i + j - 31 <-- `eqInd` for 'k'158/// expr = 16*i + j - 31, divisor = 32159/// k = (16*i + j - 31) floordiv 32160///161/// If successful, `expr` is set to dividend of the division and `divisor` is162/// set to the denominator of the division. The final division expression is163/// normalized by GCD.164static LogicalResult getDivRepr(const IntegerRelation &cst, unsigned pos,165 unsigned eqInd,166 MutableArrayRef<DynamicAPInt> expr,167 DynamicAPInt &divisor) {168 169 assert(pos <= cst.getNumVars() && "Invalid variable position");170 assert(eqInd <= cst.getNumEqualities() && "Invalid equality position");171 assert(expr.size() == cst.getNumCols() && "Invalid expression size");172 173 // Extract divisor, the divisor can be negative and hence its sign information174 // is stored in `signDiv` to reverse the sign of dividend's coefficients.175 // Equality must involve the pos-th variable and hence `tempDiv` != 0.176 DynamicAPInt tempDiv = cst.atEq(eqInd, pos);177 if (tempDiv == 0)178 return failure();179 int signDiv = tempDiv < 0 ? -1 : 1;180 181 // The divisor is always a positive integer.182 divisor = tempDiv * signDiv;183 184 for (unsigned i = 0, e = cst.getNumVars(); i < e; ++i)185 if (i != pos)186 expr[i] = -signDiv * cst.atEq(eqInd, i);187 188 expr.back() = -signDiv * cst.atEq(eqInd, cst.getNumCols() - 1);189 normalizeDivisionByGCD(expr, divisor);190 191 return success();192}193 194// Returns `false` if the constraints depends on a variable for which an195// explicit representation has not been found yet, otherwise returns `true`.196static bool checkExplicitRepresentation(const IntegerRelation &cst,197 ArrayRef<bool> foundRepr,198 ArrayRef<DynamicAPInt> dividend,199 unsigned pos) {200 // Exit to avoid circular dependencies between divisions.201 for (unsigned c = 0, e = cst.getNumVars(); c < e; ++c) {202 if (c == pos)203 continue;204 205 if (!foundRepr[c] && dividend[c] != 0) {206 // Expression can't be constructed as it depends on a yet unknown207 // variable.208 //209 // TODO: Visit/compute the variables in an order so that this doesn't210 // happen. More complex but much more efficient.211 return false;212 }213 }214 215 return true;216}217 218/// Check if the pos^th variable can be expressed as a floordiv of an affine219/// function of other variables (where the divisor is a positive constant).220/// `foundRepr` contains a boolean for each variable indicating if the221/// explicit representation for that variable has already been computed.222/// Returns the `MaybeLocalRepr` struct which contains the indices of the223/// constraints that can be expressed as a floordiv of an affine function. If224/// the representation could be computed, `dividend` and `denominator` are set.225/// If the representation could not be computed, the kind attribute in226/// `MaybeLocalRepr` is set to None.227MaybeLocalRepr presburger::computeSingleVarRepr(228 const IntegerRelation &cst, ArrayRef<bool> foundRepr, unsigned pos,229 MutableArrayRef<DynamicAPInt> dividend, DynamicAPInt &divisor) {230 assert(pos < cst.getNumVars() && "invalid position");231 assert(foundRepr.size() == cst.getNumVars() &&232 "Size of foundRepr does not match total number of variables");233 assert(dividend.size() == cst.getNumCols() && "Invalid dividend size");234 235 SmallVector<unsigned, 4> lbIndices, ubIndices, eqIndices;236 cst.getLowerAndUpperBoundIndices(pos, &lbIndices, &ubIndices, &eqIndices);237 MaybeLocalRepr repr{};238 239 for (unsigned ubPos : ubIndices) {240 for (unsigned lbPos : lbIndices) {241 // Attempt to get divison representation from ubPos, lbPos.242 if (getDivRepr(cst, pos, ubPos, lbPos, dividend, divisor).failed())243 continue;244 245 if (!checkExplicitRepresentation(cst, foundRepr, dividend, pos))246 continue;247 248 repr.kind = ReprKind::Inequality;249 repr.repr.inequalityPair = {ubPos, lbPos};250 return repr;251 }252 }253 for (unsigned eqPos : eqIndices) {254 // Attempt to get divison representation from eqPos.255 if (getDivRepr(cst, pos, eqPos, dividend, divisor).failed())256 continue;257 258 if (!checkExplicitRepresentation(cst, foundRepr, dividend, pos))259 continue;260 261 repr.kind = ReprKind::Equality;262 repr.repr.equalityIdx = eqPos;263 return repr;264 }265 return repr;266}267 268MaybeLocalRepr presburger::computeSingleVarRepr(269 const IntegerRelation &cst, ArrayRef<bool> foundRepr, unsigned pos,270 SmallVector<int64_t, 8> ÷nd, unsigned &divisor) {271 SmallVector<DynamicAPInt, 8> dividendDynamicAPInt(cst.getNumCols());272 DynamicAPInt divisorDynamicAPInt;273 MaybeLocalRepr result = computeSingleVarRepr(274 cst, foundRepr, pos, dividendDynamicAPInt, divisorDynamicAPInt);275 dividend = getInt64Vec(dividendDynamicAPInt);276 divisor = unsigned(int64_t(divisorDynamicAPInt));277 return result;278}279 280llvm::SmallBitVector presburger::getSubrangeBitVector(unsigned len,281 unsigned setOffset,282 unsigned numSet) {283 llvm::SmallBitVector vec(len, false);284 vec.set(setOffset, setOffset + numSet);285 return vec;286}287 288void presburger::mergeLocalVars(289 IntegerRelation &relA, IntegerRelation &relB,290 llvm::function_ref<bool(unsigned i, unsigned j)> merge) {291 assert(relA.getSpace().isCompatible(relB.getSpace()) &&292 "Spaces should be compatible.");293 294 // Merge local vars of relA and relB without using division information,295 // i.e. append local vars of `relB` to `relA` and insert local vars of `relA`296 // to `relB` at start of its local vars.297 unsigned initLocals = relA.getNumLocalVars();298 relA.insertVar(VarKind::Local, relA.getNumLocalVars(),299 relB.getNumLocalVars());300 relB.insertVar(VarKind::Local, 0, initLocals);301 302 // Get division representations from each rel.303 DivisionRepr divsA = relA.getLocalReprs();304 DivisionRepr divsB = relB.getLocalReprs();305 306 for (unsigned i = initLocals, e = divsB.getNumDivs(); i < e; ++i)307 divsA.setDiv(i, divsB.getDividend(i), divsB.getDenom(i));308 309 // Remove duplicate divisions from divsA. The removing duplicate divisions310 // call, calls `merge` to effectively merge divisions in relA and relB.311 divsA.removeDuplicateDivs(merge);312}313 314SmallVector<DynamicAPInt, 8>315presburger::getDivUpperBound(ArrayRef<DynamicAPInt> dividend,316 const DynamicAPInt &divisor,317 unsigned localVarIdx) {318 assert(divisor > 0 && "divisor must be positive!");319 assert(dividend[localVarIdx] == 0 &&320 "Local to be set to division must have zero coeff!");321 SmallVector<DynamicAPInt, 8> ineq(dividend);322 ineq[localVarIdx] = -divisor;323 return ineq;324}325 326SmallVector<DynamicAPInt, 8>327presburger::getDivLowerBound(ArrayRef<DynamicAPInt> dividend,328 const DynamicAPInt &divisor,329 unsigned localVarIdx) {330 assert(divisor > 0 && "divisor must be positive!");331 assert(dividend[localVarIdx] == 0 &&332 "Local to be set to division must have zero coeff!");333 SmallVector<DynamicAPInt, 8> ineq(dividend.size());334 llvm::transform(dividend, ineq.begin(), std::negate<DynamicAPInt>());335 ineq[localVarIdx] = divisor;336 ineq.back() += divisor - 1;337 return ineq;338}339 340DynamicAPInt presburger::gcdRange(ArrayRef<DynamicAPInt> range) {341 DynamicAPInt gcd(0);342 for (const DynamicAPInt &elem : range) {343 gcd = llvm::gcd(gcd, abs(elem));344 if (gcd == 1)345 return gcd;346 }347 return gcd;348}349 350DynamicAPInt presburger::normalizeRange(MutableArrayRef<DynamicAPInt> range) {351 DynamicAPInt gcd = gcdRange(range);352 if ((gcd == 0) || (gcd == 1))353 return gcd;354 for (DynamicAPInt &elem : range)355 elem /= gcd;356 return gcd;357}358 359void presburger::normalizeDiv(MutableArrayRef<DynamicAPInt> num,360 DynamicAPInt &denom) {361 assert(denom > 0 && "denom must be positive!");362 DynamicAPInt gcd = llvm::gcd(gcdRange(num), denom);363 if (gcd == 1)364 return;365 for (DynamicAPInt &coeff : num)366 coeff /= gcd;367 denom /= gcd;368}369 370SmallVector<DynamicAPInt, 8>371presburger::getNegatedCoeffs(ArrayRef<DynamicAPInt> coeffs) {372 SmallVector<DynamicAPInt, 8> negatedCoeffs;373 negatedCoeffs.reserve(coeffs.size());374 for (const DynamicAPInt &coeff : coeffs)375 negatedCoeffs.emplace_back(-coeff);376 return negatedCoeffs;377}378 379SmallVector<DynamicAPInt, 8>380presburger::getComplementIneq(ArrayRef<DynamicAPInt> ineq) {381 SmallVector<DynamicAPInt, 8> coeffs;382 coeffs.reserve(ineq.size());383 for (const DynamicAPInt &coeff : ineq)384 coeffs.emplace_back(-coeff);385 --coeffs.back();386 return coeffs;387}388 389SmallVector<std::optional<DynamicAPInt>, 4>390DivisionRepr::divValuesAt(ArrayRef<DynamicAPInt> point) const {391 assert(point.size() == getNumNonDivs() && "Incorrect point size");392 393 SmallVector<std::optional<DynamicAPInt>, 4> divValues(getNumDivs(),394 std::nullopt);395 bool changed = true;396 while (changed) {397 changed = false;398 for (unsigned i = 0, e = getNumDivs(); i < e; ++i) {399 // If division value is found, continue;400 if (divValues[i])401 continue;402 403 ArrayRef<DynamicAPInt> dividend = getDividend(i);404 DynamicAPInt divVal(0);405 406 // Check if we have all the division values required for this division.407 unsigned j, f;408 for (j = 0, f = getNumDivs(); j < f; ++j) {409 if (dividend[getDivOffset() + j] == 0)410 continue;411 // Division value required, but not found yet.412 if (!divValues[j])413 break;414 divVal += dividend[getDivOffset() + j] * *divValues[j];415 }416 417 // We have some division values that are still not found, but are required418 // to find the value of this division.419 if (j < f)420 continue;421 422 // Fill remaining values.423 divVal = std::inner_product(point.begin(), point.end(), dividend.begin(),424 divVal);425 // Add constant.426 divVal += dividend.back();427 // Take floor division with denominator.428 divVal = floorDiv(divVal, denoms[i]);429 430 // Set div value and continue.431 divValues[i] = divVal;432 changed = true;433 }434 }435 436 return divValues;437}438 439void DivisionRepr::removeDuplicateDivs(440 llvm::function_ref<bool(unsigned i, unsigned j)> merge) {441 442 // Find and merge duplicate divisions.443 // TODO: Add division normalization to support divisions that differ by444 // a constant.445 // TODO: Add division ordering such that a division representation for local446 // variable at position `i` only depends on local variables at position <447 // `i`. This would make sure that all divisions depending on other local448 // variables that can be merged, are merged.449 normalizeDivs();450 for (unsigned i = 0; i < getNumDivs(); ++i) {451 // Check if a division representation exists for the `i^th` local var.452 if (denoms[i] == 0)453 continue;454 // Check if a division exists which is a duplicate of the division at `i`.455 for (unsigned j = i + 1; j < getNumDivs(); ++j) {456 // Check if a division representation exists for the `j^th` local var.457 if (denoms[j] == 0)458 continue;459 // Check if the denominators match.460 if (denoms[i] != denoms[j])461 continue;462 // Check if the representations are equal.463 if (dividends.getRow(i) != dividends.getRow(j))464 continue;465 466 // Merge divisions at position `j` into division at position `i`. If467 // merge fails, do not merge these divs.468 bool mergeResult = merge(i, j);469 if (!mergeResult)470 continue;471 472 // Update division information to reflect merging.473 unsigned divOffset = getDivOffset();474 dividends.addToColumn(divOffset + j, divOffset + i, /*scale=*/1);475 dividends.removeColumn(divOffset + j);476 dividends.removeRow(j);477 denoms.erase(denoms.begin() + j);478 479 // Since `j` can never be zero, we do not need to worry about overflows.480 --j;481 }482 }483}484 485void DivisionRepr::normalizeDivs() {486 for (unsigned i = 0, e = getNumDivs(); i < e; ++i) {487 if (getDenom(i) == 0 || getDividend(i).empty())488 continue;489 normalizeDiv(getDividend(i), getDenom(i));490 }491}492 493void DivisionRepr::insertDiv(unsigned pos, ArrayRef<DynamicAPInt> dividend,494 const DynamicAPInt &divisor) {495 assert(pos <= getNumDivs() && "Invalid insertion position");496 assert(dividend.size() == getNumVars() + 1 && "Incorrect dividend size");497 498 dividends.appendExtraRow(dividend);499 denoms.insert(denoms.begin() + pos, divisor);500 dividends.insertColumn(getDivOffset() + pos);501}502 503void DivisionRepr::insertDiv(unsigned pos, unsigned num) {504 assert(pos <= getNumDivs() && "Invalid insertion position");505 dividends.insertColumns(getDivOffset() + pos, num);506 dividends.insertRows(pos, num);507 denoms.insert(denoms.begin() + pos, num, DynamicAPInt(0));508}509 510void DivisionRepr::print(raw_ostream &os) const {511 os << "Dividends:\n";512 dividends.print(os);513 os << "Denominators\n";514 for (const DynamicAPInt &denom : denoms)515 os << denom << " ";516 os << "\n";517}518 519void DivisionRepr::dump() const { print(llvm::errs()); }520 521SmallVector<DynamicAPInt, 8>522presburger::getDynamicAPIntVec(ArrayRef<int64_t> range) {523 SmallVector<DynamicAPInt, 8> result(range.size());524 llvm::transform(range, result.begin(), dynamicAPIntFromInt64);525 return result;526}527 528SmallVector<int64_t, 8> presburger::getInt64Vec(ArrayRef<DynamicAPInt> range) {529 SmallVector<int64_t, 8> result(range.size());530 llvm::transform(range, result.begin(), int64fromDynamicAPInt);531 return result;532}533 534Fraction presburger::dotProduct(ArrayRef<Fraction> a, ArrayRef<Fraction> b) {535 assert(a.size() == b.size() &&536 "dot product is only valid for vectors of equal sizes!");537 Fraction sum = 0;538 for (unsigned i = 0, e = a.size(); i < e; i++)539 sum += a[i] * b[i];540 return sum;541}542 543/// Find the product of two polynomials, each given by an array of544/// coefficients, by taking the convolution.545std::vector<Fraction> presburger::multiplyPolynomials(ArrayRef<Fraction> a,546 ArrayRef<Fraction> b) {547 // The length of the convolution is the sum of the lengths548 // of the two sequences. We pad the shorter one with zeroes.549 unsigned len = a.size() + b.size() - 1;550 551 // We define accessors to avoid out-of-bounds errors.552 auto getCoeff = [](ArrayRef<Fraction> arr, unsigned i) -> Fraction {553 if (i < arr.size())554 return arr[i];555 return 0;556 };557 558 std::vector<Fraction> convolution;559 convolution.reserve(len);560 for (unsigned k = 0; k < len; ++k) {561 Fraction sum(0, 1);562 for (unsigned l = 0; l <= k; ++l)563 sum += getCoeff(a, l) * getCoeff(b, k - l);564 convolution.emplace_back(sum);565 }566 return convolution;567}568 569bool presburger::isRangeZero(ArrayRef<Fraction> arr) {570 return llvm::all_of(arr, [](const Fraction &f) { return f == 0; });571}572