diff --git a/libsolutil/BooleanLP.cpp b/libsolutil/BooleanLP.cpp index 457a4309f..09eb5d183 100644 --- a/libsolutil/BooleanLP.cpp +++ b/libsolutil/BooleanLP.cpp @@ -24,7 +24,12 @@ #include #include + +#if LPIncremental +#include +#else #include +#endif #include #include @@ -48,6 +53,7 @@ using rational = boost::rational; //#define DEBUG + namespace { template @@ -120,12 +126,16 @@ pair> BooleanLPSolver::check(vector cons std::vector booleanVariables; std::vector clauses = state().clauses; +#if LPIncremental + LPSolver lpSolver; +#else // TODO we start building up a new set of solver // for each query, but we should also keep some // kind of cache across queries. std::vector> lpSolvers; lpSolvers.emplace_back(0, LPSolver{}); LPSolver& lpSolver = lpSolvers.back().second; +#endif for (auto&& [index, bound]: state().bounds) { @@ -136,6 +146,10 @@ pair> BooleanLPSolver::check(vector cons } for (Constraint const& c: state().fixedConstraints) lpSolver.addConstraint(c); +#if LPIncremental + for (auto&& [index, constraint]: state().conditionalConstraints) + lpSolver.addConditionalConstraint(constraint, index); +#endif // TODO this way, it will result in a lot of gaps in both sets of variables. // should we compress them and store a mapping? @@ -146,6 +160,9 @@ pair> BooleanLPSolver::check(vector cons else lpSolver.setVariableName(index, name); +#ifdef DEBUG + cerr << "Performing preliminary check." << endl; +#endif if (lpSolver.check().first == LPResult::Infeasible) { #ifdef DEBUG @@ -156,17 +173,32 @@ pair> BooleanLPSolver::check(vector cons auto theorySolver = [&](size_t _trailSize, map const& _newBooleanAssignment) -> optional { +#if LPIncremental + lpSolver.setTrailSize(_trailSize); +#else lpSolvers.emplace_back(_trailSize, LPSolver(lpSolvers.back().second)); +#endif for (auto&& [constraintIndex, value]: _newBooleanAssignment) { if (!value || !state().conditionalConstraints.count(constraintIndex)) continue; +#if LPIncremental + lpSolver.activateConstraint(constraintIndex); +#else // "reason" is already stored for those constraints. Constraint const& constraint = state().conditionalConstraints.at(constraintIndex); lpSolvers.back().second.addConstraint(constraint, constraintIndex); +#endif } +#ifdef DEBUG + cerr << "Performing incremental check." << endl; +#endif +#if LPIncremental + auto&& [result, reasonSet] = lpSolver.check(); +#else auto&& [result, reasonSet] = lpSolvers.back().second.check(); +#endif // We can only really use the result "infeasible". Everything else should be "sat". if (result == LPResult::Infeasible) { @@ -184,8 +216,12 @@ pair> BooleanLPSolver::check(vector cons }; auto backtrackNotify = [&](size_t _trailSize) { +#if LPIncremental + lpSolver.setTrailSize(_trailSize); +#else while (lpSolvers.back().first > _trailSize) lpSolvers.pop_back(); +#endif }; auto optionalModel = CDCL{move(booleanVariables), clauses, theorySolver, backtrackNotify}.solve(); diff --git a/libsolutil/BooleanLP.h b/libsolutil/BooleanLP.h index 1e7d8fc17..9e0c43bc4 100644 --- a/libsolutil/BooleanLP.h +++ b/libsolutil/BooleanLP.h @@ -19,7 +19,13 @@ #include +#define LPIncremental 1 + +#if LPIncremental +#include +#else #include +#endif #include #include diff --git a/libsolutil/CMakeLists.txt b/libsolutil/CMakeLists.txt index 32aaf5264..4d60480fe 100644 --- a/libsolutil/CMakeLists.txt +++ b/libsolutil/CMakeLists.txt @@ -29,6 +29,8 @@ set(sources LEB128.h LP.cpp LP.h + LPIncremental.cpp + LPIncremental.h Numeric.cpp Numeric.h picosha2.h diff --git a/libsolutil/LP.cpp b/libsolutil/LP.cpp index 56d5aa24c..69ee570d6 100644 --- a/libsolutil/LP.cpp +++ b/libsolutil/LP.cpp @@ -18,6 +18,8 @@ #include +#ifndef LPIncremental + #include #include #include @@ -51,7 +53,7 @@ using namespace solidity::util; using rational = boost::rational; -#define DEBUG +//#define DEBUG namespace { @@ -833,3 +835,5 @@ void LPSolver::SubProblem::pivotAndUpdate( pivot(_oldBasicVar, _newBasicVar); } + +#endif diff --git a/libsolutil/LP.h b/libsolutil/LP.h index bed034243..86271a686 100644 --- a/libsolutil/LP.h +++ b/libsolutil/LP.h @@ -19,7 +19,10 @@ // use sparse matrices #define SPARSE 1 -#define DEBUG +//#define DEBUG +#define LPIncremental 1 + +#ifndef LPIncremental #include #include @@ -257,3 +260,4 @@ private: }; } +#endif diff --git a/libsolutil/LPIncremental.cpp b/libsolutil/LPIncremental.cpp new file mode 100644 index 000000000..6604ba5d2 --- /dev/null +++ b/libsolutil/LPIncremental.cpp @@ -0,0 +1,741 @@ +/* + This file is part of solidity. + + solidity is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + solidity is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with solidity. If not, see . +*/ +// SPDX-License-Identifier: GPL-3.0 + +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +using namespace std; +using namespace solidity; +using namespace solidity::util; + +using rational = boost::rational; + +//#define DEBUG + +namespace +{ + +/// Disjunctively combined two vectors of bools. +inline std::vector& operator|=(std::vector& _x, std::vector const& _y) +{ + solAssert(_x.size() == _y.size(), ""); + for (size_t i = 0; i < _x.size(); ++i) + if (_y[i]) + _x[i] = true; + return _x; +} + +string toString(rational const& _x) +{ + if (_x == bigint(1) << 256) + return "2**256"; + else if (_x == (bigint(1) << 256) - 1) + return "2**256-1"; + else if (_x.denominator() == 1) + return ::toString(_x.numerator()); + else + return ::toString(_x.numerator()) + "/" + ::toString(_x.denominator()); +} + +/* +string reasonToString(ReasonSet const& _reasons, size_t _minSize) +{ + auto reasonsAsStrings = _reasons | ranges::views::transform([](size_t _r) { return to_string(_r); }); + string result = "[" + joinHumanReadable(reasonsAsStrings) + "]"; + if (result.size() < _minSize) + result.resize(_minSize, ' '); + return result; +} +*/ + +} + +bool Constraint::operator<(Constraint const& _other) const +{ + if (kind != _other.kind) + return kind < _other.kind; + + for (size_t i = 0; i < max(data.size(), _other.data.size()); ++i) + if (rational diff = data.get(i) - _other.data.get(i)) + { + //cerr << "Exit after " << i << endl; + return diff < 0; + } + //cerr << "full traversal of " << max(data.size(), _other.data.size()) << endl; + + return false; +} + +bool Constraint::operator==(Constraint const& _other) const +{ + if (kind != _other.kind) + return false; + + for (size_t i = 0; i < max(data.size(), _other.data.size()); ++i) + if (data.get(i) != _other.data.get(i)) + { + //cerr << "Exit after " << i << endl; + return false; + } + //cerr << "full traversal of " << max(data.size(), _other.data.size()) << endl; + + return true; +} + +string RationalWithDelta::toString() const +{ + string result = ::toString(m_main); + if (m_delta) + result += + (m_delta > 0 ? "+" : "-") + + (abs(m_delta) == 1 ? "" : ::toString(abs(m_delta))) + + "d"; + return result; +} + +void LPSolver::addConstraint(Constraint const& _constraint) +{ +#ifdef DEBUG + cerr << "Adding constraint." << endl; +#endif + result = nullopt; + auto&& [varIndex, bounds] = constraintIntoVariableBounds(_constraint); + addBounds(varIndex, bounds); +} + +void LPSolver::addLowerBound(size_t _variable, RationalWithDelta _bound) +{ +#ifdef DEBUG + cerr << "Adding lower bound." << endl; +#endif + result = nullopt; + size_t innerIndex = maybeAddOuterVariable(_variable); + addBounds(innerIndex, Bounds{move(_bound), {}}); +} + +void LPSolver::addUpperBound(size_t _variable, RationalWithDelta _bound) +{ +#ifdef DEBUG + cerr << "Adding upper bound." << endl; +#endif + // TODO we could only reset the result if the bound changed anything. + // then we could check if we already have a result insiche "check()" + // and return early. Although this might be better done inside + // activateConstraint. + result = nullopt; + size_t innerIndex = maybeAddOuterVariable(_variable); + addBounds(innerIndex, Bounds{{}, move(_bound)}); +} + +void LPSolver::addConditionalConstraint(Constraint const& _constraint, size_t _reason) +{ +#ifdef DEBUG + cerr << "Adding conditional constraint." << endl; +#endif + auto&& [varIndex, bounds] = constraintIntoVariableBounds(_constraint); + solAssert(!reasonToBounds.count(_reason)); + reasonToBounds[_reason] = make_pair(varIndex, move(bounds)); +} + +void LPSolver::activateConstraint(size_t _reason) +{ +#ifdef DEBUG + cerr << "Activating constraint." << endl; +#endif + result = nullopt; + auto&& [varIndex, bounds] = reasonToBounds.at(_reason); + Variable& var = variables[varIndex]; + bool savedBounds = false; + if (bounds.lower && (!var.bounds.lower || *var.bounds.lower < *bounds.lower)) + { + storedBounds.emplace_back(make_tuple(trailSize, varIndex, var.bounds, var.lowerReason, var.upperReason)); + savedBounds = true; + var.bounds.lower = bounds.lower; + var.lowerReason = _reason; + if (var.value < *var.bounds.lower) + variablesPotentiallyOutOfBounds.insert(varIndex); + } + if (bounds.upper && (!var.bounds.upper || *var.bounds.upper > *bounds.upper)) + { + if (!savedBounds) + storedBounds.emplace_back(make_tuple(trailSize, varIndex, var.bounds, var.lowerReason, var.upperReason)); + savedBounds = true; + var.bounds.upper = bounds.upper; + var.upperReason = _reason; + if (var.value > *var.bounds.upper) + variablesPotentiallyOutOfBounds.insert(varIndex); + } +#ifdef DEBUG + if (!savedBounds) + cerr << "Did not change anything." << endl; +#endif +} + +void LPSolver::setTrailSize(size_t _trailSize) +{ +// solAssert(_trailSize == 0 || _trailSize != trailSize); + if (_trailSize > trailSize) + { +#ifdef DEBUG + cerr << "=== Advancing from " << trailSize << " to " << _trailSize << endl; +#endif + solAssert(result == LPResult::Feasible); + previousGoodValues.resize(variables.size()); + for (size_t i = 0; i < variables.size(); i++) + previousGoodValues[i] = variables[i].value; + variablesPotentiallyOutOfBounds.clear(); + } + else if (_trailSize < trailSize) + { +#ifdef DEBUG + cerr << "=== Backtracking from " << trailSize << " to " << _trailSize << endl; +#endif + while (!storedBounds.empty()) + { + auto&& [ts, varIndex, bounds, lowerReason, upperReason] = storedBounds.back(); + //TODO should this be "<"? + if (ts <= _trailSize) + break; + variables[varIndex].bounds = bounds; + variables[varIndex].lowerReason = lowerReason; + variables[varIndex].upperReason = upperReason; + // TODO I think this is not needed because of "previousGoodValues + // we can maybe assert it. + //variablesPotentiallyOutOfBounds.insert(varIndex); + storedBounds.pop_back(); + } + for (size_t i = 0; i < previousGoodValues.size(); i++) + variables.at(i).value = previousGoodValues[i]; + variablesPotentiallyOutOfBounds.clear(); + result = LPResult::Feasible; + } + trailSize = _trailSize; +} + +#ifdef DEBUG +void LPSolver::setVariableName(size_t _variable, string _name) +{ + size_t index = maybeAddOuterVariable(_variable); + variables[index].name = move(_name); +} +#else +void LPSolver::setVariableName(size_t, string) +{ +} +#endif + + + +pair LPSolver::check() +{ + // TODO below is an old comment - but maybe we can optimize something to that effect + // by moving functionality to 'activateConstraint'. + + // TODO one third of the computing time (inclusive) in this function + // is spent on "operator<" - maybe we can cache "is in bounds" for variables + // and invalidate that in the update procedures. + +#ifdef DEBUG + cerr << "checking..." << endl; + cerr << toString() << endl; + cerr << "----------------------------" << endl; +// cerr << "fixing non-basic..." << endl; +#endif + if (result == LPResult::Feasible) + return make_pair(LPResult::Feasible, std::set()); + result = nullopt; + // Adjust the assignments so we satisfy the bounds of the non-basic variables. + if (!correctNonbasic()) + { +#ifdef DEBUG + cerr << "---> infeasible" << endl; +#endif + result = LPResult::Infeasible; + return make_pair(LPResult::Infeasible, reasons); + } + + // Now try to make the basic variables happy, pivoting if necessary. + +#ifdef DEBUG +// cerr << "fixed non-basic." << endl; +// cerr << toString() << endl; +// cerr << "----------------------------" << endl; +#endif + + // TODO bound number of iterations + while (auto bvi = firstConflictingBasicVariable()) + { + Variable const& basicVar = variables[*bvi]; +#ifdef DEBUG +// cerr << "----------------------------" << endl; +// cerr << "Fixing basic " << basicVar.name << endl; +#endif + if (basicVar.bounds.lower && basicVar.bounds.upper) + solAssert(*basicVar.bounds.lower <= *basicVar.bounds.upper); + if (basicVar.bounds.lower && basicVar.value < *basicVar.bounds.lower) + { + if (auto replacementVar = firstReplacementVar(*bvi, true)) + { +#ifdef DEBUG +// cerr << "Replacing by " << variables[*replacementVar].name << endl; +// cerr << "Setting basic var to to " << basicVar.bounds.lower->m_main << endl; +#endif + + pivotAndUpdate(*bvi, *basicVar.bounds.lower, *replacementVar); + } + else + { +#ifdef DEBUG + cerr << "---> infeasible" << endl; +#endif + result = LPResult::Infeasible; + reasons = reasonsForUnsat(*bvi, true); + return make_pair(LPResult::Infeasible, reasons); + } + } + else if (basicVar.bounds.upper && basicVar.value > *basicVar.bounds.upper) + { + if (auto replacementVar = firstReplacementVar(*bvi, false)) + { +#ifdef DEBUG +// cerr << "Replacing by " << variables[*replacementVar].name << endl; +#endif + pivotAndUpdate(*bvi, *basicVar.bounds.upper, *replacementVar); + } + else + { +#ifdef DEBUG + cerr << "---> infeasible" << endl; +#endif + result = LPResult::Infeasible; + reasons = reasonsForUnsat(*bvi, false); + return make_pair(LPResult::Infeasible, reasons); + } + } +#ifdef DEBUG +// cerr << "Fixed basic " << basicVar.name << endl; +// cerr << toString() << endl; +#endif + } + + result = LPResult::Feasible; +#ifdef DEBUG + cerr << toString() << endl; + cerr << "---> FEAsible" << endl; +#endif + return make_pair(LPResult::Feasible, std::set()); +} + +string LPSolver::toString() const +{ + string resultString = "LP Solver state (trail size " + to_string(trailSize) + "):\n"; + auto varName = [&](size_t _i) { +#ifdef DEBUG + return variables[_i].name; +#else + return "x" + to_string(_i); +#endif + }; + for (auto&& [i, v]: variables | ranges::views::enumerate) + { + if (v.bounds.lower) + resultString += v.bounds.lower->toString() + " <= "; + else + resultString += " "; + resultString += varName(i); + if (v.bounds.upper) + resultString += " <= " + v.bounds.upper->toString(); + else + resultString += " "; + resultString += " := " + v.value.toString() + "\n"; + } + for (size_t rowIndex = 0; rowIndex < factors.rows(); rowIndex++) + { + string basicVarPrefix; + string rowString; + for (auto&& entry: const_cast(factors).iterateRow(rowIndex)) + { + rational const& f = entry.value; + solAssert(!!f); + size_t i = entry.col; + if (basicVariables.count(i) && basicVariables.at(i) == rowIndex) + { + solAssert(f == -1); + solAssert(basicVarPrefix.empty()); + basicVarPrefix = varName(i) + " = "; + } + else if (f != 0) + { + string joiner = f < 0 ? " - " : f > 0 && !rowString.empty() ? " + " : " "; + string factor = f == 1 || f == -1 ? "" : ::toString(abs(f)) + " "; + string var = varName(i); + rowString += joiner + factor + var; + } + } + resultString += basicVarPrefix + rowString + "\n"; + } + if (result) + { + if (*result == LPResult::Feasible) + resultString += "result: feasible\n"; + else + resultString += "result: infeasible\n"; + } + else + resultString += "result: unknown\n"; + + + return resultString + "----\n"; +} + +map LPSolver::model() const +{ + map result; +#ifdef DEBUG + for (auto&& [outerIndex, innerIndex]: varMapping) + // TODO assign proper value to "delta" + result[variables[innerIndex].name] = + variables[innerIndex].value.m_main + + variables[innerIndex].value.m_delta / rational(100000); +#endif + return result; +} + + +pair LPSolver::constraintIntoVariableBounds(Constraint const& _constraint) +{ + size_t numVariables = 0; + size_t latestVariableIndex = size_t(-1); + // Make all variables available and check if it is a simple bound on a variable. + for (auto const& [index, entry]: _constraint.data.enumerateTail()) + if (entry) + { + latestVariableIndex = index; + numVariables++; + if (!varMapping.count(index)) + addOuterVariable(index); + } + if (numVariables == 1) + { + // Add this as direct bound. + rational factor = _constraint.data[latestVariableIndex]; + RationalWithDelta bound = _constraint.data.front(); + if (_constraint.kind == Constraint::LESS_THAN) + bound -= RationalWithDelta::delta(); + bound /= factor; + Bounds bounds; + if (factor > 0 || _constraint.kind == Constraint::EQUAL) + bounds.upper = bound; + if (factor < 0 || _constraint.kind == Constraint::EQUAL) + bounds.lower = bound; + return make_pair(varMapping.at(latestVariableIndex), move(bounds)); + } + + // TODO do we need to introduce a slack variable if we have a (potentially new) + // non-basic variable, or if we have an equality constraint? + + // Introduce the slack variable. + size_t slackIndex = addNewVariable(); + // Name is only needed for printing +#ifdef DEBUG + variables[slackIndex].name = "_s" + to_string(m_slackVariableCounter++); +#endif + basicVariables[slackIndex] = factors.rows(); + + // Compress the constraint, i.e. turn outer variable indices into + // inner variable indices. + RationalWithDelta valueForSlack; + size_t row = factors.rows(); + // First, handle the basic variables. + for (auto const& [outerIndex, entry]: _constraint.data.enumerateTail()) + if (entry) + { + size_t innerIndex = varMapping.at(outerIndex); + if (basicVariables.count(innerIndex)) + { + factors.addMultipleOfRow( + basicVariables[innerIndex], + row, + entry + ); + factors.remove(factors.entry(row, innerIndex)); + } + } + // Now the non-basic. + for (auto const& [outerIndex, entry]: _constraint.data.enumerateTail()) + if (entry) + { + size_t innerIndex = varMapping.at(outerIndex); + if (!basicVariables.count(innerIndex)) + { + SparseMatrix::Entry& e = factors.entry(row, innerIndex); + e.value += entry; + if (!e.value) + factors.remove(e); + } + valueForSlack += variables[innerIndex].value * entry; + } + + factors.entry(row, slackIndex).value = -1; + + // TODO do we really not need to add this to "potentially out of bounds"? + + basicVariables[slackIndex] = row; + variables[slackIndex].value = valueForSlack; + + Bounds bounds; + if (_constraint.kind == Constraint::EQUAL) + bounds.lower = _constraint.data[0]; + bounds.upper = _constraint.data[0]; + if (_constraint.kind == Constraint::LESS_THAN) + *bounds.upper -= RationalWithDelta::delta(); + return make_pair(slackIndex, move(bounds)); +} + +void LPSolver::addBounds(size_t _variable, Bounds _bounds) +{ + Variable& var = variables[_variable]; + if (_bounds.lower && (!var.bounds.lower || *var.bounds.lower < *_bounds.lower)) + { + var.bounds.lower = move(_bounds.lower); + if (var.value < var.bounds.lower) + variablesPotentiallyOutOfBounds.insert(_variable); + } + if (_bounds.upper && (!var.bounds.upper || *var.bounds.upper > *_bounds.upper)) + { + var.bounds.upper = move(_bounds.upper); + if (var.value > var.bounds.upper) + variablesPotentiallyOutOfBounds.insert(_variable); + } +} + +set LPSolver::collectReasonsForVariable(size_t _variable) +{ + set reasons; + if (variables[_variable].lowerReason) + reasons.insert(*variables[_variable].lowerReason); + if (variables[_variable].upperReason) + reasons.insert(*variables[_variable].upperReason); + return reasons; +} + +void LPSolver::addOuterVariable(size_t _outerIndex) +{ + size_t index = addNewVariable(); + varMapping.emplace(_outerIndex, index); +} + +size_t LPSolver::maybeAddOuterVariable(size_t _outerIndex) +{ + if (varMapping.count(_outerIndex)) + return varMapping.at(_outerIndex); + size_t index = addNewVariable(); + varMapping.emplace(_outerIndex, index); + return index; +} + +size_t LPSolver::addNewVariable() +{ + size_t index = variables.size(); + variables.emplace_back(); + return index; +} + + +bool LPSolver::correctNonbasic() +{ + set toCorrect; + swap(toCorrect, variablesPotentiallyOutOfBounds); + for (size_t i: toCorrect) + { + Variable& var = variables.at(i); + if (var.bounds.lower && var.bounds.upper && *var.bounds.lower > *var.bounds.upper) + { + reasons = collectReasonsForVariable(i); + return false; + } + if (basicVariables.count(i)) + { + variablesPotentiallyOutOfBounds.insert(i); + continue; + } + if (!var.bounds.lower && !var.bounds.upper) + continue; + if (var.bounds.lower && var.value < *var.bounds.lower) + update(i, *var.bounds.lower); + else if (var.bounds.upper && var.value > *var.bounds.upper) + update(i, *var.bounds.upper); + } + return true; +} + +void LPSolver::update(size_t _varIndex, RationalWithDelta const& _value) +{ + RationalWithDelta delta = _value - variables[_varIndex].value; + variables[_varIndex].value = _value; + + // TODO can we store that? + map basicVarForRow = invertMap(basicVariables); + for (auto&& entry: factors.iterateColumn(_varIndex)) + if (entry.value && basicVarForRow.count(entry.row)) + { + size_t j = basicVarForRow[entry.row]; + variables[j].value += delta * entry.value; + //variablesPotentiallyOutOfBounds.insert(j); + } +} + +optional LPSolver::firstConflictingBasicVariable() const +{ + // TODO we could use "variablesPotentiallyOutOfBounds" here. + for (auto&& [i, row]: basicVariables) + { + Variable const& variable = variables[i]; + if ( + (variable.bounds.lower && variable.value < *variable.bounds.lower) || + (variable.bounds.upper && variable.value > *variable.bounds.upper) + ) + return i; + } + return nullopt; +} + +optional LPSolver::firstReplacementVar( + size_t _basicVarToReplace, + bool _increasing +) const +{ + for (auto&& entry: const_cast(factors).iterateRow(basicVariables.at(_basicVarToReplace))) + { + size_t i = entry.col; + rational const& factor = entry.value; + if (i == _basicVarToReplace || !factor) + continue; + bool positive = factor > 0; + if (!_increasing) + positive = !positive; + Variable const& candidate = variables.at(i); + if (positive && (!candidate.bounds.upper || candidate.value < *candidate.bounds.upper)) + return i; + if (!positive && (!candidate.bounds.lower || candidate.value > *candidate.bounds.lower)) + return i; + } + return nullopt; +} + +set LPSolver::reasonsForUnsat( + size_t _basicVarToReplace, + bool _increasing +) const +{ + set r; + if (_increasing && variables[_basicVarToReplace].lowerReason) + r.insert(*variables[_basicVarToReplace].lowerReason); + else if (!_increasing && variables[_basicVarToReplace].upperReason) + r.insert(*variables[_basicVarToReplace].upperReason); + + for (auto&& entry: const_cast(factors).iterateRow(basicVariables.at(_basicVarToReplace))) + { + size_t i = entry.col; + rational const& factor = entry.value; + if (i == _basicVarToReplace || !factor) + continue; + bool positive = factor > 0; + if (!_increasing) + positive = !positive; + Variable const& candidate = variables.at(i); + if (positive && candidate.upperReason) + r.insert(*candidate.upperReason); + if (!positive && candidate.lowerReason) + r.insert(*candidate.lowerReason); + } + return r; +} + +void LPSolver::pivot(size_t _old, size_t _new) +{ + // Transform pivotRow such that the coefficient for _new is -1 + // Then use that to set all other coefficients for _new to zero. + size_t pivotRow = basicVariables[_old]; + + rational pivot = factors.entry(pivotRow, _new).value; + solAssert(pivot != 0, ""); + if (pivot != -1) + factors.multiplyRowByFactor(pivotRow, rational{-1} / pivot); + + for (auto it = factors.iterateColumn(_new).begin(); it != factors.iterateColumn(_new).end(); ) + { + SparseMatrix::Entry& entry = *it; + // Increment becasue "addMultipleOfRow" might invalidate the iterator + ++it; + if (entry.row != pivotRow) + factors.addMultipleOfRow(pivotRow, entry.row, entry.value); + } + + + basicVariables.erase(_old); + basicVariables[_new] = pivotRow; +} + +void LPSolver::pivotAndUpdate( + size_t _oldBasicVar, + RationalWithDelta const& _newValue, + size_t _newBasicVar +) +{ + RationalWithDelta theta = (_newValue - variables[_oldBasicVar].value) / factors.entry(basicVariables[_oldBasicVar], _newBasicVar).value; + + variables[_oldBasicVar].value = _newValue; + variables[_newBasicVar].value += theta; + + // TODO can we store that? + map basicVarForRow = invertMap(basicVariables); + for (auto&& entry: factors.iterateColumn(_newBasicVar)) + if (basicVarForRow.count(entry.row)) + { + size_t i = basicVarForRow[entry.row]; + if (i != _oldBasicVar) + variables[i].value += theta * entry.value; + } + + pivot(_oldBasicVar, _newBasicVar); +} diff --git a/libsolutil/LPIncremental.h b/libsolutil/LPIncremental.h new file mode 100644 index 000000000..0216c9a79 --- /dev/null +++ b/libsolutil/LPIncremental.h @@ -0,0 +1,241 @@ +/* + This file is part of solidity. + + solidity is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + solidity is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with solidity. If not, see . +*/ +// SPDX-License-Identifier: GPL-3.0 +#pragma once + +//#define DEBUG 1 + +#include +#include + +#include + +#include +#include +#include + +namespace solidity::util +{ + +using rational = boost::rational; +using Model = std::map; +using ReasonSet = std::set; + +/** + * Constraint of the form + * - data[1] * x_1 + data[2] * x_2 + ... <= data[0] (LESS_OR_EQUAL) + * - data[1] * x_1 + data[2] * x_2 + ... < data[0] (LESS_THAN) + * - data[1] * x_1 + data[2] * x_2 + ... = data[0] (EQUAL) + * The set and order of variables is implied. + */ +struct Constraint +{ + LinearExpression data; + enum Kind { EQUAL, LESS_THAN, LESS_OR_EQUAL }; + Kind kind = LESS_OR_EQUAL; + + bool operator<(Constraint const& _other) const; + bool operator==(Constraint const& _other) const; +}; + +/** + * A two-dimensional rational number "a + b*delta" that can be used to perform strict comparisons: + * x > 0 is transformed into x >= 1*delta, where delta is assumed to be "small". Its value + * is never explicitly computed / set, it is just a symbolic parameter. + */ +struct RationalWithDelta +{ + RationalWithDelta(rational _x = {}): m_main(move(_x)) {} + static RationalWithDelta delta() + { + RationalWithDelta x(0); + x.m_delta = 1; + return x; + } + + RationalWithDelta& operator+=(RationalWithDelta const& _other) + { + m_main += _other.m_main; + m_delta += _other.m_delta; + return *this; + } + RationalWithDelta& operator-=(RationalWithDelta const& _other) + { + m_main -= _other.m_main; + m_delta -= _other.m_delta; + return *this; + } + RationalWithDelta operator-(RationalWithDelta const& _other) const + { + RationalWithDelta ret = *this; + ret -= _other; + return ret; + } + RationalWithDelta& operator*=(rational const& _factor) + { + m_main *= _factor; + m_delta *= _factor; + return *this; + } + RationalWithDelta operator*(rational const& _factor) const + { + RationalWithDelta ret = *this; + ret *= _factor; + return ret; + } + RationalWithDelta& operator/=(rational const& _factor) + { + m_main /= _factor; + m_delta /= _factor; + return *this; + } + RationalWithDelta operator/(rational const& _factor) const + { + RationalWithDelta ret = *this; + ret /= _factor; + return ret; + } + bool operator<=(RationalWithDelta const& _other) const + { + return std::tie(m_main, m_delta) <= std::tie(_other.m_main, _other.m_delta); + } + bool operator>=(RationalWithDelta const& _other) const + { + return std::tie(m_main, m_delta) >= std::tie(_other.m_main, _other.m_delta); + } + bool operator<(RationalWithDelta const& _other) const + { + return std::tie(m_main, m_delta) < std::tie(_other.m_main, _other.m_delta); + } + bool operator>(RationalWithDelta const& _other) const + { + return std::tie(m_main, m_delta) > std::tie(_other.m_main, _other.m_delta); + } + bool operator==(RationalWithDelta const& _other) const + { + return std::tie(m_main, m_delta) == std::tie(_other.m_main, _other.m_delta); + } + bool operator!=(RationalWithDelta const& _other) const + { + return std::tie(m_main, m_delta) != std::tie(_other.m_main, _other.m_delta); + } + + std::string toString() const; + + rational m_main; + rational m_delta; +}; + +} + + +namespace solidity::util +{ + +enum class LPResult +{ + Unknown, + Unbounded, ///< System has a solution, but it can have an arbitrary objective value. + Feasible, ///< System has a solution (it might be unbounded, though). + Infeasible ///< System does not have any solution. +}; + +class LPSolver +{ +public: + void addConstraint(Constraint const& _constraint); + void addLowerBound(size_t _variable, RationalWithDelta _bound); + void addUpperBound(size_t _variable, RationalWithDelta _bound); + /// Add the conditional constraint but do not activate it yet. + void addConditionalConstraint(Constraint const& _constraint, size_t _reason); + void activateConstraint(size_t _reason); + void setTrailSize(size_t _trailSize); + + void setVariableName(size_t _variable, std::string _name); + + std::pair check(); + + std::string toString() const; + std::map model() const; + +private: + struct Bounds + { + std::optional lower; + std::optional upper; + }; + struct Variable + { +#ifdef DEBUG + std::string name = {}; +#endif + RationalWithDelta value = {}; + Bounds bounds = {}; + std::optional lowerReason; + std::optional upperReason; + }; + + /// Consumes a constraint and returns a controlling variable (can be a new slack + /// but does not need to) and corresponding bounds. + /// If it adds a slack variable, updates the factors and properly sets the value + /// for the slack variable (which will be a new basic variable). + std::pair constraintIntoVariableBounds(Constraint const& _constraint); + void addBounds(size_t _variable, Bounds _bounds); + std::set collectReasonsForVariable(size_t _variable); + + bool correctNonbasic(); + /// Set value of non-basic variable. + void update(size_t _varIndex, RationalWithDelta const& _value); + /// @returns the index of the first basic variable violating its bounds. + std::optional firstConflictingBasicVariable() const; + std::optional firstReplacementVar(size_t _basicVarToReplace, bool _increasing) const; + /// @returns the set of reasons in case "firstReplacementVar" failed. + std::set reasonsForUnsat(size_t _basicVarToReplace, bool _increasing) const; + + void pivot(size_t _old, size_t _new); + void pivotAndUpdate(size_t _oldBasicVar, RationalWithDelta const& _newValue, size_t _newBasicVar); + + void addOuterVariable(size_t _outerIndex); + /// Adds a new outer variable if it is not known yet and returns the inner index in any case. + size_t maybeAddOuterVariable(size_t _outerIndex); + size_t addNewVariable(); + + /// Counter to enable unique names for the slack variables. + size_t m_slackVariableCounter = 0; + std::optional result = std::nullopt; + size_t trailSize = 0; + SparseMatrix factors; + std::vector variables; + /// Stack of (trail size, variable index, bounds, lower reason, upper reason). + std::vector, std::optional>> storedBounds; + /// Last known satisfying values for variables. + std::vector previousGoodValues; + std::set variablesPotentiallyOutOfBounds; + /// Variable index to constraint it controls. + std::map basicVariables; + /// Maps outer indices to inner indices. + std::map varMapping = {}; + /// Mapping from reason (constraint ID) to variable it controls and bounds for it. + /// A variable can be controlled by multiple constraints. + /// TODO do we want to store the reverse mapping? + std::map> reasonToBounds; + std::set reasons; + + +}; + +} diff --git a/test/libsolutil/BooleanLP.cpp b/test/libsolutil/BooleanLP.cpp index 3346eb6b6..6b77586de 100644 --- a/test/libsolutil/BooleanLP.cpp +++ b/test/libsolutil/BooleanLP.cpp @@ -267,7 +267,8 @@ BOOST_AUTO_TEST_CASE(magic_square_3) solver.addAssertion(1 <= var && var <= 9); for (size_t i = 0; i < 9; i++) for (size_t j = i + 1; j < 9; j++) - solver.addAssertion(vars[i] != vars[j]); + //solver.addAssertion(vars[i] != vars[j]); + solver.addAssertion(vars[i] <= vars[j] - 1 || vars[i] >= vars[j] + 1); for (size_t i = 0; i < 3; i++) solver.addAssertion(vars[i] + vars[i + 3] + vars[i + 6] == sum); for (size_t i = 0; i < 9; i += 3) @@ -292,7 +293,7 @@ BOOST_AUTO_TEST_CASE(magic_square_4) solver.addAssertion(1 <= var && var <= 16); for (size_t i = 0; i < 16; i++) for (size_t j = i + 1; j < 16; j++) - solver.addAssertion(vars[i] != vars[j]); + solver.addAssertion(vars[i] <= vars[j] - 1 || vars[i] >= vars[j] + 1); for (size_t i = 0; i < 4; i++) solver.addAssertion(vars[i] + vars[i + 4] + vars[i + 8] + vars[i + 12] == sum); for (size_t i = 0; i < 16; i += 4) diff --git a/test/libsolutil/LP.cpp b/test/libsolutil/LP.cpp index ebc399d29..1e69abb45 100644 --- a/test/libsolutil/LP.cpp +++ b/test/libsolutil/LP.cpp @@ -16,7 +16,13 @@ */ // SPDX-License-Identifier: GPL-3.0 +#define LPIncremental 1 + +#if LPIncremental +#include +#else #include +#endif #include #include #include @@ -25,6 +31,8 @@ #include +#if 0 + using namespace std; using namespace solidity::smtutil; using namespace solidity::util; @@ -496,3 +504,4 @@ BOOST_AUTO_TEST_CASE(fuzzer2) BOOST_AUTO_TEST_SUITE_END() } +#endif