diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9fafa9849..5e971c46e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,6 +11,7 @@ Version 1.1.x
 - parametric model checking has an own binary
 - solvers can now expose requirements
 - unbounded reachability and reachability rewards now correctly respect solver requirements
+- sound (interval) value iteration
 
 ### Version 1.1.1
 - c++ api changes: Building model takes BuilderOptions instead of extended list of Booleans, does not depend on settings anymore.
diff --git a/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp b/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp
index 6c6b43be2..ceb04ef40 100644
--- a/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp
+++ b/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp
@@ -21,7 +21,7 @@ namespace storm {
             const std::string MinMaxEquationSolverSettings::valueIterationMultiplicationStyleOptionName = "vimult";
 
             MinMaxEquationSolverSettings::MinMaxEquationSolverSettings() : ModuleSettings(moduleName) {
-                std::vector<std::string> minMaxSolvingTechniques = {"vi", "value-iteration", "pi", "policy-iteration", "linear-programming", "lp", "acyclic"};
+                std::vector<std::string> minMaxSolvingTechniques = {"vi", "value-iteration", "pi", "policy-iteration", "linear-programming", "lp", "acyclic", "ratsearch"};
                 this->addOption(storm::settings::OptionBuilder(moduleName, solvingMethodOptionName, false, "Sets which min/max linear equation solving technique is preferred.")
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a min/max linear equation solving technique.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(minMaxSolvingTechniques)).setDefaultValueString("vi").build()).build());
                 
@@ -50,6 +50,8 @@ namespace storm {
                     return storm::solver::MinMaxMethod::LinearProgramming;
                 } else if (minMaxEquationSolvingTechnique == "acyclic") {
                     return storm::solver::MinMaxMethod::Acyclic;
+                } else if (minMaxEquationSolvingTechnique == "ratsearch") {
+                    return storm::solver::MinMaxMethod::RationalSearch;
                 }
                 STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown min/max equation solving technique '" << minMaxEquationSolvingTechnique << "'.");
             }
diff --git a/src/storm/solver/AbstractEquationSolver.cpp b/src/storm/solver/AbstractEquationSolver.cpp
index 18a82f498..900ebaf91 100644
--- a/src/storm/solver/AbstractEquationSolver.cpp
+++ b/src/storm/solver/AbstractEquationSolver.cpp
@@ -40,6 +40,11 @@ namespace storm {
             return *terminationCondition;
         }
         
+        template<typename ValueType>
+        std::unique_ptr<TerminationCondition<ValueType>> const& AbstractEquationSolver<ValueType>::getTerminationConditionPointer() const {
+            return terminationCondition;
+        }
+        
         template<typename ValueType>
         bool AbstractEquationSolver<ValueType>::terminateNow(std::vector<ValueType> const& values, SolverGuarantee const& guarantee) const {
             if (!this->hasCustomTerminationCondition()) {
@@ -198,10 +203,10 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void AbstractEquationSolver<ValueType>::startMeasureProgress() const {
+        void AbstractEquationSolver<ValueType>::startMeasureProgress(uint64_t startingIteration) const {
             timeOfStart = std::chrono::high_resolution_clock::now();
             timeOfLastMessage = timeOfStart;
-            iterationOfLastMessage = 0;
+            iterationOfLastMessage = startingIteration;
         }
         
         template<typename ValueType>
diff --git a/src/storm/solver/AbstractEquationSolver.h b/src/storm/solver/AbstractEquationSolver.h
index 362239360..4e27517fa 100644
--- a/src/storm/solver/AbstractEquationSolver.h
+++ b/src/storm/solver/AbstractEquationSolver.h
@@ -143,7 +143,7 @@ namespace storm {
             /*!
              * Starts to measure progress.
              */
-            void startMeasureProgress() const;
+            void startMeasureProgress(uint64_t startingIteration = 0) const;
             
             /*!
              * Shows progress if this solver is asked to do so.
@@ -157,6 +157,7 @@ namespace storm {
              * @return The custom termination condition.
              */
             TerminationCondition<ValueType> const& getTerminationCondition() const;
+            std::unique_ptr<TerminationCondition<ValueType>> const& getTerminationConditionPointer() const;
             
             void createUpperBoundsVector(std::unique_ptr<std::vector<ValueType>>& upperBoundsVector, uint64_t length) const;
             void createLowerBoundsVector(std::vector<ValueType>& lowerBoundsVector) const;
diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
index 0ddd5b5fe..27577d3ba 100644
--- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
+++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
@@ -1,5 +1,7 @@
 #include "storm/solver/IterativeMinMaxLinearEquationSolver.h"
 
+#include "storm/utility/ConstantsComparator.h"
+
 #include "storm/settings/SettingsManager.h"
 #include "storm/settings/modules/GeneralSettings.h"
 #include "storm/settings/modules/MinMaxEquationSolverSettings.h"
@@ -9,6 +11,7 @@
 #include "storm/exceptions/InvalidSettingsException.h"
 #include "storm/exceptions/InvalidStateException.h"
 #include "storm/exceptions/UnmetRequirementException.h"
+#include "storm/exceptions/NotSupportedException.h"
 
 namespace storm {
     namespace solver {
@@ -40,6 +43,7 @@ namespace storm {
                 case MinMaxMethod::ValueIteration: this->solutionMethod = SolutionMethod::ValueIteration; break;
                 case MinMaxMethod::PolicyIteration: this->solutionMethod = SolutionMethod::PolicyIteration; break;
                 case MinMaxMethod::Acyclic: this->solutionMethod = SolutionMethod::Acyclic; break;
+                case MinMaxMethod::RationalSearch: this->solutionMethod = SolutionMethod::RationalSearch; break;
                 default:
                     STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique for iterative MinMax linear equation solver.");
             }
@@ -128,6 +132,8 @@ namespace storm {
                     return solveEquationsPolicyIteration(dir, x, b);
                 case IterativeMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::Acyclic:
                     return solveEquationsAcyclic(dir, x, b);
+                case IterativeMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::RationalSearch:
+                    return solveEquationsRationalSearch(dir, x, b);
                 default:
                     STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "This solver does not implement the selected solution method");
             }
@@ -287,6 +293,46 @@ namespace storm {
             return requirements;
         }
 
+        template<typename ValueType>
+        typename IterativeMinMaxLinearEquationSolver<ValueType>::ValueIterationResult IterativeMinMaxLinearEquationSolver<ValueType>::performValueIteration(OptimizationDirection dir, std::vector<ValueType>*& currentX, std::vector<ValueType>*& newX, std::vector<ValueType> const& b, ValueType const& precision, bool relative, SolverGuarantee const& guarantee, uint64_t currentIterations) const {
+            
+            // Get handle to linear equation solver.
+            storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver = *this->linEqSolverA;
+            
+            // Allow aliased multiplications.
+            bool useGaussSeidelMultiplication = linearEquationSolver.supportsGaussSeidelMultiplication() && settings.getValueIterationMultiplicationStyle() == storm::solver::MultiplicationStyle::GaussSeidel;
+            
+            // Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations.
+            uint64_t iterations = currentIterations;
+            
+            Status status = Status::InProgress;
+            while (status == Status::InProgress) {
+                // Compute x' = min/max(A*x + b).
+                if (useGaussSeidelMultiplication) {
+                    // Copy over the current vector so we can modify it in-place.
+                    *newX = *currentX;
+                    linearEquationSolver.multiplyAndReduceGaussSeidel(dir, this->A->getRowGroupIndices(), *newX, &b);
+                } else {
+                    linearEquationSolver.multiplyAndReduce(dir, this->A->getRowGroupIndices(), *currentX, &b, *newX);
+                }
+                
+                // Determine whether the method converged.
+                if (storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, precision, relative)) {
+                    status = Status::Converged;
+                }
+                
+                // Potentially show progress.
+                this->showProgressIterative(iterations);
+                
+                // Update environment variables.
+                std::swap(currentX, newX);
+                ++iterations;
+                status = updateStatusIfNotConverged(status, *currentX, iterations, guarantee);
+            }
+            
+            return ValueIterationResult(iterations, status);
+        }
+        
         template<typename ValueType>
         bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             if (!this->linEqSolverA) {
@@ -327,44 +373,18 @@ namespace storm {
                 if (dir == storm::OptimizationDirection::Minimize) {
                     guarantee = SolverGuarantee::GreaterOrEqual;
                 }
+            } else {
+                // If no initial scheduler is given, we start from the lower bound.
+                this->createLowerBoundsVector(x);
             }
 
-            // Allow aliased multiplications.
-            bool useGaussSeidelMultiplication = this->linEqSolverA->supportsGaussSeidelMultiplication() && settings.getValueIterationMultiplicationStyle() == storm::solver::MultiplicationStyle::GaussSeidel;
-            
             std::vector<ValueType>* newX = auxiliaryRowGroupVector.get();
             std::vector<ValueType>* currentX = &x;
             
-            // Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations.
-            uint64_t iterations = 0;
-
             this->startMeasureProgress();
-            Status status = Status::InProgress;
-            while (status == Status::InProgress) {
-                // Compute x' = min/max(A*x + b).
-                if (useGaussSeidelMultiplication) {
-                    // Copy over the current vector so we can modify it in-place.
-                    *newX = *currentX;
-                    this->linEqSolverA->multiplyAndReduceGaussSeidel(dir, this->A->getRowGroupIndices(), *newX, &b);
-                } else {
-                    this->linEqSolverA->multiplyAndReduce(dir, this->A->getRowGroupIndices(), *currentX, &b, *newX);
-                }
-                
-                // Determine whether the method converged.
-                if (storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion())) {
-                    status = Status::Converged;
-                }
-                
-                // Potentially show progress.
-                this->showProgressIterative(iterations);
-                
-                // Update environment variables.
-                std::swap(currentX, newX);
-                ++iterations;
-                status = updateStatusIfNotConverged(status, *currentX, iterations, guarantee);
-            }
-            
-            reportStatus(status, iterations);
+            ValueIterationResult result = performValueIteration(dir, currentX, newX, b, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion(), guarantee, 0);
+
+            reportStatus(result.status, result.iterations);
             
             // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
             // is currently stored in currentX, but x is the output vector.
@@ -377,12 +397,12 @@ namespace storm {
                 this->schedulerChoices = std::vector<uint_fast64_t>(this->A->getRowGroupCount());
                 this->linEqSolverA->multiplyAndReduce(dir, this->A->getRowGroupIndices(), x, &b, *currentX, &this->schedulerChoices.get());
             }
-
+            
             if (!this->isCachingEnabled()) {
                 clearCache();
             }
             
-            return status == Status::Converged || status == Status::TerminatedEarly;
+            return result.status == Status::Converged || result.status == Status::TerminatedEarly;
         }
         
         template<typename ValueType>
@@ -409,6 +429,12 @@ namespace storm {
             return result;
         }
         
+        /*!
+         * This version of value iteration is sound, because it approaches the solution from below and above. This
+         * technique is due to Haddad and Monmege (Interval iteration algorithm for MDPs and IMDPs, TCS 2017) and was
+         * extended to rewards by Baier, Klein, Leuschner, Parker and Wunderlich (Ensuring the Reliability of Your
+         * Model Checker: Interval Iteration for Markov Decision Processes, CAV 2017).
+         */
         template<typename ValueType>
         bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsSoundValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             STORM_LOG_THROW(this->hasUpperBound(), storm::exceptions::UnmetRequirementException, "Solver requires upper bound, but none was given.");
@@ -580,6 +606,225 @@ namespace storm {
             return status == Status::Converged;
         }
         
+        template<typename ValueType>
+        bool IterativeMinMaxLinearEquationSolver<ValueType>::isSolution(storm::OptimizationDirection dir, storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType> const& values, std::vector<ValueType> const& b) {
+            storm::utility::ConstantsComparator<ValueType> comparator;
+            
+            auto valueIt = values.begin();
+            auto bIt = b.begin();
+            for (uint64_t group = 0; group < matrix.getRowGroupCount(); ++group, ++valueIt) {
+                ValueType groupValue = *bIt;
+                uint64_t row = matrix.getRowGroupIndices()[group];
+                groupValue += matrix.multiplyRowWithVector(row, values);
+
+                ++row;
+                ++bIt;
+
+                for (auto endRow = matrix.getRowGroupIndices()[group + 1]; row < endRow; ++row, ++bIt) {
+                    ValueType newValue = *bIt;
+                    newValue += matrix.multiplyRowWithVector(row, values);
+                    
+                    if ((dir == storm::OptimizationDirection::Minimize && newValue < groupValue) || (dir == storm::OptimizationDirection::Maximize && newValue > groupValue)) {
+                        groupValue = newValue;
+                    }
+                }
+                
+                // If the value does not match the one in the values vector, the given vector is not a solution.
+                if (!comparator.isEqual(groupValue, *valueIt)) {
+                    return false;
+                }
+            }
+            
+            // Checked all values at this point.
+            return true;
+        }
+
+        template<typename ValueType>
+        std::pair<ValueType, ValueType> findRational(ValueType const& alpha, ValueType const& beta, ValueType const& gamma, ValueType const& delta) {
+            ValueType alphaDivBeta = alpha / beta;
+            ValueType alphaDivBetaFloor = storm::utility::floor(alphaDivBeta);
+            ValueType gammaDivDeltaFloor = storm::utility::floor<ValueType>(gamma / delta);
+            
+            if (alphaDivBetaFloor == gammaDivDeltaFloor && !storm::utility::isInteger(alphaDivBeta)) {
+                std::pair<ValueType, ValueType> subresult = findRational(delta, storm::utility::mod(gamma, delta), beta, storm::utility::mod(alpha, beta));
+                auto result = std::make_pair(alphaDivBetaFloor * subresult.first + subresult.second, subresult.first);
+                
+                return result;
+            } else {
+                auto result = std::make_pair(storm::utility::ceil(alphaDivBeta), storm::utility::one<ValueType>());
+                return result;
+            }
+        }
+        
+        template<typename ValueType>
+        ValueType truncateToRational(double const& value, uint64_t precision) {
+            STORM_LOG_ASSERT(precision < 16, "Truncating via double is not sound.");
+            auto powerOfTen = std::pow(10, precision);
+            double truncated = std::trunc(value * powerOfTen);
+            return storm::utility::convertNumber<ValueType>(truncated) / storm::utility::convertNumber<ValueType>(powerOfTen);
+        }
+        
+        template<typename ValueType>
+        ValueType findRational(uint64_t precision, double const& value) {
+            ValueType truncatedFraction = truncateToRational<ValueType>(value, precision);
+
+            ValueType ten = storm::utility::convertNumber<ValueType>(10);
+            ValueType powerOfTen = storm::utility::pow<ValueType>(ten, precision);
+            ValueType multiplier = powerOfTen / storm::utility::denominator(truncatedFraction);
+            ValueType numerator = storm::utility::numerator(truncatedFraction) * multiplier;
+            
+            std::pair<ValueType, ValueType> result = findRational<ValueType>(numerator, powerOfTen, numerator + storm::utility::one<ValueType>(), powerOfTen);
+            return result.first / result.second;
+        }
+
+        template<typename ValueType>
+        bool IterativeMinMaxLinearEquationSolver<ValueType>::sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix<storm::RationalNumber> const& A, std::vector<double> const& x, std::vector<storm::RationalNumber> const& b, std::vector<storm::RationalNumber>& tmp) {
+            for (uint64_t p = 1; p < precision; ++p) {
+                for (uint64_t state = 0; state < x.size(); ++state) {
+                    storm::RationalNumber rationalX = storm::utility::convertNumber<storm::RationalNumber, double>(x[state]);
+                    storm::RationalNumber integer = storm::utility::floor(rationalX);
+                    storm::RationalNumber fraction = rationalX - integer;
+                    
+                    tmp[state] = integer + findRational<storm::RationalNumber>(p, storm::utility::convertNumber<double, storm::RationalNumber>(fraction));
+                }
+                if (IterativeMinMaxLinearEquationSolver<storm::RationalNumber>::isSolution(dir, A, tmp, b)) {
+                    return true;
+                }
+            }
+            return false;
+        }
+
+        template<typename ValueType>
+        void IterativeMinMaxLinearEquationSolver<ValueType>::createLinearEquationSolver() const {
+            this->linEqSolverA = this->linearEquationSolverFactory->create(*this->A);
+        }
+
+        template<typename ValueType>
+        template<typename ImpreciseValueType>
+        std::enable_if<std::is_same<ValueType, ImpreciseValueType>::value, bool> IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+
+        }
+
+        template<typename ValueType>
+        template<typename ImpreciseValueType>
+        std::enable_if<!std::is_same<ValueType, ImpreciseValueType>::value, bool> IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+            
+        }
+
+        template<typename ValueType>
+        bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearch(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+            if (this->getSettings().getForceSoundness()) {
+                solveEquationsRationalSearchHelper<ValueType>(dir, x, b);
+            } else {
+                solveEquationsRationalSearchHelper<double>(dir, x, b);
+            }
+            
+            // Create an imprecise solver with caching enabled.
+            IterativeMinMaxLinearEquationSolver<double> impreciseSolver(std::make_unique<storm::solver::GeneralLinearEquationSolverFactory<double>>());
+            impreciseSolver.setMatrix(this->A->template toValueType<double>());
+            impreciseSolver.createLinearEquationSolver();
+            impreciseSolver.setCachingEnabled(true);
+
+            std::vector<double> impreciseX(x.size());
+            {
+                std::vector<ValueType> tmp(x.size());
+                this->createLowerBoundsVector(tmp);
+                auto targetIt = impreciseX.begin();
+                for (auto sourceIt = tmp.begin(); targetIt != impreciseX.end(); ++targetIt, ++sourceIt) {
+                    *targetIt = storm::utility::convertNumber<double, ValueType>(*sourceIt);
+                }
+            }
+            std::vector<double> impreciseTmpX(x.size());
+            std::vector<double> impreciseB(b.size());
+            auto targetIt = impreciseB.begin();
+            for (auto sourceIt = b.begin(); targetIt != impreciseB.end(); ++targetIt, ++sourceIt) {
+                *targetIt = storm::utility::convertNumber<double, ValueType>(*sourceIt);
+            }
+
+            std::vector<double>* currentX = &impreciseX;
+            std::vector<double>* newX = &impreciseTmpX;
+
+            Status status = Status::InProgress;
+            uint64_t overallIterations = 0;
+            ValueType precision = this->getSettings().getPrecision();
+            impreciseSolver.startMeasureProgress();
+            while (status == Status::InProgress && overallIterations < this->getSettings().getMaximalNumberOfIterations()) {
+                // Perform value iteration with the current precision.
+                IterativeMinMaxLinearEquationSolver<double>::ValueIterationResult result = impreciseSolver.performValueIteration(dir, currentX, newX, impreciseB, storm::utility::convertNumber<double, ValueType>(precision), this->getSettings().getRelativeTerminationCriterion(), SolverGuarantee::LessOrEqual, overallIterations);
+                
+                // Count the iterations.
+                overallIterations += result.iterations;
+                
+                // Try to sharpen the result.
+                uint64_t p = storm::utility::convertNumber<uint64_t>(storm::utility::ceil(storm::utility::log10<ValueType>(storm::utility::one<ValueType>() / precision)));
+                bool foundSolution = sharpen(dir, p, *this->A, *currentX, b, x);
+                
+                if (foundSolution) {
+                    status = Status::Converged;
+                } else {
+                    precision = precision / 10;
+                }
+            }
+            
+            reportStatus(status, overallIterations);
+            
+            if (!this->isCachingEnabled()) {
+                clearCache();
+            }
+            
+            return status == Status::Converged || status == Status::TerminatedEarly;
+        }
+        
+        template<>
+        bool IterativeMinMaxLinearEquationSolver<double>::solveEquationsRationalSearch(OptimizationDirection dir, std::vector<double>& x, std::vector<double> const& b) const {
+            // Create a rational representation of the input so we can check for a proper solution later.
+            storm::storage::SparseMatrix<storm::RationalNumber> rationalA = this->A->toValueType<storm::RationalNumber>();
+            std::vector<storm::RationalNumber> rationalX(x.size());
+            std::vector<storm::RationalNumber> rationalB = storm::utility::vector::convertNumericVector<storm::RationalNumber>(b);
+
+            if (!this->linEqSolverA) {
+                this->linEqSolverA = this->linearEquationSolverFactory->create(*this->A);
+                this->linEqSolverA->setCachingEnabled(true);
+            }
+            
+            if (!auxiliaryRowGroupVector) {
+                auxiliaryRowGroupVector = std::make_unique<std::vector<double>>(this->A->getRowGroupCount());
+            }
+            
+            std::vector<double>* newX = auxiliaryRowGroupVector.get();
+            std::vector<double>* currentX = &x;
+            
+            Status status = Status::InProgress;
+            uint64_t overallIterations = 0;
+            double precision = this->getSettings().getPrecision();
+            this->startMeasureProgress();
+            while (status == Status::InProgress && overallIterations < this->getSettings().getMaximalNumberOfIterations()) {
+                // Perform value iteration with the current precision.
+                ValueIterationResult result = performValueIteration(dir, currentX, newX, b, precision, this->getSettings().getRelativeTerminationCriterion(), SolverGuarantee::LessOrEqual, overallIterations);
+                
+                // Count the iterations.
+                overallIterations += result.iterations;
+                
+                // Try to sharpen the result.
+                uint64_t p = storm::utility::convertNumber<uint64_t>(storm::utility::ceil(storm::utility::log10<storm::RationalNumber>(storm::utility::one<storm::RationalNumber>() / storm::utility::convertNumber<storm::RationalNumber>(precision))));
+                bool foundSolution = sharpen(dir, p, rationalA, *currentX, rationalB, rationalX);
+                
+                if (foundSolution) {
+                    status = Status::Converged;
+                } else {
+                    precision = precision / 10;
+                }
+            }
+            
+            reportStatus(status, overallIterations);
+            
+            if (!this->isCachingEnabled()) {
+                clearCache();
+            }
+            
+            return status == Status::Converged || status == Status::TerminatedEarly;
+        }
+        
         template<typename ValueType>
         bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsAcyclic(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             uint64_t numGroups = this->A->getRowGroupCount();
diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h
index fb7fee8f0..4da82d8a2 100644
--- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h
+++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h
@@ -14,7 +14,7 @@ namespace storm {
             IterativeMinMaxLinearEquationSolverSettings();
             
             enum class SolutionMethod {
-                ValueIteration, PolicyIteration, Acyclic
+                ValueIteration, PolicyIteration, Acyclic, RationalSearch
             };
             
             void setSolutionMethod(SolutionMethod const& solutionMethod);
@@ -61,19 +61,44 @@ namespace storm {
             virtual MinMaxLinearEquationSolverRequirements getRequirements(EquationSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none) const override;
             
         private:
+            static bool isSolution(storm::OptimizationDirection dir, storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType> const& values, std::vector<ValueType> const& b);
+            
             bool solveEquationsPolicyIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+            bool valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const;
+
             bool solveEquationsValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
             bool solveEquationsSoundValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
             bool solveEquationsAcyclic(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+            bool solveEquationsRationalSearch(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+            static bool sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix<storm::RationalNumber> const& A, std::vector<double> const& x, std::vector<storm::RationalNumber> const& b, std::vector<storm::RationalNumber>& tmp);
+
+            template<typename ImpreciseValueType>
+            std::enable_if<std::is_same<ValueType, ImpreciseValueType>::value, bool> solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+            template<typename ImpreciseValueType>
+            std::enable_if<!std::is_same<ValueType, ImpreciseValueType>::value, bool> solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
 
-            bool valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const;
-            
             void computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, uint_fast64_t* choice = nullptr) const;
             
             enum class Status {
                 Converged, TerminatedEarly, MaximalIterationsExceeded, InProgress
             };
             
+            struct ValueIterationResult {
+                ValueIterationResult(uint64_t iterations, Status status) : iterations(iterations), status(status) {
+                    // Intentionally left empty.
+                }
+                
+                uint64_t iterations;
+                Status status;
+            };
+            
+            template <typename ValueTypePrime>
+            friend class IterativeMinMaxLinearEquationSolver;
+            
+            ValueIterationResult performValueIteration(OptimizationDirection dir, std::vector<ValueType>*& currentX, std::vector<ValueType>*& newX, std::vector<ValueType> const& b, ValueType const& precision, bool relative, SolverGuarantee const& guarantee, uint64_t currentIterations) const;
+            
+            void createLinearEquationSolver() const;
+            
             // possibly cached data
             mutable std::unique_ptr<std::vector<ValueType>> auxiliaryRowGroupVector; // A.rowGroupCount() entries
             mutable std::unique_ptr<std::vector<ValueType>> auxiliaryRowGroupVector2; // A.rowGroupCount() entries
diff --git a/src/storm/solver/MinMaxLinearEquationSolver.cpp b/src/storm/solver/MinMaxLinearEquationSolver.cpp
index ca6366869..9a7a3f357 100644
--- a/src/storm/solver/MinMaxLinearEquationSolver.cpp
+++ b/src/storm/solver/MinMaxLinearEquationSolver.cpp
@@ -228,7 +228,7 @@ namespace storm {
         std::unique_ptr<MinMaxLinearEquationSolver<ValueType>> GeneralMinMaxLinearEquationSolverFactory<ValueType>::create() const {
             std::unique_ptr<MinMaxLinearEquationSolver<ValueType>> result;
             auto method = this->getMinMaxMethod();
-            if (method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::Acyclic) {
+            if (method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::Acyclic || method == MinMaxMethod::RationalSearch) {
                 IterativeMinMaxLinearEquationSolverSettings<ValueType> iterativeSolverSettings;
                 iterativeSolverSettings.setSolutionMethod(method);
                 result = std::make_unique<IterativeMinMaxLinearEquationSolver<ValueType>>(std::make_unique<GeneralLinearEquationSolverFactory<ValueType>>(), iterativeSolverSettings);
@@ -248,7 +248,7 @@ namespace storm {
         std::unique_ptr<MinMaxLinearEquationSolver<storm::RationalNumber>> GeneralMinMaxLinearEquationSolverFactory<storm::RationalNumber>::create() const {
             std::unique_ptr<MinMaxLinearEquationSolver<storm::RationalNumber>> result;
             auto method = this->getMinMaxMethod();
-            if (method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::Acyclic) {
+            if (method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::Acyclic || method == MinMaxMethod::RationalSearch) {
                 IterativeMinMaxLinearEquationSolverSettings<storm::RationalNumber> iterativeSolverSettings;
                 iterativeSolverSettings.setSolutionMethod(method);
                 result = std::make_unique<IterativeMinMaxLinearEquationSolver<storm::RationalNumber>>(std::make_unique<GeneralLinearEquationSolverFactory<storm::RationalNumber>>(), iterativeSolverSettings);
diff --git a/src/storm/solver/SolverSelectionOptions.cpp b/src/storm/solver/SolverSelectionOptions.cpp
index 5a516c8dd..a8d8cb793 100644
--- a/src/storm/solver/SolverSelectionOptions.cpp
+++ b/src/storm/solver/SolverSelectionOptions.cpp
@@ -14,6 +14,8 @@ namespace storm {
                     return "topological";
                 case MinMaxMethod::Acyclic:
                     return "acyclic";
+                case MinMaxMethod::RationalSearch:
+                    return "ratsearch";
             }
             return "invalid";
         }
diff --git a/src/storm/solver/SolverSelectionOptions.h b/src/storm/solver/SolverSelectionOptions.h
index 141bc8a0e..8be9f5041 100644
--- a/src/storm/solver/SolverSelectionOptions.h
+++ b/src/storm/solver/SolverSelectionOptions.h
@@ -6,7 +6,7 @@
 
 namespace storm {
     namespace solver {
-        ExtendEnumsWithSelectionField(MinMaxMethod, PolicyIteration, ValueIteration, LinearProgramming, Topological, Acyclic)
+        ExtendEnumsWithSelectionField(MinMaxMethod, PolicyIteration, ValueIteration, LinearProgramming, Topological, Acyclic, RationalSearch)
         ExtendEnumsWithSelectionField(GameMethod, PolicyIteration, ValueIteration)
         ExtendEnumsWithSelectionField(LraMethod, LinearProgramming, ValueIteration)
 
diff --git a/src/storm/solver/StandardMinMaxLinearEquationSolver.h b/src/storm/solver/StandardMinMaxLinearEquationSolver.h
index a8517109e..845466db4 100644
--- a/src/storm/solver/StandardMinMaxLinearEquationSolver.h
+++ b/src/storm/solver/StandardMinMaxLinearEquationSolver.h
@@ -39,7 +39,6 @@ namespace storm {
             // A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix
             // the reference refers to localA.
             storm::storage::SparseMatrix<ValueType> const* A;
-            
         };
      
         template<typename ValueType>
diff --git a/src/storm/storage/SparseMatrix.h b/src/storm/storage/SparseMatrix.h
index f3d613658..92cb31b00 100644
--- a/src/storm/storage/SparseMatrix.h
+++ b/src/storm/storage/SparseMatrix.h
@@ -14,6 +14,7 @@
 
 #include "storm/utility/OsDetection.h"
 #include "storm/utility/macros.h"
+#include "storm/utility/constants.h"
 #include "storm/adapters/RationalFunctionAdapter.h"
 
 // Forward declaration for adapter classes.
@@ -1070,7 +1071,7 @@ namespace storm {
 				newColumnsAndValues.resize(columnsAndValues.size());
                 std::transform(columnsAndValues.begin(), columnsAndValues.end(), newColumnsAndValues.begin(),
                                [] (MatrixEntry<SparseMatrix::index_type, ValueType> const& a) {
-                                   return MatrixEntry<SparseMatrix::index_type, NewValueType>(a.getColumn(), static_cast<NewValueType>(a.getValue()));
+                                   return MatrixEntry<SparseMatrix::index_type, NewValueType>(a.getColumn(), storm::utility::convertNumber<NewValueType, ValueType>(a.getValue()));
                                });
 
 				return SparseMatrix<NewValueType>(columnCount, std::move(newRowIndications), std::move(newColumnsAndValues), std::move(newRowGroupIndices));
diff --git a/src/storm/utility/constants.cpp b/src/storm/utility/constants.cpp
index 58771a38d..05a9de733 100644
--- a/src/storm/utility/constants.cpp
+++ b/src/storm/utility/constants.cpp
@@ -1,6 +1,7 @@
 #include "storm/utility/constants.h"
 
 #include <type_traits>
+#include <cmath>
 
 #include "storm/storage/sparse/StateType.h"
 #include "storm/storage/SparseMatrix.h"
@@ -10,7 +11,9 @@
 #include "storm/exceptions/InvalidArgumentException.h"
 
 #include "storm/adapters/RationalFunctionAdapter.h"
+
 #include "storm/utility/macros.h"
+#include "storm/exceptions/NotSupportedException.h"
 
 namespace storm {
     namespace utility {
@@ -218,6 +221,31 @@ namespace storm {
         ValueType ceil(ValueType const& number) {
             return std::ceil(number);
         }
+        
+        template<typename ValueType>
+        ValueType log(ValueType const& number) {
+            return std::log(number);
+        }
+        
+        template<typename ValueType>
+        ValueType log10(ValueType const& number) {
+            return std::log10(number);
+        }
+        
+        template<typename ValueType>
+        ValueType mod(ValueType const& first, ValueType const& second) {
+            return std::fmod(first, second);
+        }
+        
+        template<typename ValueType>
+        ValueType numerator(ValueType const& number) {
+            STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Not supported.");
+        }
+        
+        template<typename ValueType>
+        ValueType denominator(ValueType const& number) {
+            STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Not supported.");
+        }
 
         template<typename ValueType>
         std::string to_string(ValueType const& value) {
@@ -330,10 +358,36 @@ namespace storm {
             return carl::ceil(number);
         }
         
+        template<>
+        ClnRationalNumber log(ClnRationalNumber const& number) {
+            return carl::log(number);
+        }
+        
+        template<>
+        ClnRationalNumber log10(ClnRationalNumber const& number) {
+            return carl::log10(number);
+        }
+
+        template<>
+        ClnRationalNumber mod(ClnRationalNumber const& first, ClnRationalNumber const& second) {
+            STORM_LOG_ASSERT(isInteger(first) && isInteger(second), "Expecting integers in modulo computation.");
+            return carl::mod(carl::getNum(first), carl::getNum(second));
+        }
+        
         template<>
         ClnRationalNumber pow(ClnRationalNumber const& value, uint_fast64_t exponent) {
             return carl::pow(value, exponent);
         }
+        
+        template<>
+        ClnRationalNumber numerator(ClnRationalNumber const& number) {
+            return carl::getNum(number);
+        }
+        
+        template<>
+        ClnRationalNumber denominator(ClnRationalNumber const& number) {
+            return carl::getDenom(number);
+        }
 #endif
         
 #if defined(STORM_HAVE_CARL) && defined(STORM_HAVE_GMP)
@@ -460,10 +514,36 @@ namespace storm {
             return carl::ceil(number);
         }
         
+        template<>
+        GmpRationalNumber log(GmpRationalNumber const& number) {
+            return carl::log(number);
+        }
+        
+        template<>
+        GmpRationalNumber log10(GmpRationalNumber const& number) {
+            return carl::log10(number);
+        }
+        
+        template<>
+        GmpRationalNumber mod(GmpRationalNumber const& first, GmpRationalNumber const& second) {
+            STORM_LOG_ASSERT(isInteger(first) && isInteger(second), "Expecting integers in modulo computation.");
+            return carl::mod(carl::getNum(first), carl::getNum(second));
+        }
+        
         template<>
         GmpRationalNumber pow(GmpRationalNumber const& value, uint_fast64_t exponent) {
             return carl::pow(value, exponent);
         }
+        
+        template<>
+        GmpRationalNumber numerator(GmpRationalNumber const& number) {
+            return carl::getNum(number);
+        }
+        
+        template<>
+        GmpRationalNumber denominator(GmpRationalNumber const& number) {
+            return carl::getDenom(number);
+        }
 #endif
         
 #if defined(STORM_HAVE_CARL) && defined(STORM_HAVE_GMP) && defined(STORM_HAVE_CLN)
@@ -689,6 +769,11 @@ namespace storm {
         template double abs(double const& number);
         template double floor(double const& number);
         template double ceil(double const& number);
+        template double log(double const& number);
+        template double log10(double const& number);
+        template double denominator(double const& number);
+        template double numerator(double const& number);
+        template double mod(double const& first, double const& second);
         template std::string to_string(double const& value);
 
         // float
@@ -717,6 +802,7 @@ namespace storm {
         template float abs(float const& number);
         template float floor(float const& number);
         template float ceil(float const& number);
+        template float log(float const& number);
         template std::string to_string(float const& value);
 
         // int
@@ -787,6 +873,7 @@ namespace storm {
         template storm::ClnRationalNumber abs(storm::ClnRationalNumber const& number);
         template storm::ClnRationalNumber floor(storm::ClnRationalNumber const& number);
         template storm::ClnRationalNumber ceil(storm::ClnRationalNumber const& number);
+        template storm::ClnRationalNumber log(storm::ClnRationalNumber const& number);
         template std::string to_string(storm::ClnRationalNumber const& value);
 #endif
         
@@ -822,6 +909,7 @@ namespace storm {
         template storm::GmpRationalNumber abs(storm::GmpRationalNumber const& number);
         template storm::GmpRationalNumber floor(storm::GmpRationalNumber const& number);
         template storm::GmpRationalNumber ceil(storm::GmpRationalNumber const& number);
+        template storm::GmpRationalNumber log(storm::GmpRationalNumber const& number);
         template std::string to_string(storm::GmpRationalNumber const& value);
 #endif
         
diff --git a/src/storm/utility/constants.h b/src/storm/utility/constants.h
index 5c0db58bc..8c47cc490 100644
--- a/src/storm/utility/constants.h
+++ b/src/storm/utility/constants.h
@@ -101,7 +101,22 @@ namespace storm {
 
         template<typename ValueType>
         ValueType ceil(ValueType const& number);
-        
+
+        template<typename ValueType>
+        ValueType log(ValueType const& number);
+
+        template<typename ValueType>
+        ValueType log10(ValueType const& number);
+
+        template<typename ValueType>
+        ValueType mod(ValueType const& first, ValueType const& second);
+
+        template<typename ValueType>
+        ValueType numerator(ValueType const& number);
+
+        template<typename ValueType>
+        ValueType denominator(ValueType const& number);
+
         template<typename ValueType>
         std::string to_string(ValueType const& value);
     }