From 02695da9b70c9ebaea27fb226c4a065bb10f5d25 Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Wed, 16 Dec 2020 12:33:36 +0100 Subject: [PATCH] Fixed several issues regarding powers with negative exponents. --- .../3rdparty/sylvan/src/storm_wrapper.cpp | 4 +- .../expressions/ToRationalFunctionVisitor.cpp | 14 +++---- .../expressions/ToRationalNumberVisitor.cpp | 18 +++----- src/storm/utility/constants.cpp | 41 ++++++++++++++----- src/storm/utility/constants.h | 2 +- 5 files changed, 44 insertions(+), 35 deletions(-) diff --git a/resources/3rdparty/sylvan/src/storm_wrapper.cpp b/resources/3rdparty/sylvan/src/storm_wrapper.cpp index 5bfe580de..b0a10c1b8 100644 --- a/resources/3rdparty/sylvan/src/storm_wrapper.cpp +++ b/resources/3rdparty/sylvan/src/storm_wrapper.cpp @@ -215,8 +215,8 @@ storm_rational_number_ptr storm_rational_number_pow(storm_rational_number_ptr a, storm::RationalNumber const& srn_a = *(storm::RationalNumber const*)a; storm::RationalNumber const& srn_b = *(storm::RationalNumber const*)b; - carl::uint exponentAsInteger = carl::toInt(srn_b); - storm::RationalNumber* result_srn = new storm::RationalNumber(carl::pow(srn_a, exponentAsInteger)); + carl::sint exponentAsInteger = carl::toInt(srn_b); + storm::RationalNumber* result_srn = new storm::RationalNumber(storm::utility::pow(srn_a, exponentAsInteger)); return (storm_rational_number_ptr)result_srn; } diff --git a/src/storm/storage/expressions/ToRationalFunctionVisitor.cpp b/src/storm/storage/expressions/ToRationalFunctionVisitor.cpp index 21aeddce8..26a170836 100644 --- a/src/storm/storage/expressions/ToRationalFunctionVisitor.cpp +++ b/src/storm/storage/expressions/ToRationalFunctionVisitor.cpp @@ -39,25 +39,21 @@ namespace storm { boost::any ToRationalFunctionVisitor::visit(BinaryNumericalFunctionExpression const& expression, boost::any const& data) { RationalFunctionType firstOperandAsRationalFunction = boost::any_cast(expression.getFirstOperand()->accept(*this, data)); RationalFunctionType secondOperandAsRationalFunction = boost::any_cast(expression.getSecondOperand()->accept(*this, data)); - uint_fast64_t exponentAsInteger = 0; switch(expression.getOperatorType()) { case BinaryNumericalFunctionExpression::OperatorType::Plus: return firstOperandAsRationalFunction + secondOperandAsRationalFunction; - break; case BinaryNumericalFunctionExpression::OperatorType::Minus: return firstOperandAsRationalFunction - secondOperandAsRationalFunction; - break; case BinaryNumericalFunctionExpression::OperatorType::Times: return firstOperandAsRationalFunction * secondOperandAsRationalFunction; - break; case BinaryNumericalFunctionExpression::OperatorType::Divide: return firstOperandAsRationalFunction / secondOperandAsRationalFunction; - break; case BinaryNumericalFunctionExpression::OperatorType::Power: - STORM_LOG_THROW(storm::utility::isInteger(secondOperandAsRationalFunction), storm::exceptions::InvalidArgumentException, "Exponent of power operator must be a positive integer but is " << secondOperandAsRationalFunction << "."); - exponentAsInteger = storm::utility::convertNumber(secondOperandAsRationalFunction); - return storm::utility::pow(firstOperandAsRationalFunction, exponentAsInteger); - break; + { + STORM_LOG_THROW(storm::utility::isInteger(secondOperandAsRationalFunction), storm::exceptions::InvalidArgumentException, "Exponent of power operator must be an integer but is " << secondOperandAsRationalFunction << "."); + auto exponentAsInteger = storm::utility::convertNumber(secondOperandAsRationalFunction); + return storm::utility::pow(firstOperandAsRationalFunction, exponentAsInteger); + } default: STORM_LOG_ASSERT(false, "Illegal operator type."); } diff --git a/src/storm/storage/expressions/ToRationalNumberVisitor.cpp b/src/storm/storage/expressions/ToRationalNumberVisitor.cpp index 24c50e32f..0509d9677 100644 --- a/src/storm/storage/expressions/ToRationalNumberVisitor.cpp +++ b/src/storm/storage/expressions/ToRationalNumberVisitor.cpp @@ -48,38 +48,32 @@ namespace storm { RationalNumberType firstOperandAsRationalNumber = boost::any_cast(expression.getFirstOperand()->accept(*this, data)); RationalNumberType secondOperandAsRationalNumber = boost::any_cast(expression.getSecondOperand()->accept(*this, data)); RationalNumberType result; - uint_fast64_t exponentAsInteger; switch(expression.getOperatorType()) { case BinaryNumericalFunctionExpression::OperatorType::Plus: result = firstOperandAsRationalNumber + secondOperandAsRationalNumber; return result; - break; case BinaryNumericalFunctionExpression::OperatorType::Minus: result = firstOperandAsRationalNumber - secondOperandAsRationalNumber; return result; - break; case BinaryNumericalFunctionExpression::OperatorType::Times: result = firstOperandAsRationalNumber * secondOperandAsRationalNumber; return result; - break; case BinaryNumericalFunctionExpression::OperatorType::Divide: result = firstOperandAsRationalNumber / secondOperandAsRationalNumber; return result; - break; case BinaryNumericalFunctionExpression::OperatorType::Min: result = std::min(firstOperandAsRationalNumber, secondOperandAsRationalNumber); return result; - break; case BinaryNumericalFunctionExpression::OperatorType::Max: result = std::max(firstOperandAsRationalNumber, secondOperandAsRationalNumber); return result; - break; case BinaryNumericalFunctionExpression::OperatorType::Power: - STORM_LOG_THROW(storm::utility::isInteger(secondOperandAsRationalNumber), storm::exceptions::InvalidArgumentException, "Exponent of power operator must be a positive integer."); - exponentAsInteger = storm::utility::convertNumber(secondOperandAsRationalNumber); - result = storm::utility::pow(firstOperandAsRationalNumber, exponentAsInteger); - return result; - break; + { + STORM_LOG_THROW(storm::utility::isInteger(secondOperandAsRationalNumber), storm::exceptions::InvalidArgumentException, "Exponent of power operator must be an integer."); + auto exponentAsInteger = storm::utility::convertNumber(secondOperandAsRationalNumber); + result = storm::utility::pow(firstOperandAsRationalNumber, exponentAsInteger); + return result; + } default: STORM_LOG_ASSERT(false, "Illegal operator type."); } diff --git a/src/storm/utility/constants.cpp b/src/storm/utility/constants.cpp index 9172807ad..8ebb6a758 100644 --- a/src/storm/utility/constants.cpp +++ b/src/storm/utility/constants.cpp @@ -224,7 +224,7 @@ namespace storm { } template - ValueType pow(ValueType const& value, uint_fast64_t exponent) { + ValueType pow(ValueType const& value, int_fast64_t exponent) { return std::pow(value, exponent); } @@ -458,13 +458,18 @@ namespace storm { } template<> - typename NumberTraits::IntegerType pow(typename NumberTraits::IntegerType const& value, uint_fast64_t exponent) { + typename NumberTraits::IntegerType pow(typename NumberTraits::IntegerType const& value, int_fast64_t exponent) { + STORM_LOG_THROW(exponent >= 0, storm::exceptions::InvalidArgumentException, "Tried to compute the power 'x^y' as an integer, but the exponent 'y' is negative."); return carl::pow(value, exponent); } template<> - ClnRationalNumber pow(ClnRationalNumber const& value, uint_fast64_t exponent) { - return carl::pow(value, exponent); + ClnRationalNumber pow(ClnRationalNumber const& value, int_fast64_t exponent) { + if (exponent >= 0) { + return carl::pow(value, exponent); + } else { + return storm::utility::one() / carl::pow(value, -exponent); + } } template<> @@ -655,13 +660,18 @@ namespace storm { } template<> - typename NumberTraits::IntegerType pow(typename NumberTraits::IntegerType const& value, uint_fast64_t exponent) { + typename NumberTraits::IntegerType pow(typename NumberTraits::IntegerType const& value, int_fast64_t exponent) { + STORM_LOG_THROW(exponent >= 0, storm::exceptions::InvalidArgumentException, "Tried to compute the power 'x^y' as an integer, but the exponent 'y' is negative."); return carl::pow(value, exponent); } template<> - GmpRationalNumber pow(GmpRationalNumber const& value, uint_fast64_t exponent) { - return carl::pow(value, exponent); + GmpRationalNumber pow(GmpRationalNumber const& value, int_fast64_t exponent) { + if (exponent >= 0) { + return carl::pow(value, exponent); + } else { + return storm::utility::one() / carl::pow(value, -exponent); + } } template<> @@ -776,6 +786,11 @@ namespace storm { return carl::toInt(convertNumber(func)); } + template<> + carl::sint convertNumber(RationalFunction const& func){ + return carl::toInt(convertNumber(func)); + } + template<> double convertNumber(RationalFunction const& func) { return carl::toDouble(convertNumber(func)); @@ -861,8 +876,12 @@ namespace storm { } template<> - RationalFunction pow(RationalFunction const& value, uint_fast64_t exponent) { - return carl::pow(value, exponent); + RationalFunction pow(RationalFunction const& value, int_fast64_t exponent) { + if (exponent >= 0) { + return carl::pow(value, exponent); + } else { + return storm::utility::one() / carl::pow(value, -exponent); + } } template<> @@ -912,7 +931,7 @@ namespace storm { template std::pair minmax(std::map const&); template double minimum(std::map const&); template double maximum(std::map const&); - template double pow(double const& value, uint_fast64_t exponent); + template double pow(double const& value, int_fast64_t exponent); template double max(double const& first, double const& second); template double min(double const& first, double const& second); template double sqrt(double const& number); @@ -945,7 +964,7 @@ namespace storm { template std::pair minmax(std::map const&); template float minimum(std::map const&); template float maximum(std::map const&); - template float pow(float const& value, uint_fast64_t exponent); + template float pow(float const& value, int_fast64_t exponent); template float max(float const& first, float const& second); template float min(float const& first, float const& second); template float sqrt(float const& number); diff --git a/src/storm/utility/constants.h b/src/storm/utility/constants.h index 787ead089..e5fceec44 100644 --- a/src/storm/utility/constants.h +++ b/src/storm/utility/constants.h @@ -137,7 +137,7 @@ namespace storm { ValueType maximum(std::map const& values); template - ValueType pow(ValueType const& value, uint_fast64_t exponent); + ValueType pow(ValueType const& value, int_fast64_t exponent); template ValueType max(ValueType const& first, ValueType const& second);