diff --git a/resources/3rdparty/sylvan/src/storm_wrapper.cpp b/resources/3rdparty/sylvan/src/storm_wrapper.cpp
index 78d73c0c3..74d191b95 100644
--- a/resources/3rdparty/sylvan/src/storm_wrapper.cpp
+++ b/resources/3rdparty/sylvan/src/storm_wrapper.cpp
@@ -9,6 +9,7 @@
#include "storm/adapters/RationalFunctionAdapter.h"
#include "storm/utility/constants.h"
+#include "storm/utility/KwekMehlhorn.h"
#include "storm/exceptions/InvalidOperationException.h"
#include <sylvan_config.h>
@@ -300,6 +301,25 @@ storm_rational_number_ptr storm_rational_number_ceil(storm_rational_number_ptr a
return (storm_rational_number_ptr)result_srn;
}
+storm_rational_number_ptr storm_double_sharpen(double value, size_t precision) {
+#ifndef RATIONAL_NUMBER_THREAD_SAFE
+ std::lock_guard<std::mutex> lock(rationalNumberMutex);
+#endif
+
+ storm::RationalNumber* result_srn = new storm::RationalNumber(storm::utility::kwek_mehlhorn::sharpen<storm::RationalNumber, double>(precision, value));
+ return (storm_rational_number_ptr)result_srn;
+}
+
+storm_rational_number_ptr storm_rational_number_sharpen(storm_rational_number_ptr a, size_t precision) {
+#ifndef RATIONAL_NUMBER_THREAD_SAFE
+ std::lock_guard<std::mutex> lock(rationalNumberMutex);
+#endif
+
+ storm::RationalNumber const& srn_a = *(storm::RationalNumber const*)a;
+ storm::RationalNumber* result_srn = new storm::RationalNumber(storm::utility::kwek_mehlhorn::sharpen<storm::RationalNumber, storm::RationalNumber>(precision, srn_a));
+ return (storm_rational_number_ptr)result_srn;
+}
+
int storm_rational_number_equal_modulo_precision(int relative, storm_rational_number_ptr a, storm_rational_number_ptr b, storm_rational_number_ptr precision) {
#ifndef RATIONAL_NUMBER_THREAD_SAFE
std::lock_guard<std::mutex> lock(rationalNumberMutex);
diff --git a/resources/3rdparty/sylvan/src/storm_wrapper.h b/resources/3rdparty/sylvan/src/storm_wrapper.h
index 2fb6bb214..d5ee75259 100644
--- a/resources/3rdparty/sylvan/src/storm_wrapper.h
+++ b/resources/3rdparty/sylvan/src/storm_wrapper.h
@@ -8,7 +8,7 @@
#ifdef __cplusplus
extern "C" {
#endif
-
+
/***************************************************
Function-wrappers for storm::RationalNumber
****************************************************/
@@ -49,6 +49,9 @@ extern "C" {
storm_rational_number_ptr storm_rational_number_floor(storm_rational_number_ptr a);
storm_rational_number_ptr storm_rational_number_ceil(storm_rational_number_ptr a);
+ storm_rational_number_ptr storm_double_sharpen(double value, size_t precision);
+ storm_rational_number_ptr storm_rational_number_sharpen(storm_rational_number_ptr a, size_t precision);
+
// Other operations.
int storm_rational_number_equal_modulo_precision(int relative, storm_rational_number_ptr a, storm_rational_number_ptr b, storm_rational_number_ptr precision);
diff --git a/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.c b/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.c
index 39de0400b..67919ed5f 100644
--- a/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.c
+++ b/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.c
@@ -590,6 +590,37 @@ TASK_IMPL_2(MTBDD, mtbdd_op_complement, MTBDD, a, size_t, k)
(void)k; // unused variable
}
+#include "sylvan_storm_rational_number.h"
+
+TASK_IMPL_2(MTBDD, mtbdd_op_sharpen, MTBDD, a, size_t, p)
+{
+ /* We only expect double or rational number terminals, or false */
+ if (a == mtbdd_false) return mtbdd_false;
+ if (a == mtbdd_true) return mtbdd_true;
+
+ // a != constant
+ mtbddnode_t na = MTBDD_GETNODE(a);
+
+ if (mtbddnode_isleaf(na)) {
+ if (mtbddnode_gettype(na) == 1) {
+ MTBDD result = mtbdd_storm_rational_number(storm_double_sharpen(mtbdd_getdouble(a), p));
+ return result;
+ } else if (mtbddnode_gettype(na) == srn_type) {
+ return mtbdd_storm_rational_number(storm_rational_number_sharpen((storm_rational_number_ptr)mtbdd_getstorm_rational_function_ptr(a), p));
+ } else {
+ printf("ERROR: Unsupported value type in sharpen.\n");
+ assert(0);
+ }
+ }
+
+ return mtbdd_invalid;
+}
+
+TASK_IMPL_2(MTBDD, mtbdd_sharpen, MTBDD, dd, size_t, p)
+{
+ return mtbdd_uapply(dd, TASK(mtbdd_op_sharpen), p);
+}
+
TASK_IMPL_3(BDD, mtbdd_min_abstract_representative, MTBDD, a, BDD, v, BDDVAR, prev_level) {
/* Maybe perform garbage collection */
sylvan_gc_test();
diff --git a/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h b/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h
index 27066fe19..5fcd69383 100644
--- a/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h
+++ b/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h
@@ -114,7 +114,7 @@ TASK_DECL_1(MTBDD, mtbdd_bool_to_int64, MTBDD)
*/
TASK_DECL_2(double, mtbdd_non_zero_count, MTBDD, size_t)
#define mtbdd_non_zero_count(dd, nvars) CALL(mtbdd_non_zero_count, dd, nvars)
-
+
// Checks whether the given MTBDD (does represents a zero leaf.
int mtbdd_iszero(MTBDD);
int mtbdd_isnonzero(MTBDD);
@@ -127,6 +127,13 @@ int mtbdd_isnonzero(MTBDD);
/* Create a MTBDD representing just <var> or the negation of <var> */
MTBDD mtbdd_ithvar(uint32_t var);
+/**
+ * Monad that sharpens the provided vector using the Kwek-Mehlhorn algorithm.
+ */
+TASK_DECL_2(MTBDD, mtbdd_op_sharpen, MTBDD, size_t)
+TASK_DECL_2(MTBDD, mtbdd_sharpen, MTBDD, size_t)
+#define mtbdd_sharpen(dd, p) CALL(mtbdd_sharpen, dd, p)
+
/**
* Unary operation Complement.
* Supported domains: Integer, Real
diff --git a/resources/3rdparty/sylvan/src/sylvan_obj_mtbdd_storm.hpp b/resources/3rdparty/sylvan/src/sylvan_obj_mtbdd_storm.hpp
index 23ae08b1e..cab9be01b 100644
--- a/resources/3rdparty/sylvan/src/sylvan_obj_mtbdd_storm.hpp
+++ b/resources/3rdparty/sylvan/src/sylvan_obj_mtbdd_storm.hpp
@@ -30,6 +30,8 @@ Mtbdd Maximum() const;
bool EqualNorm(const Mtbdd& other, double epsilon) const;
bool EqualNormRel(const Mtbdd& other, double epsilon) const;
+Mtbdd SharpenKwekMehlhorn(size_t precision) const;
+
// Functions that operate on Mtbdds over rational numbers.
static Mtbdd stormRationalNumberTerminal(storm::RationalNumber const& value);
diff --git a/resources/3rdparty/sylvan/src/sylvan_obj_storm.cpp b/resources/3rdparty/sylvan/src/sylvan_obj_storm.cpp
index 5cbf27cc1..dcddb3b6f 100644
--- a/resources/3rdparty/sylvan/src/sylvan_obj_storm.cpp
+++ b/resources/3rdparty/sylvan/src/sylvan_obj_storm.cpp
@@ -197,6 +197,12 @@ Mtbdd::EqualNormRel(const Mtbdd& other, double epsilon) const {
return mtbdd_equal_norm_rel_d(mtbdd, other.mtbdd, epsilon);
}
+Mtbdd
+Mtbdd::SharpenKwekMehlhorn(size_t precision) const {
+ LACE_ME;
+ return mtbdd_sharpen(mtbdd, precision);
+}
+
// Functions for Mtbdds over rational numbers.
Mtbdd
Mtbdd::stormRationalNumberTerminal(storm::RationalNumber const& value)
diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
index e849f7b9e..faa04d400 100644
--- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
+++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
@@ -14,6 +14,7 @@
#include "storm/exceptions/InvalidStateException.h"
#include "storm/exceptions/UnmetRequirementException.h"
#include "storm/exceptions/NotSupportedException.h"
+#include "storm/exceptions/PrecisionExceededException.h"
namespace storm {
namespace solver {
diff --git a/src/storm/storage/dd/Add.cpp b/src/storm/storage/dd/Add.cpp
index a7129ced8..92a0f78ab 100644
--- a/src/storm/storage/dd/Add.cpp
+++ b/src/storm/storage/dd/Add.cpp
@@ -1,5 +1,7 @@
#include "storm/storage/dd/Add.h"
+#include <cstdint>
+
#include <boost/algorithm/string/join.hpp>
#include "storm/storage/dd/DdMetaVariable.h"
@@ -12,6 +14,7 @@
#include "storm/utility/macros.h"
#include "storm/exceptions/InvalidArgumentException.h"
#include "storm/exceptions/InvalidOperationException.h"
+#include "storm/exceptions/NotSupportedException.h"
#include "storm-config.h"
#include "storm/adapters/RationalFunctionAdapter.h"
@@ -141,6 +144,21 @@ namespace storm {
return Add<LibraryType, ValueType>(this->getDdManager(), internalAdd.ceil(), this->getContainedMetaVariables());
}
+ template<DdType LibraryType, typename ValueType>
+ Add<LibraryType, storm::RationalNumber> Add<LibraryType, ValueType>::sharpenKwekMehlhorn(uint64_t precision) const {
+ STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Operation not supported.");
+ }
+
+ template<>
+ Add<storm::dd::DdType::Sylvan, storm::RationalNumber> Add<storm::dd::DdType::Sylvan, double>::sharpenKwekMehlhorn(uint64_t precision) const {
+ return Add<storm::dd::DdType::Sylvan, storm::RationalNumber>(this->getDdManager(), internalAdd.sharpenKwekMehlhorn(static_cast<std::size_t>(precision)), this->getContainedMetaVariables());
+ }
+
+ template<>
+ Add<storm::dd::DdType::Sylvan, storm::RationalNumber> Add<storm::dd::DdType::Sylvan, storm::RationalNumber>::sharpenKwekMehlhorn(uint64_t precision) const {
+ return Add<storm::dd::DdType::Sylvan, storm::RationalNumber>(this->getDdManager(), internalAdd.sharpenKwekMehlhorn(static_cast<std::size_t>(precision)), this->getContainedMetaVariables());
+ }
+
template<DdType LibraryType, typename ValueType>
Add<LibraryType, ValueType> Add<LibraryType, ValueType>::minimum(Add<LibraryType, ValueType> const& other) const {
return Add<LibraryType, ValueType>(this->getDdManager(), internalAdd.minimum(other), Dd<LibraryType>::joinMetaVariables(*this, other));
diff --git a/src/storm/storage/dd/Add.h b/src/storm/storage/dd/Add.h
index bbcfb099b..1094bf0a1 100644
--- a/src/storm/storage/dd/Add.h
+++ b/src/storm/storage/dd/Add.h
@@ -237,6 +237,13 @@ namespace storm {
*/
Add<LibraryType, ValueType> ceil() const;
+ /*!
+ * Retrieves the function that sharpens all values in the current ADD with the Kwek-Mehlhorn algorithm.
+ *
+ * @return The resulting ADD.
+ */
+ Add<LibraryType, storm::RationalNumber> sharpenKwekMehlhorn(uint64_t precision) const;
+
/*!
* Retrieves the function that maps all evaluations to the minimum of the function values of the two ADDs.
*
diff --git a/src/storm/storage/dd/cudd/InternalCuddAdd.cpp b/src/storm/storage/dd/cudd/InternalCuddAdd.cpp
index 187968f23..f129177cb 100644
--- a/src/storm/storage/dd/cudd/InternalCuddAdd.cpp
+++ b/src/storm/storage/dd/cudd/InternalCuddAdd.cpp
@@ -10,6 +10,7 @@
#include "storm/utility/constants.h"
#include "storm/utility/macros.h"
#include "storm/exceptions/NotImplementedException.h"
+#include "storm/exceptions/NotSupportedException.h"
namespace storm {
namespace dd {
@@ -127,6 +128,11 @@ namespace storm {
return InternalAdd<DdType::CUDD, ValueType>(ddManager, this->getCuddAdd().Ceil());
}
+ template<typename ValueType>
+ InternalAdd<DdType::CUDD, storm::RationalNumber> InternalAdd<DdType::CUDD, ValueType>::sharpenKwekMehlhorn(size_t precision) const {
+ STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Operation not supported");
+ }
+
template<typename ValueType>
InternalAdd<DdType::CUDD, ValueType> InternalAdd<DdType::CUDD, ValueType>::minimum(InternalAdd<DdType::CUDD, ValueType> const& other) const {
return InternalAdd<DdType::CUDD, ValueType>(ddManager, this->getCuddAdd().Minimum(other.getCuddAdd()));
diff --git a/src/storm/storage/dd/cudd/InternalCuddAdd.h b/src/storm/storage/dd/cudd/InternalCuddAdd.h
index cc619cd11..49adeea69 100644
--- a/src/storm/storage/dd/cudd/InternalCuddAdd.h
+++ b/src/storm/storage/dd/cudd/InternalCuddAdd.h
@@ -6,6 +6,8 @@
#include <functional>
#include <memory>
+#include "storm/adapters/RationalNumberAdapter.h"
+
#include "storm/storage/dd/DdType.h"
#include "storm/storage/dd/InternalAdd.h"
#include "storm/storage/dd/Odd.h"
@@ -237,6 +239,13 @@ namespace storm {
*/
InternalAdd<DdType::CUDD, ValueType> ceil() const;
+ /*!
+ * Retrieves the function that sharpens all values in the current ADD with the Kwek-Mehlhorn algorithm.
+ *
+ * @return The resulting ADD.
+ */
+ InternalAdd<DdType::CUDD, storm::RationalNumber> sharpenKwekMehlhorn(size_t precision) const;
+
/*!
* Retrieves the function that maps all evaluations to the minimum of the function values of the two ADDs.
*
diff --git a/src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp b/src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp
index 46167f116..d3f0f3b85 100644
--- a/src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp
+++ b/src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp
@@ -10,6 +10,7 @@
#include "storm/utility/constants.h"
#include "storm/exceptions/NotImplementedException.h"
#include "storm/exceptions/InvalidOperationException.h"
+#include "storm/exceptions/NotSupportedException.h"
#include "storm-config.h"
@@ -394,6 +395,18 @@ namespace storm {
return InternalAdd<DdType::Sylvan, storm::RationalFunction>(ddManager, this->sylvanMtbdd.CeilRF());
}
#endif
+
+ template<typename ValueType>
+ InternalAdd<DdType::Sylvan, storm::RationalNumber> InternalAdd<DdType::Sylvan, ValueType>::sharpenKwekMehlhorn(size_t precision) const {
+ return InternalAdd<DdType::Sylvan, storm::RationalNumber>(ddManager, this->sylvanMtbdd.SharpenKwekMehlhorn(precision));
+ }
+
+#ifdef STORM_HAVE_CARL
+ template<>
+ InternalAdd<DdType::Sylvan, storm::RationalNumber> InternalAdd<DdType::Sylvan, storm::RationalFunction>::sharpenKwekMehlhorn(size_t precision) const {
+ STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Operation not supported.");
+ }
+#endif
template<typename ValueType>
InternalAdd<DdType::Sylvan, ValueType> InternalAdd<DdType::Sylvan, ValueType>::minimum(InternalAdd<DdType::Sylvan, ValueType> const& other) const {
diff --git a/src/storm/storage/dd/sylvan/InternalSylvanAdd.h b/src/storm/storage/dd/sylvan/InternalSylvanAdd.h
index 201048b6f..99f7609ef 100644
--- a/src/storm/storage/dd/sylvan/InternalSylvanAdd.h
+++ b/src/storm/storage/dd/sylvan/InternalSylvanAdd.h
@@ -197,7 +197,7 @@ namespace storm {
* Retrieves the function that represents the current ADD to the power of the given ADD.
*
* @other The exponent function (given as an ADD).
- * @retur The resulting ADD.
+ * @return The resulting ADD.
*/
InternalAdd<DdType::Sylvan, ValueType> pow(InternalAdd<DdType::Sylvan, ValueType> const& other) const;
@@ -205,7 +205,7 @@ namespace storm {
* Retrieves the function that represents the current ADD modulo the given ADD.
*
* @other The modul function (given as an ADD).
- * @retur The resulting ADD.
+ * @return The resulting ADD.
*/
InternalAdd<DdType::Sylvan, ValueType> mod(InternalAdd<DdType::Sylvan, ValueType> const& other) const;
@@ -214,24 +214,31 @@ namespace storm {
* ADD.
*
* @other The base function (given as an ADD).
- * @retur The resulting ADD.
+ * @return The resulting ADD.
*/
InternalAdd<DdType::Sylvan, ValueType> logxy(InternalAdd<DdType::Sylvan, ValueType> const& other) const;
/*!
* Retrieves the function that floors all values in the current ADD.
*
- * @retur The resulting ADD.
+ * @return The resulting ADD.
*/
InternalAdd<DdType::Sylvan, ValueType> floor() const;
/*!
* Retrieves the function that ceils all values in the current ADD.
*
- * @retur The resulting ADD.
+ * @return The resulting ADD.
*/
InternalAdd<DdType::Sylvan, ValueType> ceil() const;
+ /*!
+ * Retrieves the function that sharpens all values in the current ADD with the Kwek-Mehlhorn algorithm.
+ *
+ * @return The resulting ADD.
+ */
+ InternalAdd<DdType::Sylvan, storm::RationalNumber> sharpenKwekMehlhorn(size_t precision) const;
+
/*!
* Retrieves the function that maps all evaluations to the minimum of the function values of the two ADDs.
*
diff --git a/src/storm/utility/KwekMehlhorn.cpp b/src/storm/utility/KwekMehlhorn.cpp
new file mode 100644
index 000000000..9c76be068
--- /dev/null
+++ b/src/storm/utility/KwekMehlhorn.cpp
@@ -0,0 +1,79 @@
+#include "storm/utility/KwekMehlhorn.h"
+
+#include "storm/adapters/RationalNumberAdapter.h"
+
+#include "storm/utility/constants.h"
+#include "storm/utility/macros.h"
+
+#include "storm/exceptions/PrecisionExceededException.h"
+
+namespace storm {
+ namespace utility{
+ namespace kwek_mehlhorn {
+
+ template<typename IntegerType>
+ std::pair<IntegerType, IntegerType> findRational(IntegerType const& alpha, IntegerType const& beta, IntegerType const& gamma, IntegerType const& delta) {
+ std::pair<IntegerType, IntegerType> alphaDivBetaPair = storm::utility::divide(alpha, beta);
+ std::pair<IntegerType, IntegerType> gammaDivDeltaPair = storm::utility::divide(gamma, delta);
+
+ if (alphaDivBetaPair.first == gammaDivDeltaPair.first && !storm::utility::isZero(alphaDivBetaPair.second)) {
+ std::pair<IntegerType, IntegerType> 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(alphaDivBetaPair.second) ? alphaDivBetaPair.first : alphaDivBetaPair.first + storm::utility::one<IntegerType>(), storm::utility::one<IntegerType>());
+ return result;
+ }
+ }
+
+ template<typename RationalType, typename ImpreciseType>
+ std::pair<typename NumberTraits<RationalType>::IntegerType, typename NumberTraits<RationalType>::IntegerType> truncateToRational(ImpreciseType const& value, uint64_t precision) {
+ typedef typename NumberTraits<RationalType>::IntegerType IntegerType;
+
+ IntegerType powerOfTen = storm::utility::pow(storm::utility::convertNumber<IntegerType>(10ull), precision);
+ IntegerType truncated = storm::utility::trunc<RationalType>(value * powerOfTen);
+ return std::make_pair(truncated, powerOfTen);
+ }
+
+ template<typename RationalType>
+ std::pair<typename NumberTraits<RationalType>::IntegerType, typename NumberTraits<RationalType>::IntegerType> truncateToRational(double const& value, uint64_t precision) {
+ STORM_LOG_THROW(precision < 17, storm::exceptions::PrecisionExceededException, "Exceeded precision of double, consider switching to rational numbers.");
+
+ double powerOfTen = std::pow(10, precision);
+ double truncated = storm::utility::trunc<double>(value * powerOfTen);
+ return std::make_pair(truncated, powerOfTen);
+ }
+
+ template<typename RationalType, typename ImpreciseType>
+ RationalType findRational(uint64_t precision, ImpreciseType const& value) {
+ typedef typename NumberTraits<RationalType>::IntegerType IntegerType;
+
+ std::pair<IntegerType, IntegerType> truncatedFraction = truncateToRational<RationalType>(value, precision);
+ std::pair<IntegerType, IntegerType> result = findRational<IntegerType>(truncatedFraction.first, truncatedFraction.second, truncatedFraction.first + storm::utility::one<IntegerType>(), truncatedFraction.second);
+
+ // Convert one of the arguments to a rational type to not get integer division.
+ return storm::utility::convertNumber<RationalType>(result.first) / result.second;
+ }
+
+ template<typename RationalType, typename ImpreciseType>
+ RationalType sharpen(uint64_t precision, ImpreciseType const& value) {
+ ImpreciseType integer = storm::utility::floor(value);
+ ImpreciseType fraction = value - integer;
+ auto rational = findRational<RationalType>(precision, fraction);
+ return storm::utility::convertNumber<RationalType>(integer) + rational;
+ }
+
+ template<typename RationalType, typename ImpreciseType>
+ void sharpen(uint64_t precision, std::vector<ImpreciseType> const& input, std::vector<RationalType>& output) {
+ for (uint64_t index = 0; index < input.size(); ++index) {
+ output[index] = sharpen<RationalType, ImpreciseType>(precision, input[index]);
+ }
+ }
+
+ template void sharpen(uint64_t precision, std::vector<double> const& input, std::vector<storm::RationalNumber>& output);
+ template void sharpen(uint64_t precision, std::vector<storm::RationalNumber> const& input, std::vector<storm::RationalNumber>& output);
+
+ }
+ }
+}
diff --git a/src/storm/utility/KwekMehlhorn.h b/src/storm/utility/KwekMehlhorn.h
index 72038475d..245c875a3 100644
--- a/src/storm/utility/KwekMehlhorn.h
+++ b/src/storm/utility/KwekMehlhorn.h
@@ -1,68 +1,32 @@
#pragma once
-#include "storm/utility/constants.h"
-#include "storm/utility/macros.h"
+#include <utility>
+#include <cstdint>
+#include <vector>
-#include "storm/exceptions/PrecisionExceededException.h"
+#include "storm/utility/NumberTraits.h"
namespace storm {
namespace utility{
namespace kwek_mehlhorn {
template<typename IntegerType>
- std::pair<IntegerType, IntegerType> findRational(IntegerType const& alpha, IntegerType const& beta, IntegerType const& gamma, IntegerType const& delta) {
- std::pair<IntegerType, IntegerType> alphaDivBetaPair = storm::utility::divide(alpha, beta);
- std::pair<IntegerType, IntegerType> gammaDivDeltaPair = storm::utility::divide(gamma, delta);
-
- if (alphaDivBetaPair.first == gammaDivDeltaPair.first && !storm::utility::isZero(alphaDivBetaPair.second)) {
- std::pair<IntegerType, IntegerType> 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(alphaDivBetaPair.second) ? alphaDivBetaPair.first : alphaDivBetaPair.first + storm::utility::one<IntegerType>(), storm::utility::one<IntegerType>());
- return result;
- }
- }
+ std::pair<IntegerType, IntegerType> findRational(IntegerType const& alpha, IntegerType const& beta, IntegerType const& gamma, IntegerType const& delta);
template<typename RationalType, typename ImpreciseType>
- std::pair<typename NumberTraits<RationalType>::IntegerType, typename NumberTraits<RationalType>::IntegerType> truncateToRational(ImpreciseType const& value, uint64_t precision) {
- typedef typename NumberTraits<RationalType>::IntegerType IntegerType;
-
- IntegerType powerOfTen = storm::utility::pow(storm::utility::convertNumber<IntegerType>(10ull), precision);
- IntegerType truncated = storm::utility::trunc<RationalType>(value * powerOfTen);
- return std::make_pair(truncated, powerOfTen);
- }
+ std::pair<typename NumberTraits<RationalType>::IntegerType, typename NumberTraits<RationalType>::IntegerType> truncateToRational(ImpreciseType const& value, uint64_t precision);
template<typename RationalType>
- std::pair<typename NumberTraits<RationalType>::IntegerType, typename NumberTraits<RationalType>::IntegerType> truncateToRational(double const& value, uint64_t precision) {
- STORM_LOG_THROW(precision < 17, storm::exceptions::PrecisionExceededException, "Exceeded precision of double, consider switching to rational numbers.");
-
- double powerOfTen = std::pow(10, precision);
- double truncated = storm::utility::trunc<double>(value * powerOfTen);
- return std::make_pair(truncated, powerOfTen);
- }
+ std::pair<typename NumberTraits<RationalType>::IntegerType, typename NumberTraits<RationalType>::IntegerType> truncateToRational(double const& value, uint64_t precision);
template<typename RationalType, typename ImpreciseType>
- RationalType findRational(uint64_t precision, ImpreciseType const& value) {
- typedef typename NumberTraits<RationalType>::IntegerType IntegerType;
-
- std::pair<IntegerType, IntegerType> truncatedFraction = truncateToRational<RationalType>(value, precision);
- std::pair<IntegerType, IntegerType> result = findRational<IntegerType>(truncatedFraction.first, truncatedFraction.second, truncatedFraction.first + storm::utility::one<IntegerType>(), truncatedFraction.second);
-
- // Convert one of the arguments to a rational type to not get integer division.
- return storm::utility::convertNumber<RationalType>(result.first) / result.second;
- }
+ RationalType findRational(uint64_t precision, ImpreciseType const& value);
+
+ template<typename RationalType, typename ImpreciseType>
+ RationalType sharpen(uint64_t precision, ImpreciseType const& value);
template<typename RationalType, typename ImpreciseType>
- void sharpen(uint64_t precision, std::vector<ImpreciseType> const& input, std::vector<RationalType>& output) {
- for (uint64_t index = 0; index < input.size(); ++index) {
- ImpreciseType integer = storm::utility::floor(input[index]);
- ImpreciseType fraction = input[index] - integer;
- auto rational = findRational<RationalType>(precision, fraction);
- output[index] = storm::utility::convertNumber<RationalType>(integer) + rational;
- }
- }
+ void sharpen(uint64_t precision, std::vector<ImpreciseType> const& input, std::vector<RationalType>& output);
}
}
diff --git a/src/test/storm/storage/SylvanDdTest.cpp b/src/test/storm/storage/SylvanDdTest.cpp
index 4b8ec2f38..ca294b273 100644
--- a/src/test/storm/storage/SylvanDdTest.cpp
+++ b/src/test/storm/storage/SylvanDdTest.cpp
@@ -835,6 +835,26 @@ TEST(SylvanDd, AddOddTest) {
EXPECT_EQ(106ul, matrix.getNonzeroEntryCount());
}
+TEST(SylvanDd, AddSharpenTest) {
+ std::shared_ptr<storm::dd::DdManager<storm::dd::DdType::Sylvan>> manager(new storm::dd::DdManager<storm::dd::DdType::Sylvan>());
+ std::pair<storm::expressions::Variable, storm::expressions::Variable> x = manager->addMetaVariable("x", 1, 9);
+
+ storm::dd::Add<storm::dd::DdType::Sylvan, double> dd = manager->template getAddOne<double>();
+ ASSERT_NO_THROW(dd.setValue(x.first, 4, 1.89999999));
+ ASSERT_EQ(2ul, dd.getLeafCount());
+
+ storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalNumber> sharpened = dd.sharpenKwekMehlhorn(1);
+
+ std::map<storm::expressions::Variable, int_fast64_t> metaVariableToValueMap;
+ metaVariableToValueMap.emplace(x.first, 4);
+
+ sharpened = dd.sharpenKwekMehlhorn(1);
+ ASSERT_EQ(storm::utility::convertNumber<storm::RationalNumber>(std::string("9/5")), sharpened.getValue(metaVariableToValueMap));
+
+ sharpened = dd.sharpenKwekMehlhorn(2);
+ ASSERT_EQ(storm::utility::convertNumber<storm::RationalNumber>(std::string("19/10")), sharpened.getValue(metaVariableToValueMap));
+}
+
TEST(SylvanDd, BddOddTest) {
std::shared_ptr<storm::dd::DdManager<storm::dd::DdType::Sylvan>> manager(new storm::dd::DdManager<storm::dd::DdType::Sylvan>());
std::pair<storm::expressions::Variable, storm::expressions::Variable> a = manager->addMetaVariable("a");