From 58ca07584d5f91e1b91189c37cf58b8806428eb3 Mon Sep 17 00:00:00 2001 From: dehnert Date: Wed, 4 Oct 2017 20:25:13 +0200 Subject: [PATCH] rational search for native linear equation solver and several involved fixes --- .../IterativeMinMaxLinearEquationSolver.cpp | 57 ++-- .../IterativeMinMaxLinearEquationSolver.h | 9 +- src/storm/solver/LinearEquationSolver.cpp | 1 + .../solver/NativeLinearEquationSolver.cpp | 260 +++++++++++++++++- src/storm/solver/NativeLinearEquationSolver.h | 17 +- src/storm/storage/expressions/Expression.cpp | 1 + src/storm/utility/KwekMehlhorn.h | 13 +- src/storm/utility/constants.cpp | 19 ++ src/storm/utility/constants.h | 6 + src/storm/utility/vector.h | 2 +- 10 files changed, 352 insertions(+), 33 deletions(-) diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index 975dec3ad..8f08a586c 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -262,15 +262,26 @@ namespace storm { MinMaxLinearEquationSolverRequirements IterativeMinMaxLinearEquationSolver::getRequirements(EquationSystemType const& equationSystemType, boost::optional const& direction) const { // Start by copying the requirements of the linear equation solver. MinMaxLinearEquationSolverRequirements requirements(this->linearEquationSolverFactory->getRequirements()); - + + // General requirements. if (this->getSettings().getSolutionMethod() == IterativeMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration) { requirements.requireLowerBounds(); + } else if (this->getSettings().getSolutionMethod() == IterativeMinMaxLinearEquationSolverSettings::SolutionMethod::RationalSearch) { + // As rational search needs to approach the solution from below, we cannot start with a valid scheduler hint as we would otherwise do. + // Instead, we need to require no end components. + if (equationSystemType == EquationSystemType::ReachabilityRewards) { + if (!direction || direction.get() == OptimizationDirection::Minimize) { + requirements.requireNoEndComponents(); + } + } } - // Guide requirements by whether or not we force soundness. if (this->getSettings().getForceSoundness()) { - // Only add requirements for value iteration here as the policy iteration requirements are indifferent if (this->getSettings().getSolutionMethod() == IterativeMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration) { + // Require both bounds now. + requirements.requireBounds(); + + // If we need to also converge from above, we cannot deal with end components in the notorious cases. if (equationSystemType == EquationSystemType::UntilProbabilities) { if (!direction || direction.get() == OptimizationDirection::Maximize) { requirements.requireNoEndComponents(); @@ -280,22 +291,29 @@ namespace storm { requirements.requireNoEndComponents(); } } + } else if (this->getSettings().getSolutionMethod() == IterativeMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration) { + if (equationSystemType == EquationSystemType::UntilProbabilities) { + if (!direction || direction.get() == OptimizationDirection::Maximize) { + requirements.requireValidInitialScheduler(); + } + } else if (equationSystemType == EquationSystemType::ReachabilityRewards) { + if (!direction || direction.get() == OptimizationDirection::Minimize) { + requirements.requireValidInitialScheduler(); + } + } } - - requirements.requireBounds(); - } - - // 'Regular' requirements (even for non-sound solving techniques). - if (equationSystemType == EquationSystemType::UntilProbabilities) { - if (this->getSettings().getSolutionMethod() == IterativeMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration) { - if (!direction || direction.get() == OptimizationDirection::Maximize) { + } else { + if (equationSystemType == EquationSystemType::UntilProbabilities) { + if (this->getSettings().getSolutionMethod() == IterativeMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration) { + if (!direction || direction.get() == OptimizationDirection::Maximize) { + requirements.requireValidInitialScheduler(); + } + } + } else if (equationSystemType == EquationSystemType::ReachabilityRewards) { + if (!direction || direction.get() == OptimizationDirection::Minimize) { requirements.requireValidInitialScheduler(); } } - } else if (equationSystemType == EquationSystemType::ReachabilityRewards) { - if (!direction || direction.get() == OptimizationDirection::Minimize) { - requirements.requireValidInitialScheduler(); - } } return requirements; @@ -677,6 +695,7 @@ namespace storm { template template typename std::enable_if::value && !NumberTraits::IsExact, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + // Version for when the overall value type is imprecise. // Create a rational representation of the input so we can check for a proper solution later. storm::storage::SparseMatrix rationalA = this->A->template toValueType(); @@ -711,6 +730,7 @@ namespace storm { template template typename std::enable_if::value && NumberTraits::IsExact, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + // Version for when the overall value type is exact and the same type is to be used for the imprecise part. if (!this->linEqSolverA) { this->linEqSolverA = this->linearEquationSolverFactory->create(*this->A); @@ -734,6 +754,9 @@ namespace storm { template template typename std::enable_if::value, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + // Version for when the overall value type is exact and the imprecise one is not. We first try to solve the + // problem using the imprecise data type and fall back to the exact type as needed. + // Translate A to its imprecise version. storm::storage::SparseMatrix impreciseA = this->A->template toValueType(); @@ -769,7 +792,7 @@ namespace storm { // Forward the call to the core rational search routine. converged = solveEquationsRationalSearchHelper(dir, impreciseSolver, *this->A, x, b, impreciseA, impreciseX, impreciseB, impreciseTmpX); } catch (storm::exceptions::PrecisionExceededException const& e) { - STORM_LOG_WARN("Precision of double was exceeded, trying to recover by switching to rational arithmetic."); + STORM_LOG_WARN("Precision of value type was exceeded, trying to recover by switching to rational arithmetic."); if (!auxiliaryRowGroupVector) { auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); @@ -859,7 +882,7 @@ namespace storm { // Count the iterations. overallIterations += result.iterations; - // Try to sharpen the result. + // Compute maximal precision until which to sharpen. uint64_t p = storm::utility::convertNumber(storm::utility::ceil(storm::utility::log10(storm::utility::one() / precision))); // Make sure that currentX and rationalX are not aliased. diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h index f7d062b4f..e3c9c8d0f 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h @@ -63,8 +63,6 @@ namespace storm { virtual MinMaxLinearEquationSolverRequirements getRequirements(EquationSystemType const& equationSystemType, boost::optional const& direction = boost::none) const override; private: - static bool isSolution(storm::OptimizationDirection dir, storm::storage::SparseMatrix const& matrix, std::vector const& values, std::vector const& b); - bool solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const; bool valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const; @@ -74,17 +72,16 @@ namespace storm { bool solveEquationsRationalSearch(OptimizationDirection dir, std::vector& x, std::vector const& b) const; template - static bool sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix const& A, std::vector const& x, std::vector const& b, std::vector& tmp); - + bool solveEquationsRationalSearchHelper(OptimizationDirection dir, IterativeMinMaxLinearEquationSolver const& impreciseSolver, storm::storage::SparseMatrix const& rationalA, std::vector& rationalX, std::vector const& rationalB, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector& tmpX) const; template typename std::enable_if::value && !NumberTraits::IsExact, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const; template typename std::enable_if::value && NumberTraits::IsExact, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const; template typename std::enable_if::value, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const; - template - bool solveEquationsRationalSearchHelper(OptimizationDirection dir, IterativeMinMaxLinearEquationSolver const& impreciseSolver, storm::storage::SparseMatrix const& rationalA, std::vector& rationalX, std::vector const& rationalB, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector& tmpX) const; + static bool sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix const& A, std::vector const& x, std::vector const& b, std::vector& tmp); + static bool isSolution(storm::OptimizationDirection dir, storm::storage::SparseMatrix const& matrix, std::vector const& values, std::vector const& b); void computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector& x, std::vector const& b, uint_fast64_t* choice = nullptr) const; diff --git a/src/storm/solver/LinearEquationSolver.cpp b/src/storm/solver/LinearEquationSolver.cpp index 87b4bc45c..e5655c07d 100644 --- a/src/storm/solver/LinearEquationSolver.cpp +++ b/src/storm/solver/LinearEquationSolver.cpp @@ -208,6 +208,7 @@ namespace storm { std::unique_ptr> GeneralLinearEquationSolverFactory::create() const { EquationSolverType equationSolver = storm::settings::getModule().getEquationSolver(); switch (equationSolver) { + case EquationSolverType::Native: return std::make_unique>(); case EquationSolverType::Elimination: return std::make_unique>(); default: return std::make_unique>(); } diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index 827966c5b..131ea8c0f 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -6,11 +6,14 @@ #include "storm/settings/modules/GeneralSettings.h" #include "storm/settings/modules/NativeEquationSolverSettings.h" +#include "storm/utility/ConstantsComparator.h" +#include "storm/utility/KwekMehlhorn.h" #include "storm/utility/constants.h" #include "storm/utility/vector.h" #include "storm/exceptions/InvalidStateException.h" #include "storm/exceptions/InvalidSettingsException.h" #include "storm/exceptions/UnmetRequirementException.h" +#include "storm/exceptions/PrecisionExceededException.h" namespace storm { namespace solver { @@ -641,6 +644,252 @@ namespace storm { return converged; } + template + bool NativeLinearEquationSolver::solveEquationsRationalSearch(std::vector& x, std::vector const& b) const { + return solveEquationsRationalSearchHelper(x, b); + } + + template + struct TemporaryHelper { + static std::vector* getTemporary(std::vector& rationalX, std::vector*& currentX, std::vector*& newX) { + return &rationalX; + } + + static void swapSolutions(std::vector& rationalX, std::vector*& rationalSolution, std::vector& x, std::vector*& currentX, std::vector*& newX) { + // Nothing to do. + } + }; + + template + struct TemporaryHelper { + static std::vector* getTemporary(std::vector& rationalX, std::vector*& currentX, std::vector*& newX) { + return newX; + } + + static void swapSolutions(std::vector& rationalX, std::vector*& rationalSolution, std::vector& x, std::vector*& currentX, std::vector*& newX) { + if (&rationalX == rationalSolution) { + // In this case, the rational solution is in place. + + // However, since the rational solution is no alias to current x, the imprecise solution is stored + // in current x and and rational x is not an alias to x, we can swap the contents of currentX to x. + std::swap(x, *currentX); + } else { + // Still, we may assume that the rational solution is not current x and is therefore new x. + std::swap(rationalX, *rationalSolution); + std::swap(x, *currentX); + } + } + }; + + template + template + bool NativeLinearEquationSolver::solveEquationsRationalSearchHelper(NativeLinearEquationSolver const& impreciseSolver, storm::storage::SparseMatrix const& rationalA, std::vector& rationalX, std::vector const& rationalB, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector& tmpX) const { + + std::vector* currentX = &x; + std::vector* newX = &tmpX; + + Status status = Status::InProgress; + uint64_t overallIterations = 0; + uint64_t valueIterationInvocations = 0; + ValueType precision = this->getSettings().getPrecision(); + impreciseSolver.startMeasureProgress(); + while (status == Status::InProgress && overallIterations < this->getSettings().getMaximalNumberOfIterations()) { + // Perform value iteration with the current precision. + typename NativeLinearEquationSolver::PowerIterationResult result = impreciseSolver.performPowerIteration(currentX, newX, b, storm::utility::convertNumber(precision), this->getSettings().getRelativeTerminationCriterion(), SolverGuarantee::LessOrEqual, overallIterations); + + // At this point, the result of the imprecise value iteration is stored in the (imprecise) current x. + + ++valueIterationInvocations; + STORM_LOG_TRACE("Completed " << valueIterationInvocations << " power iteration invocations, the last one with precision " << precision << " completed in " << result.iterations << " iterations."); + + // Count the iterations. + overallIterations += result.iterations; + + // Compute maximal precision until which to sharpen. + uint64_t p = storm::utility::convertNumber(storm::utility::ceil(storm::utility::log10(storm::utility::one() / precision))); + + // Make sure that currentX and rationalX are not aliased. + std::vector* temporaryRational = TemporaryHelper::getTemporary(rationalX, currentX, newX); + + // Sharpen solution and place it in the temporary rational. + bool foundSolution = sharpen(p, rationalA, *currentX, rationalB, *temporaryRational); + + // After sharpen, if a solution was found, it is contained in the free rational. + + if (foundSolution) { + status = Status::Converged; + + TemporaryHelper::swapSolutions(rationalX, temporaryRational, x, currentX, newX); + } else { + // Increase the precision. + precision = precision / 100; + } + } + + if (status == Status::InProgress && overallIterations == this->getSettings().getMaximalNumberOfIterations()) { + status = Status::MaximalIterationsExceeded; + } + + this->logIterations(status == Status::Converged, status == Status::TerminatedEarly, overallIterations); + + return status == Status::Converged || status == Status::TerminatedEarly; + } + + template + template + typename std::enable_if::value && !NumberTraits::IsExact, bool>::type NativeLinearEquationSolver::solveEquationsRationalSearchHelper(std::vector& x, std::vector const& b) const { + // Version for when the overall value type is imprecise. + + // Create a rational representation of the input so we can check for a proper solution later. + storm::storage::SparseMatrix rationalA = this->A->template toValueType(); + std::vector rationalX(x.size()); + std::vector rationalB = storm::utility::vector::convertNumericVector(b); + + if (!this->cachedRowVector) { + this->cachedRowVector = std::make_unique>(this->A->getRowCount()); + } + + // Forward the call to the core rational search routine. + bool converged = solveEquationsRationalSearchHelper(*this, rationalA, rationalX, rationalB, *this->A, x, b, *this->cachedRowVector); + + // Translate back rational result to imprecise result. + auto targetIt = x.begin(); + for (auto it = rationalX.begin(), ite = rationalX.end(); it != ite; ++it, ++targetIt) { + *targetIt = storm::utility::convertNumber(*it); + } + + if (!this->isCachingEnabled()) { + this->clearCache(); + } + + return converged; + } + + template + template + typename std::enable_if::value && NumberTraits::IsExact, bool>::type NativeLinearEquationSolver::solveEquationsRationalSearchHelper(std::vector& x, std::vector const& b) const { + // Version for when the overall value type is exact and the same type is to be used for the imprecise part. + + if (!this->linEqSolverA) { + this->linEqSolverA = this->linearEquationSolverFactory->create(*this->A); + this->linEqSolverA->setCachingEnabled(true); + } + + if (!this->cachedRowVector) { + this->cachedRowVector = std::make_unique>(this->A->getRowCount()); + } + + // Forward the call to the core rational search routine. + bool converged = solveEquationsRationalSearchHelper(*this, *this->A, x, b, *this->A, *this->cachedRowVector, b, x); + + if (!this->isCachingEnabled()) { + this->clearCache(); + } + + return converged; + } + + template + template + typename std::enable_if::value, bool>::type NativeLinearEquationSolver::solveEquationsRationalSearchHelper(std::vector& x, std::vector const& b) const { + // Version for when the overall value type is exact and the imprecise one is not. We first try to solve the + // problem using the imprecise data type and fall back to the exact type as needed. + + // Translate A to its imprecise version. + storm::storage::SparseMatrix impreciseA = this->A->template toValueType(); + + // Translate x to its imprecise version. + std::vector impreciseX(x.size()); + { + std::vector tmp(x.size()); + this->createLowerBoundsVector(tmp); + auto targetIt = impreciseX.begin(); + for (auto sourceIt = tmp.begin(); targetIt != impreciseX.end(); ++targetIt, ++sourceIt) { + *targetIt = storm::utility::convertNumber(*sourceIt); + } + } + + // Create temporary storage for an imprecise x. + std::vector impreciseTmpX(x.size()); + + // Translate b to its imprecise version. + std::vector impreciseB(b.size()); + auto targetIt = impreciseB.begin(); + for (auto sourceIt = b.begin(); targetIt != impreciseB.end(); ++targetIt, ++sourceIt) { + *targetIt = storm::utility::convertNumber(*sourceIt); + } + + // Create imprecise solver from the imprecise data. + NativeLinearEquationSolver impreciseSolver; + impreciseSolver.setMatrix(impreciseA); + impreciseSolver.setCachingEnabled(true); + + bool converged = false; + try { + // Forward the call to the core rational search routine. + converged = solveEquationsRationalSearchHelper(impreciseSolver, *this->A, x, b, impreciseA, impreciseX, impreciseB, impreciseTmpX); + } catch (storm::exceptions::PrecisionExceededException const& e) { + STORM_LOG_WARN("Precision of value type was exceeded, trying to recover by switching to rational arithmetic."); + + if (!this->cachedRowVector) { + this->cachedRowVector = std::make_unique>(this->A->getRowGroupCount()); + } + + // Translate the imprecise value iteration result to the one we are going to use from now on. + auto targetIt = this->cachedRowVector->begin(); + for (auto it = impreciseX.begin(), ite = impreciseX.end(); it != ite; ++it, ++targetIt) { + *targetIt = storm::utility::convertNumber(*it); + } + + // Get rid of the superfluous data structures. + impreciseX = std::vector(); + impreciseTmpX = std::vector(); + impreciseB = std::vector(); + impreciseA = storm::storage::SparseMatrix(); + + // Forward the call to the core rational search routine, but now with our value type as the imprecise value type. + converged = solveEquationsRationalSearchHelper(*this, *this->A, x, b, *this->A, *this->cachedRowVector, b, x); + } + + if (!this->isCachingEnabled()) { + this->clearCache(); + } + + return converged; + } + + template + template + bool NativeLinearEquationSolver::sharpen(uint64_t precision, storm::storage::SparseMatrix const& A, std::vector const& x, std::vector const& b, std::vector& tmp) { + for (uint64_t p = 0; p <= precision; ++p) { + storm::utility::kwek_mehlhorn::sharpen(p, x, tmp); + + if (NativeLinearEquationSolver::isSolution(A, tmp, b)) { + return true; + } + } + return false; + } + + template + bool NativeLinearEquationSolver::isSolution(storm::storage::SparseMatrix const& matrix, std::vector const& values, std::vector const& b) { + storm::utility::ConstantsComparator comparator; + + auto valueIt = values.begin(); + auto bIt = b.begin(); + for (uint64_t row = 0; row < matrix.getRowCount(); ++row, ++valueIt, ++bIt) { + ValueType rowValue = *bIt + matrix.multiplyRowWithVector(row, values); + + // If the value does not match the one in the values vector, the given vector is not a solution. + if (!comparator.isEqual(rowValue, *valueIt)) { + return false; + } + } + + // Checked all values at this point. + return true; + } + template void NativeLinearEquationSolver::logIterations(bool converged, bool terminate, uint64_t iterations) const { if (converged) { @@ -666,6 +915,8 @@ namespace storm { } else { return this->solveEquationsPower(x, b); } + } else if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::RationalSearch) { + return this->solveEquationsRationalSearch(x, b); } STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unknown solving technique."); @@ -738,7 +989,7 @@ namespace storm { template LinearEquationSolverProblemFormat NativeLinearEquationSolver::getEquationProblemFormat() const { - if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::Power) { + if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::Power || this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::RationalSearch) { return LinearEquationSolverProblemFormat::FixedPointSystem; } else { return LinearEquationSolverProblemFormat::EquationSystem; @@ -804,5 +1055,12 @@ namespace storm { template class NativeLinearEquationSolverSettings; template class NativeLinearEquationSolver; template class NativeLinearEquationSolverFactory; + +#ifdef STORM_HAVE_CARL + template class NativeLinearEquationSolverSettings; + template class NativeLinearEquationSolver; + template class NativeLinearEquationSolverFactory; + +#endif } } diff --git a/src/storm/solver/NativeLinearEquationSolver.h b/src/storm/solver/NativeLinearEquationSolver.h index c62d2d61f..5617244c8 100644 --- a/src/storm/solver/NativeLinearEquationSolver.h +++ b/src/storm/solver/NativeLinearEquationSolver.h @@ -7,6 +7,8 @@ #include "storm/solver/NativeMultiplier.h" +#include "storm/utility/NumberTraits.h" + namespace storm { namespace solver { @@ -38,7 +40,7 @@ namespace storm { private: bool forceSoundness; SolutionMethod method; - double precision; + ValueType precision; bool relative; uint_fast64_t maximalNumberOfIterations; ValueType omega; @@ -104,6 +106,19 @@ namespace storm { virtual bool solveEquationsWalkerChae(std::vector& x, std::vector const& b) const; virtual bool solveEquationsPower(std::vector& x, std::vector const& b) const; virtual bool solveEquationsSoundPower(std::vector& x, std::vector const& b) const; + virtual bool solveEquationsRationalSearch(std::vector& x, std::vector const& b) const; + + template + bool solveEquationsRationalSearchHelper(NativeLinearEquationSolver const& impreciseSolver, storm::storage::SparseMatrix const& rationalA, std::vector& rationalX, std::vector const& rationalB, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector& tmpX) const; + template + typename std::enable_if::value && !NumberTraits::IsExact, bool>::type solveEquationsRationalSearchHelper(std::vector& x, std::vector const& b) const; + template + typename std::enable_if::value && NumberTraits::IsExact, bool>::type solveEquationsRationalSearchHelper(std::vector& x, std::vector const& b) const; + template + typename std::enable_if::value, bool>::type solveEquationsRationalSearchHelper(std::vector& x, std::vector const& b) const; + template + static bool sharpen(uint64_t precision, storm::storage::SparseMatrix const& A, std::vector const& x, std::vector const& b, std::vector& tmp); + static bool isSolution(storm::storage::SparseMatrix const& matrix, std::vector const& values, std::vector const& b); // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted // when the solver is destructed. diff --git a/src/storm/storage/expressions/Expression.cpp b/src/storm/storage/expressions/Expression.cpp index 7f204d0c1..95c4f69e0 100644 --- a/src/storm/storage/expressions/Expression.cpp +++ b/src/storm/storage/expressions/Expression.cpp @@ -264,6 +264,7 @@ namespace storm { STORM_LOG_THROW(first.hasBooleanType(), storm::exceptions::InvalidTypeException, "Operator requires boolean operands."); return first; } + return Expression(std::shared_ptr(new BinaryBooleanFunctionExpression(first.getBaseExpression().getManager(), first.getType().logicalConnective(second.getType()), first.getBaseExpressionPointer(), second.getBaseExpressionPointer(), BinaryBooleanFunctionExpression::OperatorType::And))); } diff --git a/src/storm/utility/KwekMehlhorn.h b/src/storm/utility/KwekMehlhorn.h index b6e73bfeb..72038475d 100644 --- a/src/storm/utility/KwekMehlhorn.h +++ b/src/storm/utility/KwekMehlhorn.h @@ -11,17 +11,16 @@ namespace storm { template std::pair findRational(IntegerType const& alpha, IntegerType const& beta, IntegerType const& gamma, IntegerType const& delta) { - IntegerType alphaDivBetaFloor = alpha / beta; - IntegerType gammaDivDeltaFloor = gamma / delta; - IntegerType alphaModBeta = storm::utility::mod(alpha, beta); + std::pair alphaDivBetaPair = storm::utility::divide(alpha, beta); + std::pair gammaDivDeltaPair = storm::utility::divide(gamma, delta); - if (alphaDivBetaFloor == gammaDivDeltaFloor && !storm::utility::isZero(alphaModBeta)) { - std::pair subresult = findRational(delta, storm::utility::mod(gamma, delta), beta, alphaModBeta); - auto result = std::make_pair(alphaDivBetaFloor * subresult.first + subresult.second, subresult.first); + if (alphaDivBetaPair.first == gammaDivDeltaPair.first && !storm::utility::isZero(alphaDivBetaPair.second)) { + std::pair subresult = findRational(delta, gammaDivDeltaPair.second, beta, alphaDivBetaPair.second); + auto result = std::make_pair(alphaDivBetaPair.first * subresult.first + subresult.second, subresult.first); return result; } else { - auto result = std::make_pair(storm::utility::isZero(alphaModBeta) ? alphaDivBetaFloor : alphaDivBetaFloor + storm::utility::one(), storm::utility::one()); + auto result = std::make_pair(storm::utility::isZero(alphaDivBetaPair.second) ? alphaDivBetaPair.first : alphaDivBetaPair.first + storm::utility::one(), storm::utility::one()); return result; } } diff --git a/src/storm/utility/constants.cpp b/src/storm/utility/constants.cpp index dd9e3cae9..7fb973653 100644 --- a/src/storm/utility/constants.cpp +++ b/src/storm/utility/constants.cpp @@ -249,6 +249,11 @@ namespace storm { return std::fmod(first, second); } + template + std::pair divide(IntegerType const& dividend, IntegerType const& divisor) { + return std::make_pair(dividend / divisor, mod(dividend, divisor)); + } + template std::string to_string(ValueType const& value) { std::stringstream ss; @@ -386,6 +391,13 @@ namespace storm { return carl::mod(first, second); } + template<> + std::pair::IntegerType, typename NumberTraits::IntegerType> divide(typename NumberTraits::IntegerType const& dividend, typename NumberTraits::IntegerType const& divisor) { + std::pair::IntegerType, typename NumberTraits::IntegerType> result; + carl::divide(dividend, divisor, result.first, result.second); + return result; + } + template<> typename NumberTraits::IntegerType pow(typename NumberTraits::IntegerType const& value, uint_fast64_t exponent) { return carl::pow(value, exponent); @@ -557,6 +569,13 @@ namespace storm { return carl::mod(first, second); } + template<> + std::pair::IntegerType, typename NumberTraits::IntegerType> divide(typename NumberTraits::IntegerType const& dividend, typename NumberTraits::IntegerType const& divisor) { + std::pair::IntegerType, typename NumberTraits::IntegerType> result; + carl::divide(dividend, divisor, result.first, result.second); + return result; + } + template<> typename NumberTraits::IntegerType pow(typename NumberTraits::IntegerType const& value, uint_fast64_t exponent) { return carl::pow(value, exponent); diff --git a/src/storm/utility/constants.h b/src/storm/utility/constants.h index f2ce01183..b56573009 100644 --- a/src/storm/utility/constants.h +++ b/src/storm/utility/constants.h @@ -120,6 +120,12 @@ namespace storm { template typename NumberTraits::IntegerType denominator(RationalType const& number); + /*! + * (Integer-)Divides the dividend by the divisor and returns the result plus the remainder. + */ + template + std::pair divide(IntegerType const& dividend, IntegerType const& divisor); + template IntegerType mod(IntegerType const& first, IntegerType const& second); diff --git a/src/storm/utility/vector.h b/src/storm/utility/vector.h index 2bab57d47..3dbfa2e80 100644 --- a/src/storm/utility/vector.h +++ b/src/storm/utility/vector.h @@ -908,7 +908,7 @@ namespace storm { auto b2It = b2.begin(); for (; b1It != b1Ite; ++b1It, ++b2It) { - result += storm::utility::pow(*b1It - *b2It, 2); + result += storm::utility::pow(*b1It - *b2It, 2); } return result;