diff --git a/libsolutil/CMakeLists.txt b/libsolutil/CMakeLists.txt index e055317be..a1d915bb3 100644 --- a/libsolutil/CMakeLists.txt +++ b/libsolutil/CMakeLists.txt @@ -23,9 +23,12 @@ set(sources Keccak256.h LazyInit.h LEB128.h + LP.cpp + LP.h Numeric.cpp Numeric.h picosha2.h + LinearExpression.h Result.h SetOnce.h StringUtils.cpp diff --git a/libsolutil/LP.cpp b/libsolutil/LP.cpp new file mode 100644 index 000000000..e7e47fb5e --- /dev/null +++ b/libsolutil/LP.cpp @@ -0,0 +1,778 @@ +/* + 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 + +using namespace std; +using namespace solidity; +using namespace solidity::util; + +using rational = boost::rational; + + +namespace +{ + +/** + * Simplex tableau. + */ +struct Tableau +{ + /// The factors of the objective function (first row of the tableau) + LinearExpression objective; + /// The tableau matrix (equational forrm). + std::vector data; +}; + + +/// Adds slack variables to remove non-equality costraints from a set of constraints. +/// The second return variable is true if new variables have been added. +pair, bool> toEquationalForm(vector _constraints) +{ + size_t varsNeeded = static_cast(ranges::count_if(_constraints, [](Constraint const& _c) { return !_c.equality; })); + if (varsNeeded > 0) + { + size_t columns = _constraints.at(0).data.size(); + size_t currentVariable = 0; + for (Constraint& constraint: _constraints) + { + solAssert(constraint.data.size() == columns, ""); + constraint.data.factors += vector(varsNeeded, rational{}); + if (!constraint.equality) + { + constraint.equality = true; + constraint.data[columns + currentVariable] = bigint(1); + currentVariable++; + } + } + } + + return make_pair(move(_constraints), varsNeeded > 0); +} + +optional findPivotColumn(Tableau const& _tableau) +{ + auto&& [maxColumn, maxValue] = ranges::max( + _tableau.objective | ranges::views::enumerate | ranges::views::tail, + {}, + [](std::pair const& _x) { return _x.second; } + ); + + if (maxValue <= rational{0}) + return nullopt; // found optimum + else + return maxColumn; +} + +optional findPivotRow(Tableau const& _tableau, size_t _pivotColumn) +{ + auto positiveColumnEntries = + ranges::views::iota(size_t(0), _tableau.data.size()) | + ranges::views::transform([&](size_t i) { + return make_pair(i, _tableau.data[i][_pivotColumn]); + }) | + ranges::views::filter([](pair const& _entry) { + return _entry.second > 0; + }); + if (positiveColumnEntries.empty()) + return nullopt; // unbounded + + return ranges::min( + positiveColumnEntries, + {}, + [&](std::pair const& _entry) { + return _tableau.data[_entry.first][0] / _entry.second; + } + ).first; +} + +/// Performs equivalence transform on @a _tableau, so that +/// the column @a _pivotColumn is all zeros except for @a _pivotRow, +/// where it is 1. +void performPivot(Tableau& _tableau, size_t _pivotRow, size_t _pivotColumn) +{ + rational pivot = _tableau.data[_pivotRow][_pivotColumn]; + solAssert(pivot != 0, ""); + if (pivot != 1) + _tableau.data[_pivotRow] /= pivot; + solAssert(_tableau.data[_pivotRow][_pivotColumn] == rational(1), ""); + + LinearExpression const& _pivotRowData = _tableau.data[_pivotRow]; + auto subtractPivotRow = [&](LinearExpression& _row) { + if (_row[_pivotColumn] == rational{1}) + _row -= _pivotRowData; + else if (_row[_pivotColumn] != rational{}) + _row -= _row[_pivotColumn] * _pivotRowData; + }; + + subtractPivotRow(_tableau.objective); + for (size_t i = 0; i < _tableau.data.size(); ++i) + if (i != _pivotRow) + subtractPivotRow(_tableau.data[i]); +} + +void selectLastVectorsAsBasis(Tableau& _tableau) +{ + // We might skip the operation for a column if it is already the correct + // unit vector and its cost coefficient is zero. + size_t columns = _tableau.objective.size(); + size_t rows = _tableau.data.size(); + for (size_t i = 0; i < rows; ++i) + performPivot(_tableau, i, columns - rows + i); +} + +/// If column @a _column inside tableau is a basis vector +/// (i.e. one entry is 1, the others are 0), returns the index +/// of the 1, otherwise nullopt. +optional basisVariable(Tableau const& _tableau, size_t _column) +{ + optional row; + for (size_t i = 0; i < _tableau.data.size(); ++i) + if (_tableau.data[i][_column] == bigint(1)) + { + if (row) + return std::nullopt; + else + row = i; + } + else if (_tableau.data[i][_column] != 0) + return std::nullopt; + return row; +} + +/// @returns a solution vector, assuming one exists. +/// The solution vector minimizes the objective function if the tableau +/// is the result of the simplex algorithm. +vector solutionVector(Tableau const& _tableau) +{ + vector result; + vector rowsSeen(_tableau.data.size(), false); + for (size_t j = 1; j < _tableau.objective.size(); j++) + { + optional row = basisVariable(_tableau, j); + if (row && rowsSeen[*row]) + row = nullopt; + result.emplace_back(row ? _tableau.data[*row][0] : rational{}); + if (row) + rowsSeen[*row] = true; + } + return result; +} + + +/// Solve the LP A x = b s.t. min c^T x +/// Here, c is _tableau.objective and the first column of _tableau.data +/// encodes b and the other columns encode A +/// Assumes the tableau has a trivial basic feasible solution. +pair simplexEq(Tableau _tableau) +{ + size_t const iterations = min(60, 50 + _tableau.objective.size() * 2); + for (size_t step = 0; step <= iterations; ++step) + { + optional pivotColumn = findPivotColumn(_tableau); + if (!pivotColumn) + return make_pair(LPResult::Feasible, move(_tableau)); + optional pivotRow = findPivotRow(_tableau, *pivotColumn); + if (!pivotRow) + return make_pair(LPResult::Unbounded, move(_tableau)); + performPivot(_tableau, *pivotRow, *pivotColumn); + } + return make_pair(LPResult::Unknown, Tableau{}); +} + +/// We add slack variables to find a basic feasible solution. +/// In particular, there is a slack variable for each row +/// which is weighted negatively. Setting the new slack +/// variables to one and all other variables to zero yields +/// a basic feasible solution. +/// If the optimal solution has all slack variables set to zero, +/// this is a basic feasible solution. Otherwise, the original +/// problem is infeasible. +/// This function returns the modified tableau with the original +/// objective function and the slack variables removed. +pair simplexPhaseI(Tableau _tableau) +{ + LinearExpression originalObjective = _tableau.objective; + + size_t rows = _tableau.data.size(); + size_t columns = _tableau.objective.size(); + for (size_t i = 0; i < rows; ++i) + { + if (_tableau.data[i][0] < 0) + _tableau.data[i] *= -1; + _tableau.data[i].factors += vector(rows, bigint{}); + _tableau.data[i][columns + i] = 1; + } + _tableau.objective.factors = + vector(columns, rational{}) + + vector(rows, rational{-1}); + + // This sets the objective factors of the slack variables + // to zero (and thus selects a basic feasible solution). + selectLastVectorsAsBasis(_tableau); + + LPResult result; + tie(result, _tableau) = simplexEq(move(_tableau)); + solAssert(result == LPResult::Feasible || result == LPResult::Unbounded, ""); + + vector optimum = solutionVector(_tableau); + + for (size_t i = columns - 1; i < optimum.size(); ++i) + if (optimum[i] != 0) + return make_pair(LPResult::Infeasible, Tableau{}); + + _tableau.objective = originalObjective; + for (auto& row: _tableau.data) + row.resize(columns); + + return make_pair(LPResult::Feasible, move(_tableau)); +} + +/// Returns true if the all-zero solution is not a solution for the tableau. +bool needsPhaseI(Tableau const& _tableau) +{ + for (auto const& row: _tableau.data) + if (row[0] < 0) + return true; + return false; +} + +/// Solve the LP Ax <= b s.t. min c^Tx +pair> simplex(vector _constraints, LinearExpression _objectives) +{ + Tableau tableau; + tableau.objective = move(_objectives); + bool hasEquations = false; + // TODO change toEquationalForm to directly return the tableau + tie(_constraints, hasEquations) = toEquationalForm(_constraints); + for (Constraint& c: _constraints) + tableau.data.emplace_back(move(c.data)); + tableau.objective.resize(tableau.data.at(0).size()); + + if (hasEquations || needsPhaseI(tableau)) + { + LPResult result; + tie(result, tableau) = simplexPhaseI(move(tableau)); + if (result == LPResult::Infeasible || result == LPResult::Unknown) + return make_pair(result, vector{}); + solAssert(result == LPResult::Feasible, ""); + } + // We know that the system is satisfiable and we know a solution, + // but it is not optimal. + LPResult result; + tie(result, tableau) = simplexEq(move(tableau)); + solAssert(result == LPResult::Feasible || result == LPResult::Unbounded, ""); + return make_pair(result, solutionVector(tableau)); +} + +/// Turns all bounds into constraints. +/// @returns false if the bounds make the state infeasible. +bool boundsToConstraints(SolvingState& _state) +{ + size_t columns = _state.variableNames.size(); + + // Turn bounds into constraints. + for (auto const& [index, bounds]: _state.bounds | ranges::views::enumerate | ranges::views::tail) + { + if (bounds[0] && bounds[1]) + { + if (*bounds[0] > *bounds[1]) + return false; + if (*bounds[0] == *bounds[1]) + { + vector c(columns); + c[0] = *bounds[0]; + c[index] = bigint(1); + _state.constraints.emplace_back(Constraint{move(c), true}); + continue; + } + } + if (bounds[0] && *bounds[0] > 0) + { + vector c(columns); + c[0] = -*bounds[0]; + c[index] = bigint(-1); + _state.constraints.emplace_back(Constraint{move(c), false}); + } + if (bounds[1]) + { + vector c(columns); + c[0] = *bounds[1]; + c[index] = bigint(1); + _state.constraints.emplace_back(Constraint{move(c), false}); + } + } + _state.bounds.clear(); + return true; +} + +template +void eraseIndices(T& _data, vector const& _indices) +{ + T result; + for (size_t i = 0; i < _data.size(); i++) + if (!_indices[i]) + result.emplace_back(move(_data[i])); + _data = move(result); +} + + +void removeColumns(SolvingState& _state, vector const& _columnsToRemove) +{ + eraseIndices(_state.bounds, _columnsToRemove); + for (Constraint& constraint: _state.constraints) + eraseIndices(constraint.data, _columnsToRemove); + eraseIndices(_state.variableNames, _columnsToRemove); +} + +bool extractDirectConstraints(SolvingState& _state, bool& _changed) +{ + // Turn constraints of the form ax <= b into an upper bound on x. + vector constraintsToRemove(_state.constraints.size(), false); + bool needsRemoval = false; + for (auto const& [index, constraint]: _state.constraints | ranges::views::enumerate) + { + auto nonzero = constraint.data | ranges::views::enumerate | ranges::views::tail | ranges::views::filter( + [](std::pair const& _x) { return !!_x.second; } + ); + // TODO we can exit early on in the loop above since we only care about zero, one or more than one nonzero entries. + // TODO could also use iterators and exit if we can advance it twice. + auto numNonzero = ranges::distance(nonzero); + if (numNonzero > 1) + continue; + constraintsToRemove[index] = true; + needsRemoval = true; + if (numNonzero == 0) + { + // 0 <= b or 0 = b + if ( + constraint.data.factors.front() < 0 || + (constraint.equality && constraint.data.factors.front() != 0) + ) + return false; // Infeasible. + } + else + { + auto&& [varIndex, factor] = nonzero.front(); + // a * x <= b + rational bound = constraint.data[0] / factor; + if ( + (factor >= 0 || constraint.equality) && + (!_state.bounds[varIndex][1] || bound < _state.bounds[varIndex][1]) + ) + _state.bounds[varIndex][1] = bound; + if ( + (factor <= 0 || constraint.equality) && + (!_state.bounds[varIndex][0] || bound > _state.bounds[varIndex][0]) + ) + // Lower bound must be at least zero. + _state.bounds[varIndex][0] = max(rational{}, bound); + } + } + if (needsRemoval) + { + _changed = true; + eraseIndices(_state.constraints, constraintsToRemove); + } + return true; +} + +bool removeFixedVariables(SolvingState& _state, map& _model, bool& _changed) +{ + // Remove variables that have equal lower and upper bound. + for (auto const& [index, bounds]: _state.bounds | ranges::views::enumerate) + { + if (!bounds[1] || (!bounds[0] && bounds[1]->numerator() > 0)) + continue; + // Lower bound must be at least zero. + rational lower = max(rational{}, bounds[0] ? *bounds[0] : rational{}); + rational upper = *bounds[1]; + if (upper < lower) + return false; // Infeasible. + if (upper != lower) + continue; + _model[_state.variableNames.at(index)] = lower; + _state.bounds[index] = {}; + _changed = true; + + // substitute variable + for (Constraint& constraint: _state.constraints) + if (constraint.data.factors.at(index) != 0) + { + constraint.data[0] -= constraint.data[index] * lower; + constraint.data[index] = 0; + } + } + + return true; +} + +bool removeEmptyColumns(SolvingState& _state, map& _model, bool& _changed) +{ + vector variablesSeen(_state.bounds.size(), false); + for (auto const& constraint: _state.constraints) + { + for (auto&& [index, factor]: constraint.data | ranges::views::enumerate | ranges::views::tail) + if (factor) + variablesSeen[index] = true; + } + + // TODO we could assert that any variable we remove does not have conflicting bounds. + // (We also remove the bounds). + + vector variablesToRemove(variablesSeen.size(), false); + bool needsRemoval = false; + for (auto&& [i, seen]: variablesSeen | ranges::views::enumerate | ranges::views::tail) + if (!seen) + { + variablesToRemove[i] = true; + needsRemoval = true; + // TODO actually it is unbounded if _state.bounds.at(i)[1] is nullopt. + if (_state.bounds.at(i)[0] || _state.bounds.at(i)[1]) + _model[_state.variableNames.at(i)] = + _state.bounds.at(i)[1] ? + *_state.bounds.at(i)[1] : + *_state.bounds.at(i)[0]; + } + if (needsRemoval) + { + _changed = true; + removeColumns(_state, variablesToRemove); + } + return true; +} + +auto nonZeroEntriesInColumn(SolvingState const& _state, size_t _column) +{ + return + _state.constraints | + ranges::views::enumerate | + ranges::views::filter([=](auto const& _entry) { return _entry.second.data[_column] != 0; }) | + ranges::views::transform([](auto const& _entry) { return _entry.first; }); +} + +pair, vector> connectedComponent(SolvingState const& _state, size_t _column) +{ + solAssert(_state.variableNames.size() >= 2, ""); + + vector includedColumns(_state.variableNames.size(), false); + vector includedRows(_state.constraints.size(), false); + stack columnsToProcess; + columnsToProcess.push(_column); + while (!columnsToProcess.empty()) + { + size_t column = columnsToProcess.top(); + columnsToProcess.pop(); + if (includedColumns[column]) + continue; + includedColumns[column] = true; + + for (size_t row: nonZeroEntriesInColumn(_state, column)) + { + if (includedRows[row]) + continue; + includedRows[row] = true; + for (auto const& [index, entry]: _state.constraints[row].data | ranges::views::enumerate | ranges::views::tail) + if (entry && !includedColumns[index]) + columnsToProcess.push(index); + } + } + return make_pair(move(includedColumns), move(includedRows)); +} + +struct ProblemSplitter +{ + ProblemSplitter(SolvingState const& _state): + state(_state), + column(1), + seenColumns(vector(state.variableNames.size(), false)) + {} + + operator bool() const + { + return column < state.variableNames.size(); + } + + SolvingState next() + { + vector includedColumns; + vector includedRows; + tie(includedColumns, includedRows) = connectedComponent(state, column); + + // Update state. + seenColumns |= includedColumns; + ++column; + while (column < state.variableNames.size() && seenColumns[column]) + ++column; + + // Happens in case of not removed empty column. + // Currently not happening because those are removed during the simplification stage. + // TODO If this is the case, we should actually also check the bounds. + if (includedRows.empty()) + return next(); + + SolvingState splitOff; + + splitOff.variableNames.emplace_back(); + splitOff.bounds.emplace_back(); + + for (auto&& [i, included]: includedColumns | ranges::views::enumerate | ranges::views::tail) + { + if (!included) + continue; + splitOff.variableNames.emplace_back(move(state.variableNames[i])); + splitOff.bounds.emplace_back(move(state.bounds[i])); + } + for (auto&& [i, included]: includedRows | ranges::views::enumerate) + { + if (!included) + continue; + Constraint splitRow{{}, state.constraints[i].equality}; + for (size_t j = 0; j < state.constraints[i].data.size(); j++) + if (j == 0 || includedColumns[j]) + splitRow.data.factors.push_back(state.constraints[i].data[j]); + splitOff.constraints.push_back(move(splitRow)); + } + + return splitOff; + } + + SolvingState const& state; + size_t column = 1; + vector seenColumns; +}; + + +/// Simplifies the solving state according to some rules (remove rows without variables, etc). +/// @returns false if the state is determined to be infeasible during this process. +bool simplifySolvingState(SolvingState& _state, map& _model) +{ + // - Constraints with exactly one nonzero coefficient represent "a x <= b" + // and thus are turned into bounds. + // - Constraints with zero nonzero coefficients are constant relations. + // If such a relation is false, answer "infeasible", otherwise remove the constraint. + // - Empty columns can be removed. + // - Variables with matching bounds can be removed from the problem by substitution. + + bool changed = true; + while (changed) + { + changed = false; + + if (!removeFixedVariables(_state, _model, changed)) + return false; + + if (!extractDirectConstraints(_state, changed)) + return false; + + if (!removeFixedVariables(_state, _model, changed)) + return false; + + if (!removeEmptyColumns(_state, _model, changed)) + return false; + } + + // TODO return the values selected for named variables in order to + // be used when returning the model. + return true; +} + +void normalizeRowLengths(SolvingState& _state) +{ + size_t vars = max(_state.variableNames.size(), _state.bounds.size()); + for (Constraint const& c: _state.constraints) + vars = max(vars, c.data.size()); + _state.variableNames.resize(vars); + _state.bounds.resize(vars); + for (Constraint& c: _state.constraints) + c.data.resize(vars); +} + +} + + +bool Constraint::operator<(Constraint const& _other) const +{ + if (equality != _other.equality) + return equality < _other.equality; + + for (size_t i = 0; i < max(data.size(), _other.data.size()); ++i) + { + rational const& a = data.get(i); + rational const& b = _other.data.get(i); + if (a != b) + return a < b; + } + return false; +} + +bool Constraint::operator==(Constraint const& _other) const +{ + if (equality != _other.equality) + return false; + + for (size_t i = 0; i < max(data.size(), _other.data.size()); ++i) + if (data.get(i) != _other.data.get(i)) + return false; + return true; +} + +bool SolvingState::operator<(SolvingState const& _other) const +{ + if (variableNames == _other.variableNames) + { + if (bounds == _other.bounds) + return constraints < _other.constraints; + else + return bounds < _other.bounds; + } + else + return variableNames < _other.variableNames; +} + +bool SolvingState::operator==(SolvingState const& _other) const +{ + return + variableNames == _other.variableNames && + bounds == _other.bounds && + constraints == _other.constraints; +} + +string SolvingState::toString() const +{ + string result; + + for (Constraint const& constraint: constraints) + { + vector line; + for (auto&& [index, multiplier]: constraint.data | ranges::views::enumerate) + if (index > 0 && multiplier != 0) + { + string mult = + multiplier == -1 ? + "-" : + multiplier == 1 ? + "" : + ::toString(multiplier) + " "; + line.emplace_back(mult + variableNames.at(index)); + } + result += + joinHumanReadable(line, " + ") + + (constraint.equality ? " = " : " <= ") + + ::toString(constraint.data.factors.front()) + + "\n"; + } + result += "Bounds:\n"; + for (auto&& [index, bounds]: bounds | ranges::views::enumerate) + { + if (!bounds[0] && !bounds[1]) + continue; + if (bounds[0]) + result += ::toString(*bounds[0]) + " <= "; + result += variableNames.at(index); + if (bounds[1]) + result += " <= " + ::toString(*bounds[1]); + result += "\n"; + } + return result; +} + + +pair> LPSolver::check(SolvingState _state) +{ + normalizeRowLengths(_state); + + map model; + + if (!simplifySolvingState(_state, model)) + return {LPResult::Infeasible, {}}; + + bool canOnlyBeUnknown = false; + ProblemSplitter splitter(_state); + while (splitter) + { + SolvingState split = splitter.next(); + solAssert(!split.constraints.empty(), ""); + solAssert(split.variableNames.size() >= 2, ""); + + LPResult lpResult; + vector solution; + auto it = m_cache.find(split); + if (it != m_cache.end()) + tie(lpResult, solution) = it->second; + else + { + SolvingState orig = split; + if (!boundsToConstraints(split)) + lpResult = LPResult::Infeasible; + else + { + LinearExpression objectives; + objectives.factors = + vector(1, rational(bigint(0))) + + vector(split.constraints.front().data.size() - 1, rational(bigint(1))); + tie(lpResult, solution) = simplex(split.constraints, move(objectives)); + } + m_cache.emplace(move(orig), make_pair(lpResult, solution)); + } + + switch (lpResult) + { + case LPResult::Feasible: + case LPResult::Unbounded: + break; + case LPResult::Infeasible: + return {LPResult::Infeasible, {}}; + case LPResult::Unknown: + // We do not stop here, because another independent query can still be infeasible. + canOnlyBeUnknown = true; + break; + } + for (auto&& [index, value]: solution | ranges::views::enumerate) + if (index + 1 < split.variableNames.size()) + model[split.variableNames.at(index + 1)] = value; + } + + if (canOnlyBeUnknown) + return {LPResult::Unknown, {}}; + + return {LPResult::Feasible, move(model)}; +} + diff --git a/libsolutil/LP.h b/libsolutil/LP.h new file mode 100644 index 000000000..e7f27ce92 --- /dev/null +++ b/libsolutil/LP.h @@ -0,0 +1,89 @@ +/* + 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 + +#include +#include + +#include + +#include +#include + +namespace solidity::util +{ + +/** + * Constraint of the form + * - data[1] * x_1 + data[2] * x_2 + ... <= data[0] (equality == false) + * - data[1] * x_1 + data[2] * x_2 + ... = data[0] (equality == true) + * The set and order of variables is implied. + */ +struct Constraint +{ + LinearExpression data; + bool equality = false; + + bool operator<(Constraint const& _other) const; + bool operator==(Constraint const& _other) const; +}; + +/** + * State used when solving an LP problem. + */ +struct SolvingState +{ + /// Names of variables, the index zero should be left empty. + /// TODO can we change that? + std::vector variableNames; + /// Lower and upper bounds for variables (in the sense of >= / <=). + std::vector>, 2>> bounds; + std::vector constraints; + + bool operator<(SolvingState const& _other) const; + bool operator==(SolvingState const& _other) const; + std::string toString() const; +}; + +enum class LPResult +{ + Unknown, + Unbounded, + Feasible, + Infeasible +}; + +/** + * LP solver for rational problems. + * + * Does not solve integer problems! + * + * Tries to split a given problem into sub-problems and utilizes a cache to quickly solve + * similar problems. + */ +class LPSolver +{ +public: + std::pair>> check(SolvingState _state); + +private: + // TODO check if the model is requested in production. If not, we do not need to cache it. + std::map>>> m_cache; +}; + +} diff --git a/libsolutil/LinearExpression.h b/libsolutil/LinearExpression.h new file mode 100644 index 000000000..0d038d260 --- /dev/null +++ b/libsolutil/LinearExpression.h @@ -0,0 +1,192 @@ +/* + 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 + +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include + + +namespace solidity::util +{ + +using rational = boost::rational; + +/** + * A linear expression of the form + * factors[0] + factors[1] * X1 + factors[2] * X2 + ... + */ +struct LinearExpression +{ + /// Creates the expression "_factor * X_index" + static LinearExpression factorForVariable(size_t _index, rational _factor) + { + LinearExpression result; + result.resizeAndSet(_index, move(_factor)); + return result; + } + + rational const& get(size_t _index) const + { + static rational const zero; + return _index < factors.size() ? factors[_index] : zero; + } + + rational const& operator[](size_t _index) const + { + return factors[_index]; + } + + rational& operator[](size_t _index) + { + return factors[_index]; + } + + auto begin() { return factors.begin(); } + auto end() { return factors.end(); } + + auto begin() const { return factors.begin(); } + auto end() const { return factors.end(); } + + void emplace_back(rational _value) { factors.emplace_back(move(_value)); } + + void resize(size_t _size) + { + factors.resize(_size); + } + + void resizeAndSet(size_t _index, rational _factor) + { + if (factors.size() <= _index) + factors.resize(_index + 1); + factors[_index] = move(_factor); + } + + bool isConstant() const + { + return ranges::all_of(factors | ranges::views::tail, [](rational const& _v) { return _v.numerator() == 0; }); + } + + size_t size() const { return factors.size(); } + + LinearExpression& operator/=(rational const& _divisor) + { + for (rational& x: factors) + if (x.numerator()) + x /= _divisor; + return *this; + } + + LinearExpression& operator*=(rational const& _factor) + { + for (rational& x: factors) + if (x.numerator()) + x *= _factor; + return *this; + } + + friend LinearExpression operator*(rational const& _factor, LinearExpression _expr) + { + for (rational& x: _expr.factors) + if (x.numerator()) + x *= _factor; + return _expr; + } + + LinearExpression& operator-=(LinearExpression const& _y) + { + if (size() < _y.size()) + factors.resize(_y.size()); + for (size_t i = 0; i < size(); ++i) + if (_y.factors[i].numerator()) + factors[i] -= _y.factors[i]; + return *this; + } + + LinearExpression operator-(LinearExpression const& _y) const + { + LinearExpression result = *this; + result -= _y; + return result; + } + + LinearExpression& operator+=(LinearExpression const& _y) + { + if (size() < _y.size()) + factors.resize(_y.size()); + for (size_t i = 0; i < size(); ++i) + if (_y.factors[i].numerator()) + factors[i] += _y.factors[i]; + return *this; + } + + LinearExpression operator+(LinearExpression const& _y) const + { + LinearExpression result = *this; + result += _y; + return result; + } + + + /// Multiply two vectors where the first element of each vector is a constant factor. + /// Only works if at most one of the vector has a nonzero element after the first. + /// If this condition is violated, returns nullopt. + static std::optional vectorProduct( + std::optional _x, + std::optional _y + ) + { + if (!_x || !_y) + return std::nullopt; + if (!_y->isConstant()) + swap(_x, _y); + if (!_y->isConstant()) + return std::nullopt; + + rational const& factor = _y->get(0); + + for (rational& element: _x->factors) + element *= factor; + return _x; + } + + std::vector factors; +}; + +// TODO + + +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; +} + +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 255218c15..c2332621a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -42,6 +42,7 @@ set(libsolutil_sources libsolutil/Keccak256.cpp libsolutil/LazyInit.cpp libsolutil/LEB128.cpp + libsolutil/LP.cpp libsolutil/StringUtils.cpp libsolutil/SwarmHash.cpp libsolutil/UTF8.cpp diff --git a/test/libsolutil/LP.cpp b/test/libsolutil/LP.cpp new file mode 100644 index 000000000..874111c46 --- /dev/null +++ b/test/libsolutil/LP.cpp @@ -0,0 +1,314 @@ +/* + 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 + +using namespace std; +using namespace solidity::smtutil; +using namespace solidity::util; + + +namespace solidity::util::test +{ + +class LPTestFramework +{ +public: + LPTestFramework() + { + m_solvingState.variableNames.emplace_back(""); + } + + LinearExpression constant(rational _value) + { + return LinearExpression::factorForVariable(0, _value); + } + + LinearExpression variable(string const& _name) + { + return LinearExpression::factorForVariable(variableIndex(_name), 1); + } + + /// Adds the constraint "_lhs <= _rhs". + void addLEConstraint(LinearExpression _lhs, LinearExpression _rhs) + { + _lhs -= _rhs; + _lhs[0] = -_lhs[0]; + m_solvingState.constraints.push_back({move(_lhs), false}); + } + + void addLEConstraint(LinearExpression _lhs, rational _rhs) + { + addLEConstraint(move(_lhs), constant(_rhs)); + } + + /// Adds the constraint "_lhs = _rhs". + void addEQConstraint(LinearExpression _lhs, LinearExpression _rhs) + { + _lhs -= _rhs; + _lhs[0] = -_lhs[0]; + m_solvingState.constraints.push_back({move(_lhs), true}); + } + + void addLowerBound(string _variable, rational _value) + { + size_t index = variableIndex(_variable); + if (index >= m_solvingState.bounds.size()) + m_solvingState.bounds.resize(index + 1); + m_solvingState.bounds.at(index)[0] = _value; + } + + void addUpperBound(string _variable, rational _value) + { + size_t index = variableIndex(_variable); + if (index >= m_solvingState.bounds.size()) + m_solvingState.bounds.resize(index + 1); + m_solvingState.bounds.at(index)[1] = _value; + } + + void feasible(vector> const& _solution) + { + auto [result, model] = m_solver.check(m_solvingState); + BOOST_REQUIRE(result == LPResult::Feasible); + for (auto const& [var, value]: _solution) + BOOST_CHECK_MESSAGE( + value == model.at(var), + var + " = "s + ::toString(model.at(var)) + " (expected " + ::toString(value) + ")" + ); + } + + void infeasible() + { + auto [result, model] = m_solver.check(m_solvingState); + BOOST_CHECK(result == LPResult::Infeasible); + } + +protected: + size_t variableIndex(string const& _name) + { + if (m_solvingState.variableNames.empty()) + m_solvingState.variableNames.emplace_back(""); + auto index = findOffset(m_solvingState.variableNames, _name); + if (!index) + { + index = m_solvingState.variableNames.size(); + m_solvingState.variableNames.emplace_back(_name); + } + return *index; + } + + LPSolver m_solver; + SolvingState m_solvingState; +}; + + +BOOST_FIXTURE_TEST_SUITE(LP, LPTestFramework, *boost::unit_test::label("nooptions")) + +BOOST_AUTO_TEST_CASE(basic) +{ + auto x = variable("x"); + addLEConstraint(2 * x, 10); + feasible({{"x", 5}}); +} + +BOOST_AUTO_TEST_CASE(not_linear_independent) +{ + addLEConstraint(2 * variable("x"), 10); + addLEConstraint(4 * variable("x"), 20); + feasible({{"x", 5}}); +} + +BOOST_AUTO_TEST_CASE(two_vars) +{ + addLEConstraint(variable("y"), 3); + addLEConstraint(variable("x"), 10); + addLEConstraint(variable("x") + variable("y"), 4); + feasible({{"x", 1}, {"y", 3}}); +} + +BOOST_AUTO_TEST_CASE(one_le_the_other) +{ + addLEConstraint(variable("x") + constant(2), variable("y") - constant(1)); + feasible({{"x", 0}, {"y", 3}}); +} + +BOOST_AUTO_TEST_CASE(factors) +{ + auto x = variable("x"); + auto y = variable("y"); + addLEConstraint(2 * y, 3); + addLEConstraint(16 * x, 10); + addLEConstraint(x + y, 4); + feasible({{"x", rational(5) / 8}, {"y", rational(3) / 2}}); +} + +BOOST_AUTO_TEST_CASE(cache) +{ + // This should use the cache already for the second part of the problem. + // We cannot really test that the cache has been used, but we can test + // that it results in the same value. + auto x = variable("x"); + auto y = variable("y"); + addLEConstraint(2 * y, 3); + addLEConstraint(2 * x, 3); + feasible({{"x", rational(3) / 2}, {"y", rational(3) / 2}}); + feasible({{"x", rational(3) / 2}, {"y", rational(3) / 2}}); +} + +BOOST_AUTO_TEST_CASE(bounds) +{ + addUpperBound("x", 200); + feasible({{"x", 200}}); + + addLEConstraint(variable("x"), 100); + feasible({{"x", 100}}); + + addLEConstraint(constant(5), variable("x")); + feasible({{"x", 100}}); + + addLowerBound("x", 20); + feasible({{"x", 100}}); + addLowerBound("x", 25); + feasible({{"x", 100}}); + + addUpperBound("x", 20); + infeasible(); +} + +BOOST_AUTO_TEST_CASE(bounds2) +{ + addLowerBound("x", 200); + addUpperBound("x", 250); + addLowerBound("y", 2); + addUpperBound("y", 3); + feasible({{"x", 250}, {"y", 3}}); + + addLEConstraint(variable("y"), variable("x")); + feasible({{"x", 250}, {"y", 3}}); + + addEQConstraint(variable("y") + constant(231), variable("x")); + feasible({{"x", 234}, {"y", 3}}); + + addEQConstraint(variable("y") + constant(10), variable("x") - variable("z")); + feasible({{"x", 234}, {"y", 3}}); + + addEQConstraint(variable("z") + variable("x"), constant(2)); + infeasible(); +} + +BOOST_AUTO_TEST_CASE(lower_bound) +{ + addLEConstraint(constant(1), variable("y")); + addLEConstraint(variable("x"), constant(10)); + addLEConstraint(2 * variable("x") + variable("y"), 2); + feasible({{"x", 0}, {"y", 2}}); +} + +BOOST_AUTO_TEST_CASE(check_infeasible) +{ + addLEConstraint(variable("x"), 3); + addLEConstraint(constant(5), variable("x")); + infeasible(); +} + +BOOST_AUTO_TEST_CASE(unbounded1) +{ + addLEConstraint(constant(2), variable("x")); + feasible({{"x", 2}}); +} + +BOOST_AUTO_TEST_CASE(unbounded2) +{ + auto x = variable("x"); + auto y = variable("y"); + addLEConstraint(constant(2), x + y); + addLEConstraint(x, 10); + feasible({{"x", 10}, {"y", 0}}); +} + +BOOST_AUTO_TEST_CASE(unbounded3) +{ + addLEConstraint(constant(0) - variable("x") - variable("y"), constant(10)); + feasible({{"x", 0}, {"y", 0}}); + + addLEConstraint(constant(0) - variable("x"), constant(10)); + feasible({{"x", 0}, {"y", 0}}); + + addEQConstraint(variable("y") + constant(3), variable("x")); + feasible({{"x", 3}, {"y", 0}}); + + addLEConstraint(variable("y") + variable("x"), constant(2)); + infeasible(); +} + + +BOOST_AUTO_TEST_CASE(equal) +{ + auto x = variable("x"); + auto y = variable("y"); + addEQConstraint(x, y + constant(10)); + addLEConstraint(x, 20); + feasible({{"x", 20}, {"y", 10}}); +} + + +BOOST_AUTO_TEST_CASE(equal_constant) +{ + auto x = variable("x"); + auto y = variable("y"); + addLEConstraint(x, y); + addEQConstraint(y, constant(5)); + feasible({{"x", 5}, {"y", 5}}); +} + +BOOST_AUTO_TEST_CASE(linear_dependent) +{ + auto x = variable("x"); + auto y = variable("y"); + auto z = variable("z"); + addLEConstraint(x, 5); + addLEConstraint(2 * y, 10); + addLEConstraint(3 * z, 15); + // Here, they should be split into three independent problems. + feasible({{"x", 5}, {"y", 5}, {"z", 5}}); + + addLEConstraint((x + y) + z, 100); + feasible({{"x", 5}, {"y", 5}, {"z", 5}}); + + addLEConstraint((x + y) + z, 2); + feasible({{"x", 2}, {"y", 0}, {"z", 0}}); + + addLEConstraint(constant(2), (x + y) + z); + feasible({{"x", 2}, {"y", 0}, {"z", 0}}); + + addEQConstraint(constant(2), (x + y) + z); + feasible({{"x", 2}, {"y", 0}, {"z", 0}}); +} + + + +BOOST_AUTO_TEST_SUITE_END() + +}