#pragma once #include "storm/utility/constants.h" #include "storm/utility/macros.h" #include "storm/exceptions/PrecisionExceededException.h" namespace storm { namespace utility{ namespace kwek_mehlhorn { 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); 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); return result; } else { auto result = std::make_pair(storm::utility::isZero(alphaModBeta) ? alphaDivBetaFloor : alphaDivBetaFloor + storm::utility::one(), storm::utility::one()); return result; } } template std::pair::IntegerType, typename NumberTraits::IntegerType> truncateToRational(ImpreciseType const& value, uint64_t precision) { typedef typename NumberTraits::IntegerType IntegerType; IntegerType powerOfTen = storm::utility::pow(storm::utility::convertNumber(10ull), precision); IntegerType truncated = storm::utility::trunc(value * powerOfTen); return std::make_pair(truncated, powerOfTen); } template std::pair::IntegerType, typename NumberTraits::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(value * powerOfTen); return std::make_pair(truncated, powerOfTen); } template RationalType findRational(uint64_t precision, ImpreciseType const& value) { typedef typename NumberTraits::IntegerType IntegerType; std::pair truncatedFraction = truncateToRational(value, precision); std::pair result = findRational(truncatedFraction.first, truncatedFraction.second, truncatedFraction.first + storm::utility::one(), truncatedFraction.second); // Convert one of the arguments to a rational type to not get integer division. return storm::utility::convertNumber(result.first) / result.second; } template void sharpen(uint64_t precision, std::vector const& input, std::vector& output) { for (uint64_t index = 0; index < input.size(); ++index) { ImpreciseType integer = storm::utility::floor(input[index]); ImpreciseType fraction = input[index] - integer; output[index] = storm::utility::convertNumber(integer) + findRational(precision, fraction); } } } } }