started on symbolic rational search

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);

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);
     

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();

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

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);
 

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)

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 {

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));

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.
              *

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()));

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.
              *

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 {

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.
              *

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);
+            
+        }
+    }
+}

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);
             
         }
     }

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");