diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6bdaa29a2..96b174c64 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -377,8 +377,10 @@ file(GLOB_RECURSE STORM_BUILDER_FILES ${PROJECT_SOURCE_DIR}/src/builder/*.h ${PR
 file(GLOB_RECURSE STORM_EXCEPTIONS_FILES ${PROJECT_SOURCE_DIR}/src/exceptions/*.h ${PROJECT_SOURCE_DIR}/src/exceptions/*.cpp)
 file(GLOB_RECURSE STORM_LOGIC_FILES ${PROJECT_SOURCE_DIR}/src/logic/*.h ${PROJECT_SOURCE_DIR}/src/logic/*.cpp)
 file(GLOB STORM_MODELCHECKER_FILES ${PROJECT_SOURCE_DIR}/src/modelchecker/*.h ${PROJECT_SOURCE_DIR}/src/modelchecker/*.cpp)
-file(GLOB_RECURSE STORM_MODELCHECKER_PRCTL_FILES ${PROJECT_SOURCE_DIR}/src/modelchecker/prctl/*.h ${PROJECT_SOURCE_DIR}/src/modelchecker/prctl/*.cpp)
-file(GLOB_RECURSE STORM_MODELCHECKER_CSL_FILES ${PROJECT_SOURCE_DIR}/src/modelchecker/csl/*.h ${PROJECT_SOURCE_DIR}/src/modelchecker/csl/*.cpp)
+file(GLOB STORM_MODELCHECKER_PRCTL_FILES ${PROJECT_SOURCE_DIR}/src/modelchecker/prctl/*.h ${PROJECT_SOURCE_DIR}/src/modelchecker/prctl/*.cpp)
+file(GLOB_RECURSE STORM_MODELCHECKER_PRCTL_HELPER_FILES ${PROJECT_SOURCE_DIR}/src/modelchecker/prctl/helper/*.h ${PROJECT_SOURCE_DIR}/src/modelchecker/prctl/helper/*.cpp)
+file(GLOB STORM_MODELCHECKER_CSL_FILES ${PROJECT_SOURCE_DIR}/src/modelchecker/csl/*.h ${PROJECT_SOURCE_DIR}/src/modelchecker/csl/*.cpp)
+file(GLOB_RECURSE STORM_MODELCHECKER_CSL_HELPER_FILES ${PROJECT_SOURCE_DIR}/src/modelchecker/csl/helper/*.h ${PROJECT_SOURCE_DIR}/src/modelchecker/csl/helper/*.cpp)
 file(GLOB_RECURSE STORM_MODELCHECKER_REACHABILITY_FILES ${PROJECT_SOURCE_DIR}/src/modelchecker/reachability/*.h ${PROJECT_SOURCE_DIR}/src/modelchecker/reachability/*.cpp)
 file(GLOB_RECURSE STORM_MODELCHECKER_PROPOSITIONAL_FILES ${PROJECT_SOURCE_DIR}/src/modelchecker/propositional/*.h ${PROJECT_SOURCE_DIR}/src/modelchecker/propositional/*.cpp)
 file(GLOB_RECURSE STORM_MODELCHECKER_RESULTS_FILES ${PROJECT_SOURCE_DIR}/src/modelchecker/results/*.h ${PROJECT_SOURCE_DIR}/src/modelchecker/results/*.cpp)
@@ -416,7 +418,9 @@ source_group(logic FILES ${STORM_LOGIC_FILES})
 source_group(generated FILES ${STORM_BUILD_HEADERS} ${STORM_BUILD_SOURCES})
 source_group(modelchecker FILES ${STORM_MODELCHECKER_FILES})
 source_group(modelchecker\\prctl FILES ${STORM_MODELCHECKER_PRCTL_FILES})
+source_group(modelchecker\\prctl\\helper FILES ${STORM_MODELCHECKER_PRCTL_HELPER_FILES})
 source_group(modelchecker\\csl FILES ${STORM_MODELCHECKER_CSL_FILES})
+source_group(modelchecker\\csl\\helper FILES ${STORM_MODELCHECKER_CSL_HELPER_FILES})
 source_group(modelchecker\\reachability FILES ${STORM_MODELCHECKER_REACHABILITY_FILES})
 source_group(modelchecker\\propositional FILES ${STORM_MODELCHECKER_PROPOSITIONAL_FILES})
 source_group(modelchecker\\results FILES ${STORM_MODELCHECKER_RESULTS_FILES})
diff --git a/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp b/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp
index 498a2fc02..428584385 100644
--- a/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp
+++ b/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp
@@ -1,5 +1,5 @@
 #include "src/modelchecker/csl/HybridCtmcCslModelChecker.h"
-#include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
+#include "src/modelchecker/csl/helper/SparseCtmcCslHelper.h"
 #include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h"
 
 #include "src/storage/dd/CuddOdd.h"
@@ -155,7 +155,7 @@ namespace storm {
                         
                         // Finally compute the transient probabilities.
                         std::vector<ValueType> values(statesWithProbabilityGreater0NonPsi.getNonZeroCount(), storm::utility::zero<ValueType>());
-                        std::vector<ValueType> subresult = SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, &explicitB, upperBound, uniformizationRate, values, *this->linearEquationSolverFactory);
+                        std::vector<ValueType> subresult = SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, &explicitB, upperBound, uniformizationRate, values, *this->linearEquationSolverFactory);
 
                         return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(),
                                                                                                       (psiStates || !statesWithProbabilityGreater0) && this->getModel().getReachableStates(),
@@ -195,7 +195,7 @@ namespace storm {
                         storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
                         
                         // Compute the transient probabilities.
-                        result = SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, result, *this->linearEquationSolverFactory);
+                        result = SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, result, *this->linearEquationSolverFactory);
                         
                         return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), !relevantStates && this->getModel().getReachableStates(), this->getModel().getManager().getAddZero(), relevantStates, odd, result));
                     } else {
@@ -221,7 +221,7 @@ namespace storm {
                             
                             // Compute the transient probabilities.
                             std::vector<ValueType> values(statesWithProbabilityGreater0NonPsi.getNonZeroCount(), storm::utility::zero<ValueType>());
-                            std::vector<ValueType> subResult = SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, &explicitB, upperBound - lowerBound, uniformizationRate, values, *this->linearEquationSolverFactory);
+                            std::vector<ValueType> subResult = SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, &explicitB, upperBound - lowerBound, uniformizationRate, values, *this->linearEquationSolverFactory);
 
                             // Transform the explicit result to a hybrid check result, so we can easily convert it to
                             // a symbolic qualitative format.
@@ -256,7 +256,7 @@ namespace storm {
                             uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates,  relevantStates, uniformizationRate);
                             explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
                             
-                            newSubresult = SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory);
+                            newSubresult = SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory);
                             
                             return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), !relevantStates && this->getModel().getReachableStates(), this->getModel().getManager().getAddZero(), relevantStates, odd, newSubresult));
                         } else {
@@ -276,7 +276,7 @@ namespace storm {
                             storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates,  statesWithProbabilityGreater0, uniformizationRate);
                             storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
                             
-                            newSubresult = SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory);
+                            newSubresult = SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory);
                             
                             return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), !statesWithProbabilityGreater0 && this->getModel().getReachableStates(), this->getModel().getManager().getAddZero(), statesWithProbabilityGreater0, odd, newSubresult));
                         }
@@ -311,7 +311,7 @@ namespace storm {
                 storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(),  this->getModel().getReachableStates(), uniformizationRate);
                 
                 storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
-                result = SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, result, *this->linearEquationSolverFactory);
+                result = SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, result, *this->linearEquationSolverFactory);
             }
             
             return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), odd, result));
@@ -353,7 +353,7 @@ namespace storm {
             std::vector<ValueType> explicitTotalRewardVector = totalRewardVector.template toVector<ValueType>(odd);
             
             // Finally, compute the transient probabilities.
-            std::vector<ValueType> result = SparseCtmcCslModelChecker<ValueType>::template computeTransientProbabilities<true>(explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, explicitTotalRewardVector, *this->linearEquationSolverFactory);
+            std::vector<ValueType> result = SparseCtmcCslHelper<ValueType>::template computeTransientProbabilities<true>(explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, explicitTotalRewardVector, *this->linearEquationSolverFactory);
             return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), std::move(odd), std::move(result)));
         }
         
@@ -370,7 +370,7 @@ namespace storm {
             storm::storage::SparseMatrix<ValueType> explicitProbabilityMatrix = probabilityMatrix.toMatrix(odd, odd);
             std::vector<ValueType> explicitExitRateVector = this->getModel().getExitRateVector().template toVector<ValueType>(odd);
             
-            std::vector<ValueType> result = SparseCtmcCslModelChecker<ValueType>::computeLongRunAverageHelper(explicitProbabilityMatrix, subResult.getTruthValuesVector().toVector(odd), &explicitExitRateVector, qualitative, *this->linearEquationSolverFactory);
+            std::vector<ValueType> result = SparseCtmcCslHelper<ValueType>::computeLongRunAverage(explicitProbabilityMatrix, subResult.getTruthValuesVector().toVector(odd), &explicitExitRateVector, qualitative, *this->linearEquationSolverFactory);
             return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), std::move(odd), std::move(result)));
         }
         
diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
index ccffd7b3d..1dba53698 100644
--- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
@@ -21,12 +21,12 @@
 namespace storm {
     namespace modelchecker {
         template <typename SparseCtmcModelType>
-        SparseCtmcCslModelChecker<SparseCtmcModelType>::SparseCtmcCslModelChecker(SparseCtmcModelType const& model) : SparsePropositionalModelChecker<ValueType>(model), linearEquationSolverFactory(new storm::utility::solver::LinearEquationSolverFactory<ValueType>()) {
+        SparseCtmcCslModelChecker<SparseCtmcModelType>::SparseCtmcCslModelChecker(SparseCtmcModelType const& model) : SparsePropositionalModelChecker<SparseCtmcModelType>(model), linearEquationSolverFactory(new storm::utility::solver::LinearEquationSolverFactory<ValueType>()) {
             // Intentionally left empty.
         }
         
         template <typename SparseCtmcModelType>
-        SparseCtmcCslModelChecker<SparseCtmcModelType>::SparseCtmcCslModelChecker(SparseCtmcModelType const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SparsePropositionalModelChecker<ValueType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
+        SparseCtmcCslModelChecker<SparseCtmcModelType>::SparseCtmcCslModelChecker(SparseCtmcModelType const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SparsePropositionalModelChecker<SparseCtmcModelType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
             // Intentionally left empty.
         }
         
@@ -58,7 +58,7 @@ namespace storm {
         std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
-            std::vector<ValueType> result = SparseDtmcPrctlModelChecker<ValueType>::computeNextProbabilitiesHelper(this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory);
+            std::vector<ValueType> result = SparseDtmcPrctlModelChecker<SparseCtmcModelType>::computeNextProbabilitiesHelper(this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory);
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(result)));
         }
         
@@ -72,12 +72,7 @@ namespace storm {
         }
         
         template <typename SparseCtmcModelType>
-        storm::models::sparse::Ctmc<ValueType> const& SparseCtmcCslModelChecker<SparseCtmcModelType>::getModel() const {
-            return this->template getModelAs<storm::models::sparse::Ctmc<ValueType>>();
-        }
-        
-        template <typename SparseCtmcModelType>
-        std::vector<ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& exitRates, bool qualitative, double lowerBound, double upperBound) const {
+        std::vector<typename SparseCtmcModelType::ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& exitRates, bool qualitative, double lowerBound, double upperBound) const {
             // If the time bounds are [0, inf], we rather call untimed reachability.
             storm::utility::ConstantsComparator<ValueType> comparator;
             if (comparator.isZero(lowerBound) && comparator.isInfinity(upperBound)) {
@@ -241,7 +236,7 @@ namespace storm {
         }
         
         template <typename SparseCtmcModelType>
-        storm::storage::SparseMatrix<ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeUniformizedMatrix(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& maybeStates, ValueType uniformizationRate, std::vector<ValueType> const& exitRates) {
+        storm::storage::SparseMatrix<typename SparseCtmcModelType::ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeUniformizedMatrix(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& maybeStates, ValueType uniformizationRate, std::vector<ValueType> const& exitRates) {
             STORM_LOG_DEBUG("Computing uniformized matrix using uniformization rate " << uniformizationRate << ".");
             STORM_LOG_DEBUG("Keeping " << maybeStates.getNumberOfSetBits() << " rows.");
             
@@ -269,7 +264,7 @@ namespace storm {
         
         template <typename SparseCtmcModelType>
         template<bool computeCumulativeReward>
-        std::vector<ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, std::vector<ValueType> const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+        std::vector<typename SparseCtmcModelType::ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, std::vector<ValueType> const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
             
             ValueType lambda = timeBound * uniformizationRate;
             
@@ -351,17 +346,17 @@ namespace storm {
         }
         
         template <typename SparseCtmcModelType>
-        std::vector<ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeUntilProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
-            return SparseDtmcPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(transitionMatrix, backwardTransitions, phiStates, psiStates, qualitative, linearEquationSolverFactory);
+        std::vector<typename SparseCtmcModelType::ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeUntilProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            return SparseDtmcPrctlModelChecker<SparseCtmcModelType>::computeUntilProbabilitiesHelper(transitionMatrix, backwardTransitions, phiStates, psiStates, qualitative, linearEquationSolverFactory);
         }
         
         template <typename SparseCtmcModelType>
-        std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+        std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeInstantaneousRewardsHelper(rewardPathFormula.getContinuousTimeBound())));
         }
         
         template <typename SparseCtmcModelType>
-        std::vector<ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeInstantaneousRewardsHelper(double timeBound) const {
+        std::vector<typename SparseCtmcModelType::ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeInstantaneousRewardsHelper(double timeBound) const {
             // Only compute the result if the model has a state-based reward this->getModel().
             STORM_LOG_THROW(this->getModel().hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
             
@@ -385,12 +380,12 @@ namespace storm {
         }
         
         template <typename SparseCtmcModelType>
-        std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+        std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeCumulativeRewardsHelper(rewardPathFormula.getContinuousTimeBound())));
         }
         
         template <typename SparseCtmcModelType>
-        std::vector<ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeCumulativeRewardsHelper(double timeBound) const {
+        std::vector<typename SparseCtmcModelType::ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeCumulativeRewardsHelper(double timeBound) const {
             // Only compute the result if the model has a state-based reward this->getModel().
             STORM_LOG_THROW(this->getModel().hasStateRewards() || this->getModel().hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
             
@@ -427,7 +422,7 @@ namespace storm {
         }
         
         template <typename SparseCtmcModelType>
-        storm::storage::SparseMatrix<ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeProbabilityMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates) {
+        storm::storage::SparseMatrix<typename SparseCtmcModelType::ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeProbabilityMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates) {
             // Turn the rates into probabilities by scaling each row with the exit rate of the state.
             storm::storage::SparseMatrix<ValueType> result(rateMatrix);
             for (uint_fast64_t row = 0; row < result.getRowCount(); ++row) {
@@ -439,7 +434,7 @@ namespace storm {
         }
         
         template <typename SparseCtmcModelType>
-        storm::storage::SparseMatrix<ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeGeneratorMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates) {
+        storm::storage::SparseMatrix<typename SparseCtmcModelType::ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeGeneratorMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates) {
             storm::storage::SparseMatrix<ValueType> generatorMatrix(rateMatrix, true);
             
             // Place the negative exit rate on the diagonal.
@@ -455,7 +450,7 @@ namespace storm {
         }
         
         template <typename SparseCtmcModelType>
-        std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+        std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
             storm::storage::SparseMatrix<ValueType> probabilityMatrix = computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector());
@@ -482,8 +477,8 @@ namespace storm {
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(computeLongRunAverageHelper(probabilityMatrix, subResult.getTruthValuesVector(), &this->getModel().getExitRateVector(), qualitative, *linearEquationSolverFactory)));
         }
         
-        template<typename ValueType>
-        std::vector<ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+        template <typename SparseCtmcModelType>
+        std::vector<typename SparseCtmcModelType::ValueType> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
             // If there are no goal states, we avoid the computation and directly return zero.
             uint_fast64_t numOfStates = transitionMatrix.getRowCount();
             if (psiStates.empty()) {
diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/modelchecker/csl/SparseCtmcCslModelChecker.h
index 498eea05b..d5fadf772 100644
--- a/src/modelchecker/csl/SparseCtmcCslModelChecker.h
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.h
@@ -20,9 +20,7 @@ namespace storm {
         class SparseCtmcCslModelChecker : public SparsePropositionalModelChecker<SparseCtmcModelType> {
         public:
             typedef typename SparseCtmcModelType::ValueType ValueType;
-            
-            friend class HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, ValueType>;
-            friend class SparseDtmcPrctlModelChecker<SparseCtmcModelType>;
+            typedef typename SparseCtmcModelType::RewardModelType RewardModelType;
             
             explicit SparseCtmcCslModelChecker(SparseCtmcModelType const& model);
             explicit SparseCtmcCslModelChecker(SparseCtmcModelType const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
@@ -32,17 +30,16 @@ namespace storm {
             virtual std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
 
         private:
-            // The methods that perform the actual checking.
             std::vector<ValueType> computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& exitRates, bool qualitative, double lowerBound, double upperBound) const;
             static std::vector<ValueType> computeUntilProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
-            std::vector<ValueType> computeInstantaneousRewardsHelper(double timeBound) const;
-            std::vector<ValueType> computeCumulativeRewardsHelper(double timeBound) const;
+            std::vector<ValueType> computeInstantaneousRewardsHelper(RewardModelType const& rewardModel, double timeBound) const;
+            std::vector<ValueType> computeCumulativeRewardsHelper(RewardModelType const& rewardModel, double timeBound) const;
             static std::vector<ValueType> computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
             /*!
              * Computes the matrix representing the transitions of the uniformized CTMC.
diff --git a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/modelchecker/csl/helper/SparseCtmcCslHelper.h b/src/modelchecker/csl/helper/SparseCtmcCslHelper.h
new file mode 100644
index 000000000..d2edcbd18
--- /dev/null
+++ b/src/modelchecker/csl/helper/SparseCtmcCslHelper.h
@@ -0,0 +1,25 @@
+#ifndef STORM_MODELCHECKER_SPARSE_CTMC_CSL_MODELCHECKER_HELPER_H_
+#define STORM_MODELCHECKER_SPARSE_CTMC_CSL_MODELCHECKER_HELPER_H_
+
+#include <vector>
+
+#include "src/models/sparse/StandardRewardModel.h"
+
+#include "src/storage/SparseMatrix.h"
+#include "src/storage/BitVector.h"
+
+#include "src/utility/solver.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            template <typename ValueType, typename RewardModelType = storm::models::sparse::StandardRewardModel<ValueType>>
+            class SparseCtmcCslHelper {
+            public:
+
+            }
+        }
+    }
+}
+
+#endif /* STORM_MODELCHECKER_SPARSE_CTMC_CSL_MODELCHECKER_HELPER_H_ */
\ No newline at end of file
diff --git a/src/modelchecker/csl/helper/SymbolicCtmcCslHelper.cpp b/src/modelchecker/csl/helper/SymbolicCtmcCslHelper.cpp
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/modelchecker/csl/helper/SymbolicCtmcCslHelper.h b/src/modelchecker/csl/helper/SymbolicCtmcCslHelper.h
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h
index 8964d29d7..016e2898a 100644
--- a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h
+++ b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h
@@ -23,9 +23,9 @@ namespace storm {
             virtual std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
 
         protected:
diff --git a/src/modelchecker/prctl/HybridMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/HybridMdpPrctlModelChecker.cpp
index f54bed203..a8af6665c 100644
--- a/src/modelchecker/prctl/HybridMdpPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/HybridMdpPrctlModelChecker.cpp
@@ -28,71 +28,7 @@ namespace storm {
         bool HybridMdpPrctlModelChecker<DdType, ValueType>::canHandle(storm::logic::Formula const& formula) const {
             return formula.isPctlStateFormula() || formula.isPctlPathFormula() || formula.isRewardPathFormula();
         }
-        
-        template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeUntilProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
-            // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
-            // probability 0 and 1 of satisfying the until-formula.
-            std::pair<storm::dd::Bdd<DdType>, storm::dd::Bdd<DdType>> statesWithProbability01;
-            if (minimize) {
-                statesWithProbability01 = storm::utility::graph::performProb01Min(model, phiStates, psiStates);
-            } else {
-                statesWithProbability01 = storm::utility::graph::performProb01Max(model, phiStates, psiStates);
-            }
-            storm::dd::Bdd<DdType> maybeStates = !statesWithProbability01.first && !statesWithProbability01.second && model.getReachableStates();
-            
-            // Perform some logging.
-            STORM_LOG_INFO("Found " << statesWithProbability01.first.getNonZeroCount() << " 'no' states.");
-            STORM_LOG_INFO("Found " << statesWithProbability01.second.getNonZeroCount() << " 'yes' states.");
-            STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
-            
-            // Check whether we need to compute exact probabilities for some states.
-            if (qualitative) {
-                // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1.
-                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd() + maybeStates.toAdd() * model.getManager().getConstant(0.5)));
-            } else {
-                // If there are maybe states, we need to solve an equation system.
-                if (!maybeStates.isZero()) {
-                    // Create the ODD for the translation between symbolic and explicit storage.
-                    storm::dd::Odd<DdType> odd(maybeStates);
-                    
-                    // Create the matrix and the vector for the equation system.
-                    storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
-                    
-                    // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
-                    // non-maybe states in the matrix.
-                    storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
-                    
-                    // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
-                    // maybe states.
-                    storm::dd::Add<DdType> prob1StatesAsColumn = statesWithProbability01.second.toAdd();
-                    prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs());
-                    storm::dd::Add<DdType> subvector = submatrix * prob1StatesAsColumn;
-                    subvector = subvector.sumAbstract(model.getColumnVariables());
-                    
-                    // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
-                    std::vector<uint_fast64_t> rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).toAdd().sumAbstract(model.getNondeterminismVariables()).template toVector<uint_fast64_t>(odd);
-                    
-                    // Finally cut away all columns targeting non-maybe states.
-                    submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
-                    
-                    // Create the solution vector.
-                    std::vector<ValueType> x(maybeStates.getNonZeroCount(), ValueType(0.5));
-                    
-                    // Translate the symbolic matrix/vector to their explicit representations and solve the equation system.
-                    std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
-                    
-                    std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitRepresentation.first);
-                    solver->solveEquationSystem(minimize, x, explicitRepresentation.second);
-                    
-                    // Return a hybrid check result that stores the numerical values explicitly.
-                    return std::unique_ptr<CheckResult>(new storm::modelchecker::HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getReachableStates() && !maybeStates, statesWithProbability01.second.toAdd(), maybeStates, odd, x));
-                } else {
-                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd()));
-                }
-            }
-        }
-        
+                
         template<storm::dd::DdType DdType, typename ValueType>
         std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
@@ -110,13 +46,7 @@ namespace storm {
             SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
             return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->computeNextProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector())));
         }
-        
-        template<storm::dd::DdType DdType, typename ValueType>
-        storm::dd::Add<DdType> HybridMdpPrctlModelChecker<DdType, ValueType>::computeNextProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates) {
-            storm::dd::Add<DdType> result = transitionMatrix * nextStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd();
-            return result.sumAbstract(model.getColumnVariables());
-        }
-        
+                
         template<storm::dd::DdType DdType, typename ValueType>
         std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
@@ -129,197 +59,27 @@ namespace storm {
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
-            // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
-            // probability 0 or 1 of satisfying the until-formula.
-            storm::dd::Bdd<DdType> statesWithProbabilityGreater0;
-            if (minimize) {
-                statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0A(model, transitionMatrix.notZero(), phiStates, psiStates);
-            } else {
-                statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0E(model, transitionMatrix.notZero(), phiStates, psiStates);
-            }
-            storm::dd::Bdd<DdType> maybeStates = statesWithProbabilityGreater0 && !psiStates && model.getReachableStates();
-                        
-            // If there are maybe states, we need to perform matrix-vector multiplications.
-            if (!maybeStates.isZero()) {
-                // Create the ODD for the translation between symbolic and explicit storage.
-                storm::dd::Odd<DdType> odd(maybeStates);
-                
-                // Create the matrix and the vector for the equation system.
-                storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
-                
-                // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
-                // non-maybe states in the matrix.
-                storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
-                
-                // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
-                // maybe states.
-                storm::dd::Add<DdType> prob1StatesAsColumn = psiStates.toAdd().swapVariables(model.getRowColumnMetaVariablePairs());
-                storm::dd::Add<DdType> subvector = (submatrix * prob1StatesAsColumn).sumAbstract(model.getColumnVariables());
-                
-                // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
-                std::vector<uint_fast64_t> rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).toAdd().sumAbstract(model.getNondeterminismVariables()).template toVector<uint_fast64_t>(odd);
-                
-                // Finally cut away all columns targeting non-maybe states.
-                submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
-                
-                // Create the solution vector.
-                std::vector<ValueType> x(maybeStates.getNonZeroCount(), storm::utility::zero<ValueType>());
-                
-                // Translate the symbolic matrix/vector to their explicit representations.
-                std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
-                
-                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitRepresentation.first);
-                solver->performMatrixVectorMultiplication(minimize, x, &explicitRepresentation.second, stepBound);
-                
-                // Return a hybrid check result that stores the numerical values explicitly.
-                return std::unique_ptr<CheckResult>(new storm::modelchecker::HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getReachableStates() && !maybeStates, psiStates.toAdd(), maybeStates, odd, x));
-            } else {
-                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd()));
-            }
-        }
-        
-        template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
             return this->computeCumulativeRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
         }
-        
-        template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeCumulativeRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
-            // Only compute the result if the model has at least one reward this->getModel().
-            STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
-            
-            // Compute the reward vector to add in each step based on the available reward models.
-            storm::dd::Add<DdType> totalRewardVector = model.hasStateRewards() ? model.getStateRewardVector() : model.getManager().getAddZero();
-            if (model.hasTransitionRewards()) {
-                totalRewardVector += (transitionMatrix * model.getTransitionRewardMatrix()).sumAbstract(model.getColumnVariables());
-            }
-            
-            // Create the ODD for the translation between symbolic and explicit storage.
-            storm::dd::Odd<DdType> odd(model.getReachableStates());
-            
-            // Create the solution vector.
-            std::vector<ValueType> x(model.getNumberOfStates(), storm::utility::zero<ValueType>());
-            
-            // Translate the symbolic matrix/vector to their explicit representations.
-            storm::storage::SparseMatrix<ValueType> explicitMatrix = transitionMatrix.toMatrix(model.getNondeterminismVariables(), odd, odd);
-            std::vector<ValueType> b = totalRewardVector.template toVector<ValueType>(model.getNondeterminismVariables(), odd, explicitMatrix.getRowGroupIndices());
-            
-            // Perform the matrix-vector multiplication.
-            std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitMatrix);
-            solver->performMatrixVectorMultiplication(minimize, x, &b, stepBound);
-            
-            // Return a hybrid check result that stores the numerical values explicitly.
-            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), odd, x));
-        }
-        
+                
         template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
             return this->computeInstantaneousRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
         }
-        
-        template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
-            // Only compute the result if the model has at least one reward this->getModel().
-            STORM_LOG_THROW(model.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
-            
-            // Create the ODD for the translation between symbolic and explicit storage.
-            storm::dd::Odd<DdType> odd(model.getReachableStates());
-
-            // Translate the symbolic matrix to its explicit representations.
-            storm::storage::SparseMatrix<ValueType> explicitMatrix = transitionMatrix.toMatrix(model.getNondeterminismVariables(), odd, odd);
-
-            // Create the solution vector (and initialize it to the state rewards of the model).
-            std::vector<ValueType> x = model.getStateRewardVector().template toVector<ValueType>(model.getNondeterminismVariables(), odd, explicitMatrix.getRowGroupIndices());
-            
-            // Perform the matrix-vector multiplication.
-            std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitMatrix);
-            solver->performMatrixVectorMultiplication(minimize, x, nullptr, stepBound);
-            
-            // Return a hybrid check result that stores the numerical values explicitly.
-            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), odd, x));
-        }
-        
+                
         template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
             SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
             return this->computeReachabilityRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative);
         }
         
-        template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeReachabilityRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative) {
-            
-            // Only compute the result if there is at least one reward model.
-            STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
-            
-            // Determine which states have a reward of infinity by definition.
-            storm::dd::Bdd<DdType> infinityStates;
-            storm::dd::Bdd<DdType> transitionMatrixBdd = transitionMatrix.notZero();
-            if (minimize) {
-                infinityStates = storm::utility::graph::performProb1A(model, transitionMatrixBdd, model.getReachableStates(), targetStates, storm::utility::graph::performProbGreater0A(model, transitionMatrixBdd, model.getReachableStates(), targetStates));
-            } else {
-                infinityStates = storm::utility::graph::performProb1E(model, transitionMatrixBdd, model.getReachableStates(), targetStates, storm::utility::graph::performProbGreater0E(model, transitionMatrixBdd, model.getReachableStates(), targetStates));
-            }
-            infinityStates = !infinityStates && model.getReachableStates();
-            storm::dd::Bdd<DdType> maybeStates = (!targetStates && !infinityStates) && model.getReachableStates();
-            STORM_LOG_INFO("Found " << infinityStates.getNonZeroCount() << " 'infinity' states.");
-            STORM_LOG_INFO("Found " << targetStates.getNonZeroCount() << " 'target' states.");
-            STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
-            
-            // Check whether we need to compute exact rewards for some states.
-            if (qualitative) {
-                // Set the values for all maybe-states to 1 to indicate that their reward values
-                // are neither 0 nor infinity.
-                return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + maybeStates.toAdd() * model.getManager().getConstant(storm::utility::one<ValueType>())));
-            } else {
-                // If there are maybe states, we need to solve an equation system.
-                if (!maybeStates.isZero()) {
-                    // Create the ODD for the translation between symbolic and explicit storage.
-                    storm::dd::Odd<DdType> odd(maybeStates);
-                    
-                    // Create the matrix and the vector for the equation system.
-                    storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
-                    
-                    // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
-                    // non-maybe states in the matrix.
-                    storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
-                    
-                    // Then compute the state reward vector to use in the computation.
-                    storm::dd::Add<DdType> subvector = stateRewardVector ? maybeStatesAdd * stateRewardVector.get() : model.getManager().getAddZero();
-                    if (transitionRewardMatrix) {
-                        subvector += (submatrix * transitionRewardMatrix.get()).sumAbstract(model.getColumnVariables());
-                    }
-                    
-                    // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
-                    std::vector<uint_fast64_t> rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).toAdd().sumAbstract(model.getNondeterminismVariables()).template toVector<uint_fast64_t>(odd);
-                    
-                    // Finally cut away all columns targeting non-maybe states.
-                    submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
-                    
-                    // Create the solution vector.
-                    std::vector<ValueType> x(maybeStates.getNonZeroCount(), ValueType(0.5));
-                    
-                    // Translate the symbolic matrix/vector to their explicit representations.
-                    std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
-                    
-                    // Now solve the resulting equation system.
-                    std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitRepresentation.first);
-                    solver->solveEquationSystem(minimize, x, explicitRepresentation.second);
-                                        
-                    // Return a hybrid check result that stores the numerical values explicitly.
-                    return std::unique_ptr<CheckResult>(new storm::modelchecker::HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getReachableStates() && !maybeStates, infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()), maybeStates, odd, x));
-                } else {
-                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>())));
-                }
-            }
-        }
-        
         template<storm::dd::DdType DdType, typename ValueType>
         storm::models::symbolic::Mdp<DdType> const& HybridMdpPrctlModelChecker<DdType, ValueType>::getModel() const {
             return this->template getModelAs<storm::models::symbolic::Mdp<DdType>>();
diff --git a/src/modelchecker/prctl/HybridMdpPrctlModelChecker.h b/src/modelchecker/prctl/HybridMdpPrctlModelChecker.h
index 24a33ee79..ee22ff766 100644
--- a/src/modelchecker/prctl/HybridMdpPrctlModelChecker.h
+++ b/src/modelchecker/prctl/HybridMdpPrctlModelChecker.h
@@ -18,22 +18,14 @@ namespace storm {
             virtual std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             
         protected:
             storm::models::symbolic::Mdp<DdType> const& getModel() const override;
             
         private:
-            // The methods that perform the actual checking.
-            static std::unique_ptr<CheckResult> computeBoundedUntilProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
-            static storm::dd::Add<DdType> computeNextProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates);
-            static std::unique_ptr<CheckResult> computeUntilProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
-            static std::unique_ptr<CheckResult> computeCumulativeRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
-            static std::unique_ptr<CheckResult> computeInstantaneousRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
-            static std::unique_ptr<CheckResult> computeReachabilityRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative);
-            
             // An object that is used for retrieving linear equation solvers.
             std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
         };
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
index d8d2f1378..ec9743ffe 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
@@ -4,17 +4,11 @@
 #include <memory>
 
 #include "src/utility/macros.h"
-#include "src/utility/vector.h"
-#include "src/utility/graph.h"
-#include "src/utility/solver.h"
 
-#include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
 #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 
-#include "src/exceptions/InvalidStateException.h"
-
-#include "src/exceptions/InvalidPropertyException.h"
+#include "src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h"
 
 namespace storm {
     namespace modelchecker {
@@ -33,58 +27,15 @@ namespace storm {
             return formula.isPctlStateFormula() || formula.isPctlPathFormula() || formula.isRewardPathFormula();
         }
         
-        template<typename SparseDtmcModelType>
-        std::vector<typename SparseDtmcPrctlModelChecker<SparseDtmcModelType>::ValueType> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound) const {
-            std::vector<ValueType> result(this->getModel().getNumberOfStates(), storm::utility::zero<ValueType>());
-            
-            // If we identify the states that have probability 0 of reaching the target states, we can exclude them in the further analysis.
-            storm::storage::BitVector maybeStates = storm::utility::graph::performProbGreater0(this->getModel().getBackwardTransitions(), phiStates, psiStates, true, stepBound);
-            maybeStates &= ~psiStates;
-            STORM_LOG_INFO("Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
-            
-            if (!maybeStates.empty()) {
-                // We can eliminate the rows and columns from the original transition probability matrix that have probability 0.
-                storm::storage::SparseMatrix<ValueType> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(true, maybeStates, maybeStates, true);
-                
-                // Create the vector of one-step probabilities to go to target states.
-                std::vector<ValueType> b = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(maybeStates, psiStates);
-                
-                // Create the vector with which to multiply.
-                std::vector<ValueType> subresult(maybeStates.getNumberOfSetBits());
-                
-                // Perform the matrix vector multiplication as often as required by the formula bound.
-                STORM_LOG_THROW(linearEquationSolverFactory != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
-                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(submatrix);
-                solver->performMatrixVectorMultiplication(subresult, &b, stepBound);
-                
-                // Set the values of the resulting vector accordingly.
-                storm::utility::vector::setVectorValues(result, maybeStates, subresult);
-            }
-            storm::utility::vector::setVectorValues<ValueType>(result, psiStates, storm::utility::one<ValueType>());
-            
-            return result;
-        }
-        
         template<typename SparseDtmcModelType>
         std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(pathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
             std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
-            ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();;
+            ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();
             ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();
-            std::unique_ptr<CheckResult> result = std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeBoundedUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound())));
-            return result;
-        }
-        
-        template<typename SparseDtmcModelType>
-        std::vector<typename SparseDtmcPrctlModelChecker<SparseDtmcModelType>::ValueType> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeNextProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
-            // Create the vector with which to multiply and initialize it correctly.
-            std::vector<ValueType> result(transitionMatrix.getRowCount());
-            storm::utility::vector::setVectorValues(result, nextStates, storm::utility::one<ValueType>());
-            
-            // Perform one single matrix-vector multiplication.
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix);
-            solver->performMatrixVectorMultiplication(result);
+            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper<ValueType>::computeBoundedUntilProbabilities(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *linearEquationSolverFactory);
+            std::unique_ptr<CheckResult> result = std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
             return result;
         }
         
@@ -92,189 +43,50 @@ namespace storm {
         std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeNextProbabilitiesHelper(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory)));
-        }
-        
-        template<typename SparseDtmcModelType>
-        std::vector<typename SparseDtmcPrctlModelChecker<SparseDtmcModelType>::ValueType> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeUntilProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
-            // We need to identify the states which have to be taken out of the matrix, i.e.
-            // all states that have probability 0 and 1 of satisfying the until-formula.
-            std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(backwardTransitions, phiStates, psiStates);
-            storm::storage::BitVector statesWithProbability0 = std::move(statesWithProbability01.first);
-            storm::storage::BitVector statesWithProbability1 = std::move(statesWithProbability01.second);
-            
-            // Perform some logging.
-            storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
-            STORM_LOG_INFO("Found " << statesWithProbability0.getNumberOfSetBits() << " 'no' states.");
-            STORM_LOG_INFO("Found " << statesWithProbability1.getNumberOfSetBits() << " 'yes' states.");
-            STORM_LOG_INFO("Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
-            
-            // Create resulting vector.
-            std::vector<ValueType> result(transitionMatrix.getRowCount());
-            
-            // Check whether we need to compute exact probabilities for some states.
-            if (qualitative) {
-                // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1.
-                storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, ValueType(0.5));
-            } else {
-                if (!maybeStates.empty()) {
-                    // In this case we have have to compute the probabilities.
-                    
-                    // We can eliminate the rows and columns from the original transition probability matrix.
-                    storm::storage::SparseMatrix<ValueType> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, true);
-                    
-                    // Converting the matrix from the fixpoint notation to the form needed for the equation
-                    // system. That is, we go from x = A*x + b to (I-A)x = b.
-                    submatrix.convertToEquationSystem();
-                    
-                    // Initialize the x vector with 0.5 for each element. This is the initial guess for
-                    // the iterative solvers. It should be safe as for all 'maybe' states we know that the
-                    // probability is strictly larger than 0.
-                    std::vector<ValueType> x(maybeStates.getNumberOfSetBits(), ValueType(0.5));
-                    
-                    // Prepare the right-hand side of the equation system. For entry i this corresponds to
-                    // the accumulated probability of going from state i to some 'yes' state.
-                    std::vector<ValueType> b = transitionMatrix.getConstrainedRowSumVector(maybeStates, statesWithProbability1);
-                    
-                    // Now solve the created system of linear equations.
-                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(submatrix);
-                    solver->solveEquationSystem(x, b);
-                    
-                    // Set values of resulting vector according to result.
-                    storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
-                }
-            }
-            
-            // Set values of resulting vector that are known exactly.
-            storm::utility::vector::setVectorValues<ValueType>(result, statesWithProbability0, storm::utility::zero<ValueType>());
-            storm::utility::vector::setVectorValues<ValueType>(result, statesWithProbability1, storm::utility::one<ValueType>());
-            
-            return result;
+            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper<ValueType>::computeNextProbabilities(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), *linearEquationSolverFactory);
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
         }
         
         template<typename SparseDtmcModelType>
         std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
-            ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();;
-            ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();;
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory)));
-        }
-        
-        template<typename SparseDtmcModelType>
-        std::vector<typename SparseDtmcPrctlModelChecker<SparseDtmcModelType>::ValueType> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeCumulativeRewardsHelper(RewardModelType const& rewardModel, uint_fast64_t stepBound) const {
-            // Compute the reward vector to add in each step based on the available reward models.
-            std::vector<ValueType> totalRewardVector = rewardModel.getTotalRewardVector(this->getModel().getTransitionMatrix());
-            
-            // Initialize result to either the state rewards of the model or the null vector.
-            std::vector<ValueType> result = rewardModel.getTotalStateActionRewardVector(this->getModel().getNumberOfStates(), this->getModel().getTransitionMatrix().getRowGroupIndices());
-            
-            // Perform the matrix vector multiplication as often as required by the formula bound.
-            STORM_LOG_THROW(linearEquationSolverFactory != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(this->getModel().getTransitionMatrix());
-            solver->performMatrixVectorMultiplication(result, &totalRewardVector, stepBound);
-            
-            return result;
+            ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();
+            ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();
+            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper<ValueType>::computeUntilProbabilities(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory);
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
         }
         
         template<typename SparseDtmcModelType>
         std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeCumulativeRewardsHelper(rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), rewardPathFormula.getDiscreteTimeBound())));
-        }
-        
-        template<typename SparseDtmcModelType>
-        std::vector<typename SparseDtmcPrctlModelChecker<SparseDtmcModelType>::ValueType> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeInstantaneousRewardsHelper(RewardModelType const& rewardModel, uint_fast64_t stepCount) const {
-            // Only compute the result if the model has a state-based reward this->getModel().
-            STORM_LOG_THROW(rewardModel.hasStateRewards() || rewardModel.hasStateActionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
-            
-            // Initialize result to state rewards of the model.
-            std::vector<ValueType> result = rewardModel.getTotalStateActionRewardVector(this->getModel().getNumberOfStates(), this->getModel().getTransitionMatrix().getRowGroupIndices());
-            
-            // Perform the matrix vector multiplication as often as required by the formula bound.
-            STORM_LOG_THROW(linearEquationSolverFactory != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(this->getModel().getTransitionMatrix());
-            solver->performMatrixVectorMultiplication(result, nullptr, stepCount);
-            
-            return result;
+            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper<ValueType>::computeCumulativeRewards(this->getModel().getTransitionMatrix(), rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), rewardPathFormula.getDiscreteTimeBound(), *linearEquationSolverFactory);
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
         }
         
         template<typename SparseDtmcModelType>
         std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeInstantaneousRewardsHelper(rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), rewardPathFormula.getDiscreteTimeBound())));
-        }
-        
-        template<typename SparseDtmcModelType>
-        std::vector<typename SparseDtmcPrctlModelChecker<SparseDtmcModelType>::ValueType> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeReachabilityRewardsHelper(RewardModelType const& rewardModel, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative) {
-            
-            // Determine which states have a reward of infinity by definition.
-            storm::storage::BitVector trueStates(transitionMatrix.getRowCount(), true);
-            storm::storage::BitVector infinityStates = storm::utility::graph::performProb1(backwardTransitions, trueStates, targetStates);
-            infinityStates.complement();
-            storm::storage::BitVector maybeStates = ~targetStates & ~infinityStates;
-            STORM_LOG_INFO("Found " << infinityStates.getNumberOfSetBits() << " 'infinity' states.");
-            STORM_LOG_INFO("Found " << targetStates.getNumberOfSetBits() << " 'target' states.");
-            STORM_LOG_INFO("Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
-            
-            // Create resulting vector.
-            std::vector<ValueType> result(transitionMatrix.getRowCount(), storm::utility::zero<ValueType>());
-            
-            // Check whether we need to compute exact rewards for some states.
-            if (qualitative) {
-                // Set the values for all maybe-states to 1 to indicate that their reward values
-                // are neither 0 nor infinity.
-                storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, storm::utility::one<ValueType>());
-            } else {
-                // In this case we have to compute the reward values for the remaining states.
-                // We can eliminate the rows and columns from the original transition probability matrix.
-                storm::storage::SparseMatrix<ValueType> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, true);
-                
-                // Converting the matrix from the fixpoint notation to the form needed for the equation
-                // system. That is, we go from x = A*x + b to (I-A)x = b.
-                submatrix.convertToEquationSystem();
-                                
-                // Initialize the x vector with 1 for each element. This is the initial guess for
-                // the iterative solvers.
-                std::vector<ValueType> x(submatrix.getColumnCount(), storm::utility::one<ValueType>());
-                
-                // Prepare the right-hand side of the equation system.
-                std::vector<ValueType> b = rewardModel.getTotalRewardVector(submatrix.getRowCount(), transitionMatrix, maybeStates);
-                
-                // Now solve the resulting equation system.
-                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(submatrix);
-                solver->solveEquationSystem(x, b);
-                
-                // Set values of resulting vector according to result.
-                storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
-            }
-            
-            // Set values of resulting vector that are known exactly.
-            storm::utility::vector::setVectorValues(result, infinityStates, storm::utility::infinity<ValueType>());
-            
-            return result;
+            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper<ValueType>::computeInstantaneousRewards(this->getModel().getTransitionMatrix(), rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), rewardPathFormula.getDiscreteTimeBound(), *linearEquationSolverFactory);
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
         }
         
         template<typename SparseDtmcModelType>
         std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative)));
+            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper<ValueType>::computeReachabilityRewards(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory);
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
         }
 
         template<typename SparseDtmcModelType>
         std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
-            
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(computeLongRunAverageHelper(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory)));
+            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseCtmcCslHelper<ValueType>::computeLongRunAverage(this->getModel().getTransitionMatrix(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory);
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
         }
-        
-		template<typename SparseDtmcModelType>
-        std::vector<typename SparseDtmcPrctlModelChecker<SparseDtmcModelType>::ValueType> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
-            return SparseCtmcCslModelChecker<ValueType>::computeLongRunAverageHelper(transitionMatrix, psiStates, nullptr, qualitative, linearEquationSolverFactory);
-		}
-        
+                
         template class SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>>;
     }
 }
\ No newline at end of file
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
index f1d9744fd..cc574d289 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
@@ -12,9 +12,6 @@ namespace storm {
         template<storm::dd::DdType DdType, typename ValueType>
         class HybridDtmcPrctlModelChecker;
         
-        // Forward-declare CTMC model checker so we can make it a friend.
-        template<typename SparseModelType> class SparseCtmcCslModelChecker;
-        
         template<class SparseDtmcModelType>
         class SparseDtmcPrctlModelChecker : public SparsePropositionalModelChecker<SparseDtmcModelType> {
         public:
@@ -22,7 +19,6 @@ namespace storm {
             typedef typename SparseDtmcModelType::RewardModelType RewardModelType;
             
             friend class HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, ValueType>;
-            friend class SparseCtmcCslModelChecker<ValueType>;
             
             explicit SparseDtmcPrctlModelChecker(SparseDtmcModelType const& model);
             explicit SparseDtmcPrctlModelChecker(SparseDtmcModelType const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
@@ -32,21 +28,12 @@ namespace storm {
             virtual std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>());
-            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>());
-            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>());
+            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             
         private:
-            // The methods that perform the actual checking.
-            std::vector<ValueType> computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound) const;
-            static std::vector<ValueType> computeNextProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
-            static std::vector<ValueType> computeUntilProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
-            std::vector<ValueType> computeInstantaneousRewardsHelper(RewardModelType const& rewardModel, uint_fast64_t stepCount) const;
-            std::vector<ValueType> computeCumulativeRewardsHelper(RewardModelType const& rewardModel, uint_fast64_t stepBound) const;
-            static std::vector<ValueType> computeReachabilityRewardsHelper(RewardModelType const& rewardModel, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative);
-            static std::vector<ValueType> computeLongRunAverageHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
-
             // An object that is used for retrieving linear equation solvers.
             std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
         };
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
index b079ff63c..adceed6f8 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
@@ -36,40 +36,6 @@ namespace storm {
             return formula.isPctlStateFormula() || formula.isPctlPathFormula() || formula.isRewardPathFormula();
         }
         
-        template<typename SparseMdpModelType>
-        std::vector<typename SparseMdpModelType::ValueType> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeBoundedUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound) const {
-            std::vector<ValueType> result(this->getModel().getNumberOfStates(), storm::utility::zero<ValueType>());
-            
-            // Determine the states that have 0 probability of reaching the target states.
-            storm::storage::BitVector maybeStates;
-            if (minimize) {
-                maybeStates = storm::utility::graph::performProbGreater0A(this->getModel().getTransitionMatrix(), this->getModel().getTransitionMatrix().getRowGroupIndices(), this->getModel().getBackwardTransitions(), phiStates, psiStates, true, stepBound);
-            } else {
-                maybeStates = storm::utility::graph::performProbGreater0E(this->getModel().getTransitionMatrix(), this->getModel().getTransitionMatrix().getRowGroupIndices(), this->getModel().getBackwardTransitions(), phiStates, psiStates, true, stepBound);
-            }
-            maybeStates &= ~psiStates;
-            STORM_LOG_INFO("Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
-            
-            if (!maybeStates.empty()) {
-                // We can eliminate the rows and columns from the original transition probability matrix that have probability 0.
-                storm::storage::SparseMatrix<ValueType> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(true, maybeStates, maybeStates, false);
-                std::vector<ValueType> b = this->getModel().getTransitionMatrix().getConstrainedRowGroupSumVector(maybeStates, psiStates);
-                
-                // Create the vector with which to multiply.
-                std::vector<ValueType> subresult(maybeStates.getNumberOfSetBits());
-            
-                STORM_LOG_THROW(MinMaxLinearEquationSolverFactory != nullptr, storm::exceptions::InvalidStateException, "No valid equation solver available.");
-                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory->create(submatrix);
-                solver->performMatrixVectorMultiplication(minimize, subresult, &b, stepBound);
-                
-                // Set the values of the resulting vector accordingly.
-                storm::utility::vector::setVectorValues(result, maybeStates, subresult);
-            }
-            storm::utility::vector::setVectorValues(result, psiStates, storm::utility::one<ValueType>());
-            
-            return result;
-        }
-        
         template<typename SparseMdpModelType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
@@ -81,19 +47,6 @@ namespace storm {
             return result;
         }
         
-        template<typename SparseMdpModelType>
-        std::vector<typename SparseMdpModelType::ValueType> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeNextProbabilitiesHelper(bool minimize, storm::storage::BitVector const& nextStates) {
-            // Create the vector with which to multiply and initialize it correctly.
-            std::vector<ValueType> result(this->getModel().getNumberOfStates());
-            storm::utility::vector::setVectorValues(result, nextStates, storm::utility::one<ValueType>());
-            
-            STORM_LOG_THROW(MinMaxLinearEquationSolverFactory != nullptr, storm::exceptions::InvalidStateException, "No valid equation solver available.");
-            std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory->create(this->getModel().getTransitionMatrix());
-            solver->performMatrixVectorMultiplication(minimize, result);
-            
-            return result;
-        }
-        
         template<typename SparseMdpModelType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
@@ -102,68 +55,6 @@ namespace storm {
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeNextProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector())));
         }
         
-        template<typename SparseMdpModelType>
-        std::vector<typename SparseMdpModelType::ValueType> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const {
-            return computeUntilProbabilitiesHelper(minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), phiStates, psiStates, *MinMaxLinearEquationSolverFactory, qualitative);
-        }
-        
-        template<typename SparseMdpModelType>
-        std::vector<typename SparseMdpModelType::ValueType> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeUntilProbabilitiesHelper(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& MinMaxLinearEquationSolverFactory, bool qualitative) {
-            size_t numberOfStates = phiStates.size();
-            
-            // We need to identify the states which have to be taken out of the matrix, i.e.
-            // all states that have probability 0 and 1 of satisfying the until-formula.
-            std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01;
-            if (minimize) {
-                statesWithProbability01 = storm::utility::graph::performProb01Min(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, phiStates, psiStates);
-            } else {
-                statesWithProbability01 = storm::utility::graph::performProb01Max(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, phiStates, psiStates);
-            }
-            storm::storage::BitVector statesWithProbability0 = std::move(statesWithProbability01.first);
-            storm::storage::BitVector statesWithProbability1 = std::move(statesWithProbability01.second);
-            storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
-            LOG4CPLUS_INFO(logger, "Found " << statesWithProbability0.getNumberOfSetBits() << " 'no' states.");
-            LOG4CPLUS_INFO(logger, "Found " << statesWithProbability1.getNumberOfSetBits() << " 'yes' states.");
-            LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
-            
-            // Create resulting vector.
-            std::vector<ValueType> result(numberOfStates);
-            
-            // Check whether we need to compute exact probabilities for some states.
-            if (qualitative) {
-                // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1.
-                storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, ValueType(0.5));
-            } else {
-                if (!maybeStates.empty()) {
-                    // In this case we have have to compute the probabilities.
-
-                    // First, we can eliminate the rows and columns from the original transition probability matrix for states
-                    // whose probabilities are already known.
-                    storm::storage::SparseMatrix<ValueType> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, false);
-                    
-                    // Prepare the right-hand side of the equation system. For entry i this corresponds to
-                    // the accumulated probability of going from state i to some 'yes' state.
-                    std::vector<ValueType> b = transitionMatrix.getConstrainedRowGroupSumVector(maybeStates, statesWithProbability1);
-                    
-                    // Create vector for results for maybe states.
-                    std::vector<ValueType> x(maybeStates.getNumberOfSetBits());
-                    
-                    // Solve the corresponding system of equations.
-                    std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory.create(submatrix);
-                    solver->solveEquationSystem(minimize, x, b);
-                    
-                    // Set values of resulting vector according to result.
-                    storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
-                }
-            }
-            
-            // Set values of resulting vector that are known exactly.
-            storm::utility::vector::setVectorValues<ValueType>(result, statesWithProbability0, storm::utility::zero<ValueType>());
-            storm::utility::vector::setVectorValues<ValueType>(result, statesWithProbability1, storm::utility::one<ValueType>());
-            
-            return result;
-        }
-        
         template<typename SparseMdpModelType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
@@ -174,28 +65,6 @@ namespace storm {
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(SparseMdpPrctlModelChecker<SparseMdpModelType>::computeUntilProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), *MinMaxLinearEquationSolverFactory, qualitative)));
         }
         
-        template<typename SparseMdpModelType>
-        std::vector<typename SparseMdpModelType::ValueType> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeCumulativeRewardsHelper(RewardModelType const& rewardModel, bool minimize, uint_fast64_t stepBound) const {
-            // Only compute the result if the model has at least one reward this->getModel().
-            STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
-            
-            // Compute the reward vector to add in each step based on the available reward models.
-            std::vector<ValueType> totalRewardVector = rewardModel.getTotalRewardVector(this->getModel().getTransitionMatrix());
-            
-            // Initialize result to either the state rewards of the model or the null vector.
-            std::vector<ValueType> result;
-            if (rewardModel.hasStateRewards()) {
-                result = rewardModel.getStateRewardVector();
-            } else {
-                result.resize(this->getModel().getNumberOfStates());
-            }
-            
-            std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory->create(this->getModel().getTransitionMatrix());
-            solver->performMatrixVectorMultiplication(minimize, result, &totalRewardVector, stepBound);
-            
-            return result;
-        }
-        
         template<typename SparseMdpModelType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
@@ -203,85 +72,13 @@ namespace storm {
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeCumulativeRewardsHelper(rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), optimalityType.get() == storm::logic::OptimalityType::Minimize, rewardPathFormula.getDiscreteTimeBound())));
         }
         
-        template<typename SparseMdpModelType>
-        std::vector<typename SparseMdpModelType::ValueType> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeInstantaneousRewardsHelper(RewardModelType const& rewardModel, bool minimize, uint_fast64_t stepCount) const {
-            // Only compute the result if the model has a state-based reward this->getModel().
-            STORM_LOG_THROW(rewardModel.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
-            
-            // Initialize result to state rewards of the this->getModel().
-            std::vector<ValueType> result(rewardModel.getStateRewardVector());
-            
-            STORM_LOG_THROW(MinMaxLinearEquationSolverFactory != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
-            std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory->create(this->getModel().getTransitionMatrix());
-            solver->performMatrixVectorMultiplication(minimize, result, nullptr, stepCount);
-            
-            return result;
-        }
-        
         template<typename SparseMdpModelType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeInstantaneousRewardsHelper(rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), optimalityType.get() == storm::logic::OptimalityType::Minimize, rewardPathFormula.getDiscreteTimeBound())));
         }
-        
-        template<typename SparseMdpModelType>
-        std::vector<typename SparseMdpModelType::ValueType> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeReachabilityRewardsHelper(RewardModelType const& rewardModel, bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& MinMaxLinearEquationSolverFactory, bool qualitative) const {
-            // Only compute the result if the model has at least one reward this->getModel().
-            STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
-            
-            // Determine which states have a reward of infinity by definition.
-            storm::storage::BitVector infinityStates;
-            storm::storage::BitVector trueStates(transitionMatrix.getRowCount(), true);
-            if (minimize) {
-                infinityStates = std::move(storm::utility::graph::performProb1A(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, trueStates, targetStates));
-            } else {
-                infinityStates = std::move(storm::utility::graph::performProb1E(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, trueStates, targetStates));
-            }
-            infinityStates.complement();
-            storm::storage::BitVector maybeStates = ~targetStates & ~infinityStates;
-            LOG4CPLUS_INFO(logger, "Found " << infinityStates.getNumberOfSetBits() << " 'infinity' states.");
-            LOG4CPLUS_INFO(logger, "Found " << targetStates.getNumberOfSetBits() << " 'target' states.");
-            LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
-            
-            // Create resulting vector.
-            std::vector<ValueType> result(transitionMatrix.getRowCount());
-            
-            // Check whether we need to compute exact rewards for some states.
-            if (qualitative) {
-                LOG4CPLUS_INFO(logger, "The rewards for the initial states were determined in a preprocessing step. No exact rewards were computed.");
-                // Set the values for all maybe-states to 1 to indicate that their reward values
-                // are neither 0 nor infinity.
-                storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, storm::utility::one<ValueType>());
-            } else {
-                // In this case we have to compute the reward values for the remaining states.
-                
-                // We can eliminate the rows and columns from the original transition probability matrix for states
-                // whose reward values are already known.
-                storm::storage::SparseMatrix<ValueType> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, false);
                 
-                // Prepare the right-hand side of the equation system.
-                std::vector<ValueType> b = rewardModel.getTotalRewardVector(submatrix.getRowCount(), transitionMatrix, maybeStates);
-                
-                // Create vector for results for maybe states.
-                std::vector<ValueType> x(maybeStates.getNumberOfSetBits());
-                
-                // Solve the corresponding system of equations.
-                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory.create(submatrix);
-                solver->solveEquationSystem(minimize, x, b);
-                
-                
-                // Set values of resulting vector according to result.
-                storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
-            }
-            
-            // Set values of resulting vector that are known exactly.
-            storm::utility::vector::setVectorValues(result, targetStates, storm::utility::zero<ValueType>());
-            storm::utility::vector::setVectorValues(result, infinityStates, storm::utility::infinity<ValueType>());
-            
-            return result;
-        }
-        
         template<typename SparseMdpModelType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
@@ -289,158 +86,6 @@ namespace storm {
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->MinMaxLinearEquationSolverFactory, qualitative)));
         }
-        
-
-		template<typename SparseMdpModelType>
-		std::vector<typename SparseMdpModelType::ValueType> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const {
-			// If there are no goal states, we avoid the computation and directly return zero.
-			auto numOfStates = this->getModel().getNumberOfStates();
-			if (psiStates.empty()) {
-				return std::vector<ValueType>(numOfStates, storm::utility::zero<ValueType>());
-			}
-
-			// Likewise, if all bits are set, we can avoid the computation and set.
-			if ((~psiStates).empty()) {
-				return std::vector<ValueType>(numOfStates, storm::utility::one<ValueType>());
-			}
-
-			// Start by decomposing the MDP into its MECs.
-			storm::storage::MaximalEndComponentDecomposition<double> mecDecomposition(this->getModel());
-
-			// Get some data members for convenience.
-			typename storm::storage::SparseMatrix<ValueType> const& transitionMatrix = this->getModel().getTransitionMatrix();
-			std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
-			ValueType zero = storm::utility::zero<ValueType>();
-
-			//first calculate LRA for the Maximal End Components.
-			storm::storage::BitVector statesInMecs(numOfStates);
-			std::vector<uint_fast64_t> stateToMecIndexMap(transitionMatrix.getColumnCount());
-			std::vector<ValueType> lraValuesForEndComponents(mecDecomposition.size(), zero);
-
-			for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) {
-				storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex];
-
-				lraValuesForEndComponents[currentMecIndex] = computeLraForMaximalEndComponent(minimize, transitionMatrix, psiStates, mec);
-
-				// Gather information for later use.
-				for (auto const& stateChoicesPair : mec) {
-					statesInMecs.set(stateChoicesPair.first);
-					stateToMecIndexMap[stateChoicesPair.first] = currentMecIndex;
-				}
-			}
-
-			// For fast transition rewriting, we build some auxiliary data structures.
-			storm::storage::BitVector statesNotContainedInAnyMec = ~statesInMecs;
-			uint_fast64_t firstAuxiliaryStateIndex = statesNotContainedInAnyMec.getNumberOfSetBits();
-			uint_fast64_t lastStateNotInMecs = 0;
-			uint_fast64_t numberOfStatesNotInMecs = 0;
-			std::vector<uint_fast64_t> statesNotInMecsBeforeIndex;
-			statesNotInMecsBeforeIndex.reserve(this->getModel().getNumberOfStates());
-			for (auto state : statesNotContainedInAnyMec) {
-				while (lastStateNotInMecs <= state) {
-					statesNotInMecsBeforeIndex.push_back(numberOfStatesNotInMecs);
-					++lastStateNotInMecs;
-				}
-				++numberOfStatesNotInMecs;
-			}
-
-			// Finally, we are ready to create the SSP matrix and right-hand side of the SSP.
-			std::vector<ValueType> b;
-			typename storm::storage::SparseMatrixBuilder<ValueType> sspMatrixBuilder(0, 0, 0, false, true, numberOfStatesNotInMecs + mecDecomposition.size());
-
-			// If the source state is not contained in any MEC, we copy its choices (and perform the necessary modifications).
-			uint_fast64_t currentChoice = 0;
-			for (auto state : statesNotContainedInAnyMec) {
-				sspMatrixBuilder.newRowGroup(currentChoice);
-
-				for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice, ++currentChoice) {
-					std::vector<ValueType> auxiliaryStateToProbabilityMap(mecDecomposition.size());
-					b.push_back(storm::utility::zero<ValueType>());
-
-					for (auto element : transitionMatrix.getRow(choice)) {
-						if (statesNotContainedInAnyMec.get(element.getColumn())) {
-							// If the target state is not contained in an MEC, we can copy over the entry.
-							sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue());
-						} else {
-							// If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector
-							// so that we are able to write the cumulative probability to the MEC into the matrix.
-							auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue();
-						}
-					}
-
-					// Now insert all (cumulative) probability values that target an MEC.
-					for (uint_fast64_t mecIndex = 0; mecIndex < auxiliaryStateToProbabilityMap.size(); ++mecIndex) {
-						if (auxiliaryStateToProbabilityMap[mecIndex] != 0) {
-							sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + mecIndex, auxiliaryStateToProbabilityMap[mecIndex]);
-						}
-					}
-				}
-			}
-
-			// Now we are ready to construct the choices for the auxiliary states.
-			for (uint_fast64_t mecIndex = 0; mecIndex < mecDecomposition.size(); ++mecIndex) {
-				storm::storage::MaximalEndComponent const& mec = mecDecomposition[mecIndex];
-				sspMatrixBuilder.newRowGroup(currentChoice);
-
-				for (auto const& stateChoicesPair : mec) {
-					uint_fast64_t state = stateChoicesPair.first;
-					boost::container::flat_set<uint_fast64_t> const& choicesInMec = stateChoicesPair.second;
-
-					for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) {
-						std::vector<ValueType> auxiliaryStateToProbabilityMap(mecDecomposition.size());
-
-						// If the choice is not contained in the MEC itself, we have to add a similar distribution to the auxiliary state.
-						if (choicesInMec.find(choice) == choicesInMec.end()) {
-							b.push_back(storm::utility::zero<ValueType>());
-
-							for (auto element : transitionMatrix.getRow(choice)) {
-								if (statesNotContainedInAnyMec.get(element.getColumn())) {
-									// If the target state is not contained in an MEC, we can copy over the entry.
-									sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue());
-								} else {
-									// If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector
-									// so that we are able to write the cumulative probability to the MEC into the matrix.
-									auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue();
-								}
-							}
-
-							// Now insert all (cumulative) probability values that target an MEC.
-							for (uint_fast64_t targetMecIndex = 0; targetMecIndex < auxiliaryStateToProbabilityMap.size(); ++targetMecIndex) {
-								if (auxiliaryStateToProbabilityMap[targetMecIndex] != 0) {
-									sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]);
-								}
-							}
-
-							++currentChoice;
-						}
-					}
-				}
-
-				// For each auxiliary state, there is the option to achieve the reward value of the LRA associated with the MEC.
-				++currentChoice;
-				b.push_back(lraValuesForEndComponents[mecIndex]);
-			}
-
-			// Finalize the matrix and solve the corresponding system of equations.
-			storm::storage::SparseMatrix<ValueType> sspMatrix = sspMatrixBuilder.build(currentChoice);
-
-			std::vector<ValueType> sspResult(numberOfStatesNotInMecs + mecDecomposition.size());
-			std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory->create(sspMatrix);
-			solver->solveEquationSystem(minimize, sspResult, b);
-
-			// Prepare result vector.
-			std::vector<ValueType> result(this->getModel().getNumberOfStates());
-
-			// Set the values for states not contained in MECs.
-			storm::utility::vector::setVectorValues(result, statesNotContainedInAnyMec, sspResult);
-
-			// Set the values for all states in MECs.
-			for (auto state : statesInMecs) {
-				result[state] = sspResult[firstAuxiliaryStateIndex + stateToMecIndexMap[state]];
-			}
-
-			return result;
-		}
 
 		template<typename SparseMdpModelType>
 		std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<SparseMdpModelType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
@@ -451,50 +96,6 @@ namespace storm {
 			
 			return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeLongRunAverageHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector(), qualitative)));
 		}
-
-		template<typename SparseMdpModelType>
-		typename SparseMdpModelType::ValueType SparseMdpPrctlModelChecker<SparseMdpModelType>::computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::MaximalEndComponent const& mec) {
-			std::shared_ptr<storm::solver::LpSolver> solver = storm::utility::solver::getLpSolver("LRA for MEC");
-			solver->setModelSense(minimize ? storm::solver::LpSolver::ModelSense::Maximize : storm::solver::LpSolver::ModelSense::Minimize);
-
-			//// First, we need to create the variables for the problem.
-			std::map<uint_fast64_t, storm::expressions::Variable> stateToVariableMap;
-			for (auto const& stateChoicesPair : mec) {
-				std::string variableName = "h" + std::to_string(stateChoicesPair.first);
-				stateToVariableMap[stateChoicesPair.first] = solver->addUnboundedContinuousVariable(variableName);
-			}
-			storm::expressions::Variable lambda = solver->addUnboundedContinuousVariable("L", 1);
-			solver->update();
-
-			// Now we encode the problem as constraints.
-			for (auto const& stateChoicesPair : mec) {
-				uint_fast64_t state = stateChoicesPair.first;
-
-				// Now, based on the type of the state, create a suitable constraint.
-				for (auto choice : stateChoicesPair.second) {
-					storm::expressions::Expression constraint = -lambda;
-					ValueType r = 0;
-
-					for (auto element : transitionMatrix.getRow(choice)) {
-						constraint = constraint + stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue());
-						if (psiStates.get(element.getColumn())) {
-							r += element.getValue();
-						}
-					}
-					constraint = solver->getConstant(r) + constraint;
-
-					if (minimize) {
-						constraint = stateToVariableMap.at(state) <= constraint;
-					} else {
-						constraint = stateToVariableMap.at(state) >= constraint;
-					}
-					solver->addConstraint("state" + std::to_string(state) + "," + std::to_string(choice), constraint);
-				}
-			}
-
-			solver->optimize();
-			return solver->getContinuousValue(lambda);
-		}
                 
         template class SparseMdpPrctlModelChecker<storm::models::sparse::Mdp<double>>;
     }
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
index 4223c7e62..2bd9c091a 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
@@ -5,8 +5,6 @@
 #include "src/models/sparse/Mdp.h"
 #include "src/utility/solver.h"
 #include "src/solver/MinMaxLinearEquationSolver.h"
-#include "src/storage/MaximalEndComponent.h"
-
 
 namespace storm {
     namespace counterexamples {
@@ -18,18 +16,12 @@ namespace storm {
     }
     
     namespace modelchecker {
-        
-        // Forward-declare other model checkers to make them friend classes.
-        template<typename SparseMarkovAutomatonModelType>
-        class SparseMarkovAutomatonCslModelChecker;
-        
         template<class SparseMdpModelType>
         class SparseMdpPrctlModelChecker : public SparsePropositionalModelChecker<SparseMdpModelType> {
         public:
             typedef typename SparseMdpModelType::ValueType ValueType;
             typedef typename SparseMdpModelType::RewardModelType RewardModelType;
             
-            friend class SparseMarkovAutomatonCslModelChecker<SparseMdpModelType>;
             friend class storm::counterexamples::SMTMinimalCommandSetGenerator<ValueType>;
             friend class storm::counterexamples::MILPMinimalLabelSetGenerator<ValueType>;
             
@@ -47,20 +39,8 @@ namespace storm {
 			virtual std::unique_ptr<CheckResult> computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             
         private:
-            // The methods that perform the actual checking.
-            std::vector<ValueType> computeBoundedUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound) const;
-            std::vector<ValueType> computeNextProbabilitiesHelper(bool minimize, storm::storage::BitVector const& nextStates);
-            std::vector<ValueType> computeUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const;
-            static std::vector<ValueType> computeUntilProbabilitiesHelper(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& MinMaxLinearEquationSolverFactory, bool qualitative);
-            std::vector<ValueType> computeInstantaneousRewardsHelper(RewardModelType const& rewardModel, bool minimize, uint_fast64_t stepCount) const;
-            std::vector<ValueType> computeCumulativeRewardsHelper(RewardModelType const& rewardModel, bool minimize, uint_fast64_t stepBound) const;
-            std::vector<ValueType> computeReachabilityRewardsHelper(RewardModelType const& rewardModel, bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& MinMaxLinearEquationSolverFactory, bool qualitative) const;
-			std::vector<ValueType> computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const;
-            
-			static ValueType computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::MaximalEndComponent const& mec);
-
             // An object that is used for retrieving solvers for systems of linear equations that are the result of nondeterministic choices.
-            std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>> MinMaxLinearEquationSolverFactory;
+            std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>> minMaxLinearEquationSolverFactory;
         };
     } // namespace modelchecker
 } // namespace storm
diff --git a/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.cpp
index 022eca0af..891061edc 100644
--- a/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.cpp
@@ -1,13 +1,13 @@
 #include "src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h"
 
-#include "src/storage/dd/CuddOdd.h"
+#include "src/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.h"
+
+#include "src/storage/dd/Add.h"
 
 #include "src/utility/macros.h"
-#include "src/utility/graph.h"
 
 #include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
 #include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
-#include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
 
 #include "src/exceptions/InvalidStateException.h"
 #include "src/exceptions/InvalidPropertyException.h"
@@ -29,80 +29,24 @@ namespace storm {
             return formula.isPctlStateFormula() || formula.isPctlPathFormula() || formula.isRewardPathFormula();
         }
         
-        template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
-            // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
-            // probability 0 and 1 of satisfying the until-formula.
-            std::pair<storm::dd::Bdd<DdType>, storm::dd::Bdd<DdType>> statesWithProbability01 = storm::utility::graph::performProb01(model, transitionMatrix, phiStates, psiStates);
-            storm::dd::Bdd<DdType> maybeStates = !statesWithProbability01.first && !statesWithProbability01.second && model.getReachableStates();
-            
-            // Perform some logging.
-            STORM_LOG_INFO("Found " << statesWithProbability01.first.getNonZeroCount() << " 'no' states.");
-            STORM_LOG_INFO("Found " << statesWithProbability01.second.getNonZeroCount() << " 'yes' states.");
-            STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
-            
-            // Check whether we need to compute exact probabilities for some states.
-            if (qualitative) {
-                // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1.
-                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd() + maybeStates.toAdd() * model.getManager().getConstant(0.5)));
-            } else {
-                // If there are maybe states, we need to solve an equation system.
-                if (!maybeStates.isZero()) {
-                    // Create the ODD for the translation between symbolic and explicit storage.
-                    storm::dd::Odd<DdType> odd(maybeStates);
-                    
-                    // Create the matrix and the vector for the equation system.
-                    storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
-                    
-                    // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
-                    // non-maybe states in the matrix.
-                    storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
-                    
-                    // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
-                    // maybe states.
-                    storm::dd::Add<DdType> prob1StatesAsColumn = statesWithProbability01.second.toAdd();
-                    prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs());
-                    storm::dd::Add<DdType> subvector = submatrix * prob1StatesAsColumn;
-                    subvector = subvector.sumAbstract(model.getColumnVariables());
-                    
-                    // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
-                    // for solving the equation system (i.e. compute (I-A)).
-                    submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
-                    submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
-                    
-                    // Solve the equation system.
-                    std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
-                    storm::dd::Add<DdType> result = solver->solveEquationSystem(model.getManager().getConstant(0.5) * maybeStatesAdd, subvector);
-                    
-                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd() + result));
-                } else {
-                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd()));
-                }
-            }
-        }
-        
         template<storm::dd::DdType DdType, typename ValueType>
         std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>();
             SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>();
-            return this->computeUntilProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory);
+            storm::dd::Add<DdType> numericResult = storm::modelchecker::helper::SymbolicDtmcPrctlHelper<DdType, ValueType>::computeUntilProbabilities(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory);
+            return std::unique_ptr<SymbolicQuantitativeCheckResult<DdType>>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), numericResult));
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
         std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula());
             SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
-            return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->computeNextProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector())));
-        }
-        
-        template<storm::dd::DdType DdType, typename ValueType>
-        storm::dd::Add<DdType> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeNextProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates) {
-            storm::dd::Add<DdType> result = transitionMatrix * nextStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd();
-            return result.sumAbstract(model.getColumnVariables());
+            storm::dd::Add<DdType> numericResult = storm::modelchecker::helper::SymbolicDtmcPrctlHelper<DdType, ValueType>::computeNextProbabilities(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector());
+            return std::unique_ptr<SymbolicQuantitativeCheckResult<DdType>>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), numericResult));
         }
-        
+                
         template<storm::dd::DdType DdType, typename ValueType>
         std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(pathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
@@ -110,147 +54,30 @@ namespace storm {
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>();
             SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>();
-            return this->computeBoundedUntilProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+            storm::dd::Add<DdType> numericResult = storm::modelchecker::helper::SymbolicDtmcPrctlHelper<DdType, ValueType>::computeBoundedUntilProbabilities(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+            return std::unique_ptr<SymbolicQuantitativeCheckResult<DdType>>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), numericResult));
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
-            // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
-            // probability 0 or 1 of satisfying the until-formula.
-            storm::dd::Bdd<DdType> statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(model, transitionMatrix.notZero(), phiStates, psiStates, stepBound);
-            storm::dd::Bdd<DdType> maybeStates = statesWithProbabilityGreater0 && !psiStates && model.getReachableStates();
-            
-            // If there are maybe states, we need to perform matrix-vector multiplications.
-            if (!maybeStates.isZero()) {
-                // Create the ODD for the translation between symbolic and explicit storage.
-                storm::dd::Odd<DdType> odd(maybeStates);
-                
-                // Create the matrix and the vector for the equation system.
-                storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
-                
-                // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
-                // non-maybe states in the matrix.
-                storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
-                
-                // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
-                // maybe states.
-                storm::dd::Add<DdType> prob1StatesAsColumn = psiStates.toAdd().swapVariables(model.getRowColumnMetaVariablePairs());
-                storm::dd::Add<DdType> subvector = (submatrix * prob1StatesAsColumn).sumAbstract(model.getColumnVariables());
-                
-                // Finally cut away all columns targeting non-maybe states.
-                submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
-                
-                // Perform the matrix-vector multiplication.
-                std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
-                storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getManager().getAddZero(), &subvector, stepBound);
-                
-                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd() + result));
-            } else {
-                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd()));
-            }
-        }
-        
-        template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
-            return this->computeCumulativeRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
-        }
-        
-        template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeCumulativeRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
-            // Only compute the result if the model has at least one reward this->getModel().
-            STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
-            
-            // Compute the reward vector to add in each step based on the available reward models.
-            storm::dd::Add<DdType> totalRewardVector = model.hasStateRewards() ? model.getStateRewardVector() : model.getManager().getAddZero();
-            if (model.hasTransitionRewards()) {
-                totalRewardVector += (transitionMatrix * model.getTransitionRewardMatrix()).sumAbstract(model.getColumnVariables());
-            }
-            
-            // Perform the matrix-vector multiplication.
-            std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix, model.getReachableStates(), model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
-            storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getManager().getAddZero(), &totalRewardVector, stepBound);
-            
-            return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), result));
+            storm::dd::Add<DdType> numericResult = storm::modelchecker::helper::SymbolicDtmcPrctlHelper<DdType, ValueType>::computeCumulativeRewards(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+            return std::unique_ptr<SymbolicQuantitativeCheckResult<DdType>>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), numericResult));
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
-            return this->computeInstantaneousRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+            storm::dd::Add<DdType> numericResult = storm::modelchecker::helper::SymbolicDtmcPrctlHelper<DdType, ValueType>::computeInstantaneousRewards(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+            return std::unique_ptr<SymbolicQuantitativeCheckResult<DdType>>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), numericResult));
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
-            // Only compute the result if the model has at least one reward this->getModel().
-            STORM_LOG_THROW(model.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
-            
-            // Perform the matrix-vector multiplication.
-            std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix, model.getReachableStates(), model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
-            storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getStateRewardVector(), nullptr, stepBound);
-            
-            return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), result));
-        }
-        
-        template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
             SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
-            return this->computeReachabilityRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative);
-        }
-        
-        template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory, bool qualitative) {
-            
-            // Only compute the result if there is at least one reward model.
-            STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
-            
-            // Determine which states have a reward of infinity by definition.
-            storm::dd::Bdd<DdType> infinityStates = storm::utility::graph::performProb1(model, transitionMatrix.notZero(), model.getReachableStates(), targetStates);
-            infinityStates = !infinityStates && model.getReachableStates();
-            storm::dd::Bdd<DdType> maybeStates = (!targetStates && !infinityStates) && model.getReachableStates();
-            STORM_LOG_INFO("Found " << infinityStates.getNonZeroCount() << " 'infinity' states.");
-            STORM_LOG_INFO("Found " << targetStates.getNonZeroCount() << " 'target' states.");
-            STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
-            
-            // Check whether we need to compute exact rewards for some states.
-            if (qualitative) {
-                // Set the values for all maybe-states to 1 to indicate that their reward values
-                // are neither 0 nor infinity.
-                return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + maybeStates.toAdd() * model.getManager().getConstant(storm::utility::one<ValueType>())));
-            } else {
-                // If there are maybe states, we need to solve an equation system.
-                if (!maybeStates.isZero()) {
-                    // Create the ODD for the translation between symbolic and explicit storage.
-                    storm::dd::Odd<DdType> odd(maybeStates);
-                    
-                    // Create the matrix and the vector for the equation system.
-                    storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
-                    
-                    // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
-                    // non-maybe states in the matrix.
-                    storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
-                    
-                    // Then compute the state reward vector to use in the computation.
-                    storm::dd::Add<DdType> subvector = stateRewardVector ? maybeStatesAdd * stateRewardVector.get() : model.getManager().getAddZero();
-                    if (transitionRewardMatrix) {
-                        subvector += (submatrix * transitionRewardMatrix.get()).sumAbstract(model.getColumnVariables());
-                    }
-                    
-                    // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
-                    // for solving the equation system (i.e. compute (I-A)).
-                    submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
-                    submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
-                    
-                    // Solve the equation system.
-                    std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
-                    storm::dd::Add<DdType> result = solver->solveEquationSystem(model.getManager().getConstant(0.5) * maybeStatesAdd, subvector);
-                    
-                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + result));
-                } else {
-                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>())));
-                }
-            }
+            storm::dd::Add<DdType> numericResult = storm::modelchecker::helper::SymbolicDtmcPrctlHelper<DdType, ValueType>::computeReachabilityRewards(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory);
+            return std::unique_ptr<SymbolicQuantitativeCheckResult<DdType>>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), numericResult));
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
diff --git a/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h
index 9dca0335b..761427f5e 100644
--- a/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h
@@ -7,14 +7,9 @@
 
 namespace storm {
     namespace modelchecker {
-        template<storm::dd::DdType DdType, typename ValueType>
-        class SymbolicCtmcCslModelChecker;
-        
         template<storm::dd::DdType DdType, typename ValueType>
         class SymbolicDtmcPrctlModelChecker : public SymbolicPropositionalModelChecker<DdType> {
         public:
-            friend class SymbolicCtmcCslModelChecker<DdType, ValueType>;
-            
             explicit SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model);
             explicit SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory);
             
@@ -23,22 +18,14 @@ namespace storm {
             virtual std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
-            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional<std::string> const& rewardModelName = boost::optional<std::string>(), bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             
         protected:
             storm::models::symbolic::Dtmc<DdType> const& getModel() const override;
             
         private:
-            // The methods that perform the actual checking.
-            static std::unique_ptr<CheckResult> computeBoundedUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
-            static storm::dd::Add<DdType> computeNextProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates);
-            static std::unique_ptr<CheckResult> computeUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
-            static std::unique_ptr<CheckResult> computeCumulativeRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
-            static std::unique_ptr<CheckResult> computeInstantaneousRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
-            static std::unique_ptr<CheckResult> computeReachabilityRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory, bool qualitative);
-            
             // An object that is used for retrieving linear equation solvers.
             std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>> linearEquationSolverFactory;
         };
diff --git a/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp b/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h b/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp
new file mode 100644
index 000000000..23cebc26a
--- /dev/null
+++ b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp
@@ -0,0 +1,248 @@
+#include "src/modelchecker/prctl/helper/HybridMdpPrctlHelper.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeUntilProbabilities(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
+                // probability 0 and 1 of satisfying the until-formula.
+                std::pair<storm::dd::Bdd<DdType>, storm::dd::Bdd<DdType>> statesWithProbability01;
+                if (minimize) {
+                    statesWithProbability01 = storm::utility::graph::performProb01Min(model, phiStates, psiStates);
+                } else {
+                    statesWithProbability01 = storm::utility::graph::performProb01Max(model, phiStates, psiStates);
+                }
+                storm::dd::Bdd<DdType> maybeStates = !statesWithProbability01.first && !statesWithProbability01.second && model.getReachableStates();
+                
+                // Perform some logging.
+                STORM_LOG_INFO("Found " << statesWithProbability01.first.getNonZeroCount() << " 'no' states.");
+                STORM_LOG_INFO("Found " << statesWithProbability01.second.getNonZeroCount() << " 'yes' states.");
+                STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
+                
+                // Check whether we need to compute exact probabilities for some states.
+                if (qualitative) {
+                    // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1.
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd() + maybeStates.toAdd() * model.getManager().getConstant(0.5)));
+                } else {
+                    // If there are maybe states, we need to solve an equation system.
+                    if (!maybeStates.isZero()) {
+                        // Create the ODD for the translation between symbolic and explicit storage.
+                        storm::dd::Odd<DdType> odd(maybeStates);
+                        
+                        // Create the matrix and the vector for the equation system.
+                        storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                        
+                        // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                        // non-maybe states in the matrix.
+                        storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                        
+                        // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
+                        // maybe states.
+                        storm::dd::Add<DdType> prob1StatesAsColumn = statesWithProbability01.second.toAdd();
+                        prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs());
+                        storm::dd::Add<DdType> subvector = submatrix * prob1StatesAsColumn;
+                        subvector = subvector.sumAbstract(model.getColumnVariables());
+                        
+                        // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
+                        std::vector<uint_fast64_t> rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).toAdd().sumAbstract(model.getNondeterminismVariables()).template toVector<uint_fast64_t>(odd);
+                        
+                        // Finally cut away all columns targeting non-maybe states.
+                        submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                        
+                        // Create the solution vector.
+                        std::vector<ValueType> x(maybeStates.getNonZeroCount(), ValueType(0.5));
+                        
+                        // Translate the symbolic matrix/vector to their explicit representations and solve the equation system.
+                        std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
+                        
+                        std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitRepresentation.first);
+                        solver->solveEquationSystem(minimize, x, explicitRepresentation.second);
+                        
+                        // Return a hybrid check result that stores the numerical values explicitly.
+                        return std::unique_ptr<CheckResult>(new storm::modelchecker::HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getReachableStates() && !maybeStates, statesWithProbability01.second.toAdd(), maybeStates, odd, x));
+                    } else {
+                        return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd()));
+                    }
+                }
+            }
+
+            template<storm::dd::DdType DdType, typename ValueType>
+            storm::dd::Add<DdType> HybridMdpPrctlModelChecker<DdType, ValueType>::computeNextProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates) {
+                storm::dd::Add<DdType> result = transitionMatrix * nextStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd();
+                return result.sumAbstract(model.getColumnVariables());
+            }
+
+            template<storm::dd::DdType DdType, typename ValueType>
+            std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
+                // probability 0 or 1 of satisfying the until-formula.
+                storm::dd::Bdd<DdType> statesWithProbabilityGreater0;
+                if (minimize) {
+                    statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0A(model, transitionMatrix.notZero(), phiStates, psiStates);
+                } else {
+                    statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0E(model, transitionMatrix.notZero(), phiStates, psiStates);
+                }
+                storm::dd::Bdd<DdType> maybeStates = statesWithProbabilityGreater0 && !psiStates && model.getReachableStates();
+                
+                // If there are maybe states, we need to perform matrix-vector multiplications.
+                if (!maybeStates.isZero()) {
+                    // Create the ODD for the translation between symbolic and explicit storage.
+                    storm::dd::Odd<DdType> odd(maybeStates);
+                    
+                    // Create the matrix and the vector for the equation system.
+                    storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                    
+                    // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                    // non-maybe states in the matrix.
+                    storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                    
+                    // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
+                    // maybe states.
+                    storm::dd::Add<DdType> prob1StatesAsColumn = psiStates.toAdd().swapVariables(model.getRowColumnMetaVariablePairs());
+                    storm::dd::Add<DdType> subvector = (submatrix * prob1StatesAsColumn).sumAbstract(model.getColumnVariables());
+                    
+                    // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
+                    std::vector<uint_fast64_t> rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).toAdd().sumAbstract(model.getNondeterminismVariables()).template toVector<uint_fast64_t>(odd);
+                    
+                    // Finally cut away all columns targeting non-maybe states.
+                    submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                    
+                    // Create the solution vector.
+                    std::vector<ValueType> x(maybeStates.getNonZeroCount(), storm::utility::zero<ValueType>());
+                    
+                    // Translate the symbolic matrix/vector to their explicit representations.
+                    std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
+                    
+                    std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitRepresentation.first);
+                    solver->performMatrixVectorMultiplication(minimize, x, &explicitRepresentation.second, stepBound);
+                    
+                    // Return a hybrid check result that stores the numerical values explicitly.
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getReachableStates() && !maybeStates, psiStates.toAdd(), maybeStates, odd, x));
+                } else {
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd()));
+                }
+            }
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                // Only compute the result if the model has at least one reward this->getModel().
+                STORM_LOG_THROW(model.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // Create the ODD for the translation between symbolic and explicit storage.
+                storm::dd::Odd<DdType> odd(model.getReachableStates());
+                
+                // Translate the symbolic matrix to its explicit representations.
+                storm::storage::SparseMatrix<ValueType> explicitMatrix = transitionMatrix.toMatrix(model.getNondeterminismVariables(), odd, odd);
+                
+                // Create the solution vector (and initialize it to the state rewards of the model).
+                std::vector<ValueType> x = model.getStateRewardVector().template toVector<ValueType>(model.getNondeterminismVariables(), odd, explicitMatrix.getRowGroupIndices());
+                
+                // Perform the matrix-vector multiplication.
+                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitMatrix);
+                solver->performMatrixVectorMultiplication(minimize, x, nullptr, stepBound);
+                
+                // Return a hybrid check result that stores the numerical values explicitly.
+                return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), odd, x));
+            }
+
+            template<storm::dd::DdType DdType, typename ValueType>
+            std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeCumulativeRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                // Only compute the result if the model has at least one reward this->getModel().
+                STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // Compute the reward vector to add in each step based on the available reward models.
+                storm::dd::Add<DdType> totalRewardVector = model.hasStateRewards() ? model.getStateRewardVector() : model.getManager().getAddZero();
+                if (model.hasTransitionRewards()) {
+                    totalRewardVector += (transitionMatrix * model.getTransitionRewardMatrix()).sumAbstract(model.getColumnVariables());
+                }
+                
+                // Create the ODD for the translation between symbolic and explicit storage.
+                storm::dd::Odd<DdType> odd(model.getReachableStates());
+                
+                // Create the solution vector.
+                std::vector<ValueType> x(model.getNumberOfStates(), storm::utility::zero<ValueType>());
+                
+                // Translate the symbolic matrix/vector to their explicit representations.
+                storm::storage::SparseMatrix<ValueType> explicitMatrix = transitionMatrix.toMatrix(model.getNondeterminismVariables(), odd, odd);
+                std::vector<ValueType> b = totalRewardVector.template toVector<ValueType>(model.getNondeterminismVariables(), odd, explicitMatrix.getRowGroupIndices());
+                
+                // Perform the matrix-vector multiplication.
+                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitMatrix);
+                solver->performMatrixVectorMultiplication(minimize, x, &b, stepBound);
+                
+                // Return a hybrid check result that stores the numerical values explicitly.
+                return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), odd, x));
+            }
+
+            template<storm::dd::DdType DdType, typename ValueType>
+            std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeReachabilityRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative) {
+                
+                // Only compute the result if there is at least one reward model.
+                STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // Determine which states have a reward of infinity by definition.
+                storm::dd::Bdd<DdType> infinityStates;
+                storm::dd::Bdd<DdType> transitionMatrixBdd = transitionMatrix.notZero();
+                if (minimize) {
+                    infinityStates = storm::utility::graph::performProb1A(model, transitionMatrixBdd, model.getReachableStates(), targetStates, storm::utility::graph::performProbGreater0A(model, transitionMatrixBdd, model.getReachableStates(), targetStates));
+                } else {
+                    infinityStates = storm::utility::graph::performProb1E(model, transitionMatrixBdd, model.getReachableStates(), targetStates, storm::utility::graph::performProbGreater0E(model, transitionMatrixBdd, model.getReachableStates(), targetStates));
+                }
+                infinityStates = !infinityStates && model.getReachableStates();
+                storm::dd::Bdd<DdType> maybeStates = (!targetStates && !infinityStates) && model.getReachableStates();
+                STORM_LOG_INFO("Found " << infinityStates.getNonZeroCount() << " 'infinity' states.");
+                STORM_LOG_INFO("Found " << targetStates.getNonZeroCount() << " 'target' states.");
+                STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
+                
+                // Check whether we need to compute exact rewards for some states.
+                if (qualitative) {
+                    // Set the values for all maybe-states to 1 to indicate that their reward values
+                    // are neither 0 nor infinity.
+                    return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + maybeStates.toAdd() * model.getManager().getConstant(storm::utility::one<ValueType>())));
+                } else {
+                    // If there are maybe states, we need to solve an equation system.
+                    if (!maybeStates.isZero()) {
+                        // Create the ODD for the translation between symbolic and explicit storage.
+                        storm::dd::Odd<DdType> odd(maybeStates);
+                        
+                        // Create the matrix and the vector for the equation system.
+                        storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                        
+                        // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                        // non-maybe states in the matrix.
+                        storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                        
+                        // Then compute the state reward vector to use in the computation.
+                        storm::dd::Add<DdType> subvector = stateRewardVector ? maybeStatesAdd * stateRewardVector.get() : model.getManager().getAddZero();
+                        if (transitionRewardMatrix) {
+                            subvector += (submatrix * transitionRewardMatrix.get()).sumAbstract(model.getColumnVariables());
+                        }
+                        
+                        // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
+                        std::vector<uint_fast64_t> rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).toAdd().sumAbstract(model.getNondeterminismVariables()).template toVector<uint_fast64_t>(odd);
+                        
+                        // Finally cut away all columns targeting non-maybe states.
+                        submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                        
+                        // Create the solution vector.
+                        std::vector<ValueType> x(maybeStates.getNonZeroCount(), ValueType(0.5));
+                        
+                        // Translate the symbolic matrix/vector to their explicit representations.
+                        std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
+                        
+                        // Now solve the resulting equation system.
+                        std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitRepresentation.first);
+                        solver->solveEquationSystem(minimize, x, explicitRepresentation.second);
+                        
+                        // Return a hybrid check result that stores the numerical values explicitly.
+                        return std::unique_ptr<CheckResult>(new storm::modelchecker::HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getReachableStates() && !maybeStates, infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()), maybeStates, odd, x));
+                    } else {
+                        return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>())));
+                    }
+                }
+            }
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.h b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.h
new file mode 100644
index 000000000..e55b5762b
--- /dev/null
+++ b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.h
@@ -0,0 +1,34 @@
+#ifndef STORM_MODELCHECKER_SPARSE_MDP_PRCTL_MODELCHECKER_HELPER_H_
+#define STORM_MODELCHECKER_SPARSE_MDP_PRCTL_MODELCHECKER_HELPER_H_
+
+#include "src/models/symbolic/Model.h"
+
+#include "src/storage/dd/Add.h"
+#include "src/storage/dd/Bdd.h"
+
+#include "src/utility/solver.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            class HybridMdpPrctlHelper {
+            public:
+                static std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static storm::dd::Add<DdType> computeNextProbabilities(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates);
+                
+                static std::unique_ptr<CheckResult> computeUntilProbabilities(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> computeCumulativeRewards(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> computeInstantaneousRewards(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> computeReachabilityRewards(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+            }
+        }
+    }
+}
+
+#endif /* STORM_MODELCHECKER_SPARSE_MDP_PRCTL_MODELCHECKER_HELPER_H_ */
diff --git a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
new file mode 100644
index 000000000..0cb95666f
--- /dev/null
+++ b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
@@ -0,0 +1,196 @@
+#include "src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h"
+
+#include "src/utility/macros.h"
+#include "src/utility/vector.h"
+#include "src/utility/graph.h"
+
+#include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/InvalidPropertyException.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeBoundedUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                std::vector<ValueType> result(transitionMatrix.getRowCount(), storm::utility::zero<ValueType>());
+                
+                // If we identify the states that have probability 0 of reaching the target states, we can exclude them in the further analysis.
+                storm::storage::BitVector maybeStates = storm::utility::graph::performProbGreater0(backwardTransitions, phiStates, psiStates, true, stepBound);
+                maybeStates &= ~psiStates;
+                STORM_LOG_INFO("Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
+                
+                if (!maybeStates.empty()) {
+                    // We can eliminate the rows and columns from the original transition probability matrix that have probability 0.
+                    storm::storage::SparseMatrix<ValueType> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, true);
+                    
+                    // Create the vector of one-step probabilities to go to target states.
+                    std::vector<ValueType> b = transitionMatrix.getConstrainedRowSumVector(maybeStates, psiStates);
+                    
+                    // Create the vector with which to multiply.
+                    std::vector<ValueType> subresult(maybeStates.getNumberOfSetBits());
+                    
+                    // Perform the matrix vector multiplication as often as required by the formula bound.
+                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(submatrix);
+                    solver->performMatrixVectorMultiplication(subresult, &b, stepBound);
+                    
+                    // Set the values of the resulting vector accordingly.
+                    storm::utility::vector::setVectorValues(result, maybeStates, subresult);
+                }
+                storm::utility::vector::setVectorValues<ValueType>(result, psiStates, storm::utility::one<ValueType>());
+                
+                return result;
+            }
+        
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                // We need to identify the states which have to be taken out of the matrix, i.e.
+                // all states that have probability 0 and 1 of satisfying the until-formula.
+                std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(backwardTransitions, phiStates, psiStates);
+                storm::storage::BitVector statesWithProbability0 = std::move(statesWithProbability01.first);
+                storm::storage::BitVector statesWithProbability1 = std::move(statesWithProbability01.second);
+                
+                // Perform some logging.
+                storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
+                STORM_LOG_INFO("Found " << statesWithProbability0.getNumberOfSetBits() << " 'no' states.");
+                STORM_LOG_INFO("Found " << statesWithProbability1.getNumberOfSetBits() << " 'yes' states.");
+                STORM_LOG_INFO("Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
+                
+                // Create resulting vector.
+                std::vector<ValueType> result(transitionMatrix.getRowCount());
+                
+                // Check whether we need to compute exact probabilities for some states.
+                if (qualitative) {
+                    // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1.
+                    storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, ValueType(0.5));
+                } else {
+                    if (!maybeStates.empty()) {
+                        // In this case we have have to compute the probabilities.
+                        
+                        // We can eliminate the rows and columns from the original transition probability matrix.
+                        storm::storage::SparseMatrix<ValueType> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, true);
+                        
+                        // Converting the matrix from the fixpoint notation to the form needed for the equation
+                        // system. That is, we go from x = A*x + b to (I-A)x = b.
+                        submatrix.convertToEquationSystem();
+                        
+                        // Initialize the x vector with 0.5 for each element. This is the initial guess for
+                        // the iterative solvers. It should be safe as for all 'maybe' states we know that the
+                        // probability is strictly larger than 0.
+                        std::vector<ValueType> x(maybeStates.getNumberOfSetBits(), ValueType(0.5));
+                        
+                        // Prepare the right-hand side of the equation system. For entry i this corresponds to
+                        // the accumulated probability of going from state i to some 'yes' state.
+                        std::vector<ValueType> b = transitionMatrix.getConstrainedRowSumVector(maybeStates, statesWithProbability1);
+                        
+                        // Now solve the created system of linear equations.
+                        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(submatrix);
+                        solver->solveEquationSystem(x, b);
+                        
+                        // Set values of resulting vector according to result.
+                        storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
+                    }
+                }
+                
+                // Set values of resulting vector that are known exactly.
+                storm::utility::vector::setVectorValues<ValueType>(result, statesWithProbability0, storm::utility::zero<ValueType>());
+                storm::utility::vector::setVectorValues<ValueType>(result, statesWithProbability1, storm::utility::one<ValueType>());
+                
+                return result;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeNextProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                // Create the vector with which to multiply and initialize it correctly.
+                std::vector<ValueType> result(transitionMatrix.getRowCount());
+                storm::utility::vector::setVectorValues(result, nextStates, storm::utility::one<ValueType>());
+                
+                // Perform one single matrix-vector multiplication.
+                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix);
+                solver->performMatrixVectorMultiplication(result);
+                return result;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeCumulativeRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                // Compute the reward vector to add in each step based on the available reward models.
+                std::vector<ValueType> totalRewardVector = rewardModel.getTotalRewardVector(transitionMatrix);
+                
+                // Initialize result to either the state rewards of the model or the null vector.
+                std::vector<ValueType> result = rewardModel.getTotalStateActionRewardVector(transitionMatrix.getRowCount(), transitionMatrix.getRowGroupIndices());
+                
+                // Perform the matrix vector multiplication as often as required by the formula bound.
+                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix);
+                solver->performMatrixVectorMultiplication(result, &totalRewardVector, stepBound);
+                
+                return result;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeInstantaneousRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepCount, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                // Only compute the result if the model has a state-based reward this->getModel().
+                STORM_LOG_THROW(rewardModel.hasStateRewards() || rewardModel.hasStateActionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // Initialize result to state rewards of the model.
+                std::vector<ValueType> result = rewardModel.getTotalStateActionRewardVector(transitionMatrix.getRowCount(), transitionMatrix.getRowGroupIndices());
+                
+                // Perform the matrix vector multiplication as often as required by the formula bound.
+                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix);
+                solver->performMatrixVectorMultiplication(result, nullptr, stepCount);
+                
+                return result;
+            }
+
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+                
+                // Determine which states have a reward of infinity by definition.
+                storm::storage::BitVector trueStates(transitionMatrix.getRowCount(), true);
+                storm::storage::BitVector infinityStates = storm::utility::graph::performProb1(backwardTransitions, trueStates, targetStates);
+                infinityStates.complement();
+                storm::storage::BitVector maybeStates = ~targetStates & ~infinityStates;
+                STORM_LOG_INFO("Found " << infinityStates.getNumberOfSetBits() << " 'infinity' states.");
+                STORM_LOG_INFO("Found " << targetStates.getNumberOfSetBits() << " 'target' states.");
+                STORM_LOG_INFO("Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
+                
+                // Create resulting vector.
+                std::vector<ValueType> result(transitionMatrix.getRowCount(), storm::utility::zero<ValueType>());
+                
+                // Check whether we need to compute exact rewards for some states.
+                if (qualitative) {
+                    // Set the values for all maybe-states to 1 to indicate that their reward values
+                    // are neither 0 nor infinity.
+                    storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, storm::utility::one<ValueType>());
+                } else {
+                    // In this case we have to compute the reward values for the remaining states.
+                    // We can eliminate the rows and columns from the original transition probability matrix.
+                    storm::storage::SparseMatrix<ValueType> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, true);
+                    
+                    // Converting the matrix from the fixpoint notation to the form needed for the equation
+                    // system. That is, we go from x = A*x + b to (I-A)x = b.
+                    submatrix.convertToEquationSystem();
+                    
+                    // Initialize the x vector with 1 for each element. This is the initial guess for
+                    // the iterative solvers.
+                    std::vector<ValueType> x(submatrix.getColumnCount(), storm::utility::one<ValueType>());
+                    
+                    // Prepare the right-hand side of the equation system.
+                    std::vector<ValueType> b = rewardModel.getTotalRewardVector(submatrix.getRowCount(), transitionMatrix, maybeStates);
+                    
+                    // Now solve the resulting equation system.
+                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(submatrix);
+                    solver->solveEquationSystem(x, b);
+                    
+                    // Set values of resulting vector according to result.
+                    storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
+                }
+                
+                // Set values of resulting vector that are known exactly.
+                storm::utility::vector::setVectorValues(result, infinityStates, storm::utility::infinity<ValueType>());
+                
+                return result;
+            }
+         
+            template class SparseDtmcPrctlHelper<double>;
+        }
+    }
+}
diff --git a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h
new file mode 100644
index 000000000..57d2a9de7
--- /dev/null
+++ b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h
@@ -0,0 +1,39 @@
+#ifndef STORM_MODELCHECKER_SPARSE_DTMC_PRCTL_MODELCHECKER_HELPER_H_
+#define STORM_MODELCHECKER_SPARSE_DTMC_PRCTL_MODELCHECKER_HELPER_H_
+
+#include <vector>
+
+#include "src/models/sparse/StandardRewardModel.h"
+
+#include "src/storage/SparseMatrix.h"
+#include "src/storage/BitVector.h"
+
+#include "src/utility/solver.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            
+            template <typename ValueType, typename RewardModelType = storm::models::sparse::StandardRewardModel<ValueType>>
+            class SparseDtmcPrctlHelper {
+            public:
+                static std::vector<ValueType> computeBoundedUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeNextProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeCumulativeRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeInstantaneousRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepCount, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeLongRunAverage(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+            };
+            
+        }
+    }
+}
+
+#endif /* STORM_MODELCHECKER_SPARSE_DTMC_PRCTL_MODELCHECKER_HELPER_H_ */
\ No newline at end of file
diff --git a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
new file mode 100644
index 000000000..ff278e20e
--- /dev/null
+++ b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
@@ -0,0 +1,407 @@
+#include "src/modelchecker/prctl/helper/SparseMdpPrctlHelper.h"
+
+#include "src/storage/MaximalEndComponentDecomposition.h"
+
+#include "src/utility/macros.h"
+#include "src/utility/vector.h"
+#include "src/utility/graph.h"
+
+#include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/InvalidPropertyException.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMdpPrctlHelper<ValueType, RewardModelType>::computeBoundedUntilProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                std::vector<ValueType> result(transitionMatrix.getRowCount(), storm::utility::zero<ValueType>());
+                
+                // Determine the states that have 0 probability of reaching the target states.
+                storm::storage::BitVector maybeStates;
+                if (minimize) {
+                    maybeStates = storm::utility::graph::performProbGreater0A(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, phiStates, psiStates, true, stepBound);
+                } else {
+                    maybeStates = storm::utility::graph::performProbGreater0E(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, phiStates, psiStates, true, stepBound);
+                }
+                maybeStates &= ~psiStates;
+                STORM_LOG_INFO("Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
+                
+                if (!maybeStates.empty()) {
+                    // We can eliminate the rows and columns from the original transition probability matrix that have probability 0.
+                    storm::storage::SparseMatrix<ValueType> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, false);
+                    std::vector<ValueType> b = transitionMatrix.getConstrainedRowGroupSumVector(maybeStates, psiStates);
+                    
+                    // Create the vector with which to multiply.
+                    std::vector<ValueType> subresult(maybeStates.getNumberOfSetBits());
+                    
+                    std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(submatrix);
+                    solver->performMatrixVectorMultiplication(minimize, subresult, &b, stepBound);
+                    
+                    // Set the values of the resulting vector accordingly.
+                    storm::utility::vector::setVectorValues(result, maybeStates, subresult);
+                }
+                storm::utility::vector::setVectorValues(result, psiStates, storm::utility::one<ValueType>());
+                
+                return result;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMdpPrctlHelper<ValueType, RewardModelType>::computeNextProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                // Create the vector with which to multiply and initialize it correctly.
+                std::vector<ValueType> result(transitionMatrix.getRowCount());
+                storm::utility::vector::setVectorValues(result, nextStates, storm::utility::one<ValueType>());
+                
+                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix);
+                solver->performMatrixVectorMultiplication(minimize, result);
+                
+                return result;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMdpPrctlHelper<ValueType, RewardModelType>::computeUntilProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                uint_fast64_t numberOfStates = transitionMatrix.getRowCount();
+                
+                // We need to identify the states which have to be taken out of the matrix, i.e.
+                // all states that have probability 0 and 1 of satisfying the until-formula.
+                std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01;
+                if (minimize) {
+                    statesWithProbability01 = storm::utility::graph::performProb01Min(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, phiStates, psiStates);
+                } else {
+                    statesWithProbability01 = storm::utility::graph::performProb01Max(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, phiStates, psiStates);
+                }
+                storm::storage::BitVector statesWithProbability0 = std::move(statesWithProbability01.first);
+                storm::storage::BitVector statesWithProbability1 = std::move(statesWithProbability01.second);
+                storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
+                LOG4CPLUS_INFO(logger, "Found " << statesWithProbability0.getNumberOfSetBits() << " 'no' states.");
+                LOG4CPLUS_INFO(logger, "Found " << statesWithProbability1.getNumberOfSetBits() << " 'yes' states.");
+                LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
+                
+                // Create resulting vector.
+                std::vector<ValueType> result(numberOfStates);
+                
+                // Check whether we need to compute exact probabilities for some states.
+                if (qualitative) {
+                    // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1.
+                    storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, ValueType(0.5));
+                } else {
+                    if (!maybeStates.empty()) {
+                        // In this case we have have to compute the probabilities.
+                        
+                        // First, we can eliminate the rows and columns from the original transition probability matrix for states
+                        // whose probabilities are already known.
+                        storm::storage::SparseMatrix<ValueType> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, false);
+                        
+                        // Prepare the right-hand side of the equation system. For entry i this corresponds to
+                        // the accumulated probability of going from state i to some 'yes' state.
+                        std::vector<ValueType> b = transitionMatrix.getConstrainedRowGroupSumVector(maybeStates, statesWithProbability1);
+                        
+                        // Create vector for results for maybe states.
+                        std::vector<ValueType> x(maybeStates.getNumberOfSetBits());
+                        
+                        // Solve the corresponding system of equations.
+                        std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(submatrix);
+                        solver->solveEquationSystem(minimize, x, b);
+                        
+                        // Set values of resulting vector according to result.
+                        storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
+                    }
+                }
+                
+                // Set values of resulting vector that are known exactly.
+                storm::utility::vector::setVectorValues<ValueType>(result, statesWithProbability0, storm::utility::zero<ValueType>());
+                storm::utility::vector::setVectorValues<ValueType>(result, statesWithProbability1, storm::utility::one<ValueType>());
+                
+                return result;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMdpPrctlHelper<ValueType, RewardModelType>::computeInstantaneousRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, bool minimize, uint_fast64_t stepCount, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                // Only compute the result if the model has a state-based reward this->getModel().
+                STORM_LOG_THROW(rewardModel.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // Initialize result to state rewards of the this->getModel().
+                std::vector<ValueType> result(rewardModel.getStateRewardVector());
+                
+                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix);
+                solver->performMatrixVectorMultiplication(minimize, result, nullptr, stepCount);
+                
+                return result;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMdpPrctlHelper<ValueType, RewardModelType>::computeCumulativeRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, bool minimize, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                // Only compute the result if the model has at least one reward this->getModel().
+                STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // Compute the reward vector to add in each step based on the available reward models.
+                std::vector<ValueType> totalRewardVector = rewardModel.getTotalRewardVector(transitionMatrix);
+                
+                // Initialize result to either the state rewards of the model or the null vector.
+                std::vector<ValueType> result;
+                if (rewardModel.hasStateRewards()) {
+                    result = rewardModel.getStateRewardVector();
+                } else {
+                    result.resize(transitionMatrix.getRowCount());
+                }
+                
+                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix);
+                solver->performMatrixVectorMultiplication(minimize, result, &totalRewardVector, stepBound);
+                
+                return result;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMdpPrctlHelper<ValueType, RewardModelType>::computeReachabilityRewards(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                // Only compute the result if the model has at least one reward this->getModel().
+                STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // Determine which states have a reward of infinity by definition.
+                storm::storage::BitVector infinityStates;
+                storm::storage::BitVector trueStates(transitionMatrix.getRowCount(), true);
+                if (minimize) {
+                    infinityStates = std::move(storm::utility::graph::performProb1A(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, trueStates, targetStates));
+                } else {
+                    infinityStates = std::move(storm::utility::graph::performProb1E(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, trueStates, targetStates));
+                }
+                infinityStates.complement();
+                storm::storage::BitVector maybeStates = ~targetStates & ~infinityStates;
+                LOG4CPLUS_INFO(logger, "Found " << infinityStates.getNumberOfSetBits() << " 'infinity' states.");
+                LOG4CPLUS_INFO(logger, "Found " << targetStates.getNumberOfSetBits() << " 'target' states.");
+                LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
+                
+                // Create resulting vector.
+                std::vector<ValueType> result(transitionMatrix.getRowCount());
+                
+                // Check whether we need to compute exact rewards for some states.
+                if (qualitative) {
+                    LOG4CPLUS_INFO(logger, "The rewards for the initial states were determined in a preprocessing step. No exact rewards were computed.");
+                    // Set the values for all maybe-states to 1 to indicate that their reward values
+                    // are neither 0 nor infinity.
+                    storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, storm::utility::one<ValueType>());
+                } else {
+                    // In this case we have to compute the reward values for the remaining states.
+                    
+                    // We can eliminate the rows and columns from the original transition probability matrix for states
+                    // whose reward values are already known.
+                    storm::storage::SparseMatrix<ValueType> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, false);
+                    
+                    // Prepare the right-hand side of the equation system.
+                    std::vector<ValueType> b = rewardModel.getTotalRewardVector(submatrix.getRowCount(), transitionMatrix, maybeStates);
+                    
+                    // Create vector for results for maybe states.
+                    std::vector<ValueType> x(maybeStates.getNumberOfSetBits());
+                    
+                    // Solve the corresponding system of equations.
+                    std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(submatrix);
+                    solver->solveEquationSystem(minimize, x, b);
+                    
+                    
+                    // Set values of resulting vector according to result.
+                    storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
+                }
+                
+                // Set values of resulting vector that are known exactly.
+                storm::utility::vector::setVectorValues(result, targetStates, storm::utility::zero<ValueType>());
+                storm::utility::vector::setVectorValues(result, infinityStates, storm::utility::infinity<ValueType>());
+                
+                return result;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseMdpPrctlHelper<ValueType, RewardModelType>::computeLongRunAverage(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
+                // If there are no goal states, we avoid the computation and directly return zero.
+                uint_fast64_t numberOfStates = transitionMatrix.getRowCount();
+                if (psiStates.empty()) {
+                    return std::vector<ValueType>(numberOfStates, storm::utility::zero<ValueType>());
+                }
+                
+                // Likewise, if all bits are set, we can avoid the computation and set.
+                if ((~psiStates).empty()) {
+                    return std::vector<ValueType>(numberOfStates, storm::utility::one<ValueType>());
+                }
+                
+                // Start by decomposing the MDP into its MECs.
+                storm::storage::MaximalEndComponentDecomposition<double> mecDecomposition(transitionMatrix, backwardTransitions);
+                
+                // Get some data members for convenience.
+                std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = transitionMatrix.getRowGroupIndices();
+                ValueType zero = storm::utility::zero<ValueType>();
+                
+                //first calculate LRA for the Maximal End Components.
+                storm::storage::BitVector statesInMecs(numberOfStates);
+                std::vector<uint_fast64_t> stateToMecIndexMap(transitionMatrix.getColumnCount());
+                std::vector<ValueType> lraValuesForEndComponents(mecDecomposition.size(), zero);
+                
+                for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) {
+                    storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex];
+                    
+                    lraValuesForEndComponents[currentMecIndex] = computeLraForMaximalEndComponent(minimize, transitionMatrix, psiStates, mec);
+                    
+                    // Gather information for later use.
+                    for (auto const& stateChoicesPair : mec) {
+                        statesInMecs.set(stateChoicesPair.first);
+                        stateToMecIndexMap[stateChoicesPair.first] = currentMecIndex;
+                    }
+                }
+                
+                // For fast transition rewriting, we build some auxiliary data structures.
+                storm::storage::BitVector statesNotContainedInAnyMec = ~statesInMecs;
+                uint_fast64_t firstAuxiliaryStateIndex = statesNotContainedInAnyMec.getNumberOfSetBits();
+                uint_fast64_t lastStateNotInMecs = 0;
+                uint_fast64_t numberOfStatesNotInMecs = 0;
+                std::vector<uint_fast64_t> statesNotInMecsBeforeIndex;
+                statesNotInMecsBeforeIndex.reserve(numberOfStates);
+                for (auto state : statesNotContainedInAnyMec) {
+                    while (lastStateNotInMecs <= state) {
+                        statesNotInMecsBeforeIndex.push_back(numberOfStatesNotInMecs);
+                        ++lastStateNotInMecs;
+                    }
+                    ++numberOfStatesNotInMecs;
+                }
+                
+                // Finally, we are ready to create the SSP matrix and right-hand side of the SSP.
+                std::vector<ValueType> b;
+                typename storm::storage::SparseMatrixBuilder<ValueType> sspMatrixBuilder(0, 0, 0, false, true, numberOfStatesNotInMecs + mecDecomposition.size());
+                
+                // If the source state is not contained in any MEC, we copy its choices (and perform the necessary modifications).
+                uint_fast64_t currentChoice = 0;
+                for (auto state : statesNotContainedInAnyMec) {
+                    sspMatrixBuilder.newRowGroup(currentChoice);
+                    
+                    for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice, ++currentChoice) {
+                        std::vector<ValueType> auxiliaryStateToProbabilityMap(mecDecomposition.size());
+                        b.push_back(storm::utility::zero<ValueType>());
+                        
+                        for (auto element : transitionMatrix.getRow(choice)) {
+                            if (statesNotContainedInAnyMec.get(element.getColumn())) {
+                                // If the target state is not contained in an MEC, we can copy over the entry.
+                                sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue());
+                            } else {
+                                // If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector
+                                // so that we are able to write the cumulative probability to the MEC into the matrix.
+                                auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue();
+                            }
+                        }
+                        
+                        // Now insert all (cumulative) probability values that target an MEC.
+                        for (uint_fast64_t mecIndex = 0; mecIndex < auxiliaryStateToProbabilityMap.size(); ++mecIndex) {
+                            if (auxiliaryStateToProbabilityMap[mecIndex] != 0) {
+                                sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + mecIndex, auxiliaryStateToProbabilityMap[mecIndex]);
+                            }
+                        }
+                    }
+                }
+                
+                // Now we are ready to construct the choices for the auxiliary states.
+                for (uint_fast64_t mecIndex = 0; mecIndex < mecDecomposition.size(); ++mecIndex) {
+                    storm::storage::MaximalEndComponent const& mec = mecDecomposition[mecIndex];
+                    sspMatrixBuilder.newRowGroup(currentChoice);
+                    
+                    for (auto const& stateChoicesPair : mec) {
+                        uint_fast64_t state = stateChoicesPair.first;
+                        boost::container::flat_set<uint_fast64_t> const& choicesInMec = stateChoicesPair.second;
+                        
+                        for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) {
+                            std::vector<ValueType> auxiliaryStateToProbabilityMap(mecDecomposition.size());
+                            
+                            // If the choice is not contained in the MEC itself, we have to add a similar distribution to the auxiliary state.
+                            if (choicesInMec.find(choice) == choicesInMec.end()) {
+                                b.push_back(storm::utility::zero<ValueType>());
+                                
+                                for (auto element : transitionMatrix.getRow(choice)) {
+                                    if (statesNotContainedInAnyMec.get(element.getColumn())) {
+                                        // If the target state is not contained in an MEC, we can copy over the entry.
+                                        sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue());
+                                    } else {
+                                        // If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector
+                                        // so that we are able to write the cumulative probability to the MEC into the matrix.
+                                        auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue();
+                                    }
+                                }
+                                
+                                // Now insert all (cumulative) probability values that target an MEC.
+                                for (uint_fast64_t targetMecIndex = 0; targetMecIndex < auxiliaryStateToProbabilityMap.size(); ++targetMecIndex) {
+                                    if (auxiliaryStateToProbabilityMap[targetMecIndex] != 0) {
+                                        sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]);
+                                    }
+                                }
+                                
+                                ++currentChoice;
+                            }
+                        }
+                    }
+                    
+                    // For each auxiliary state, there is the option to achieve the reward value of the LRA associated with the MEC.
+                    ++currentChoice;
+                    b.push_back(lraValuesForEndComponents[mecIndex]);
+                }
+                
+                // Finalize the matrix and solve the corresponding system of equations.
+                storm::storage::SparseMatrix<ValueType> sspMatrix = sspMatrixBuilder.build(currentChoice);
+                
+                std::vector<ValueType> sspResult(numberOfStatesNotInMecs + mecDecomposition.size());
+                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(sspMatrix);
+                solver->solveEquationSystem(minimize, sspResult, b);
+                
+                // Prepare result vector.
+                std::vector<ValueType> result(numberOfStates, zero);
+                
+                // Set the values for states not contained in MECs.
+                storm::utility::vector::setVectorValues(result, statesNotContainedInAnyMec, sspResult);
+                
+                // Set the values for all states in MECs.
+                for (auto state : statesInMecs) {
+                    result[state] = sspResult[firstAuxiliaryStateIndex + stateToMecIndexMap[state]];
+                }
+                
+                return result;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            ValueType SparseMdpPrctlHelper<ValueType, RewardModelType>::computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::MaximalEndComponent const& mec) {
+                std::shared_ptr<storm::solver::LpSolver> solver = storm::utility::solver::getLpSolver("LRA for MEC");
+                solver->setModelSense(minimize ? storm::solver::LpSolver::ModelSense::Maximize : storm::solver::LpSolver::ModelSense::Minimize);
+                
+                //// First, we need to create the variables for the problem.
+                std::map<uint_fast64_t, storm::expressions::Variable> stateToVariableMap;
+                for (auto const& stateChoicesPair : mec) {
+                    std::string variableName = "h" + std::to_string(stateChoicesPair.first);
+                    stateToVariableMap[stateChoicesPair.first] = solver->addUnboundedContinuousVariable(variableName);
+                }
+                storm::expressions::Variable lambda = solver->addUnboundedContinuousVariable("L", 1);
+                solver->update();
+                
+                // Now we encode the problem as constraints.
+                for (auto const& stateChoicesPair : mec) {
+                    uint_fast64_t state = stateChoicesPair.first;
+                    
+                    // Now, based on the type of the state, create a suitable constraint.
+                    for (auto choice : stateChoicesPair.second) {
+                        storm::expressions::Expression constraint = -lambda;
+                        ValueType r = 0;
+                        
+                        for (auto element : transitionMatrix.getRow(choice)) {
+                            constraint = constraint + stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue());
+                            if (psiStates.get(element.getColumn())) {
+                                r += element.getValue();
+                            }
+                        }
+                        constraint = solver->getConstant(r) + constraint;
+                        
+                        if (minimize) {
+                            constraint = stateToVariableMap.at(state) <= constraint;
+                        } else {
+                            constraint = stateToVariableMap.at(state) >= constraint;
+                        }
+                        solver->addConstraint("state" + std::to_string(state) + "," + std::to_string(choice), constraint);
+                    }
+                }
+                
+                solver->optimize();
+                return solver->getContinuousValue(lambda);
+            }
+            
+            template class SparseMdpPrctlHelper<double>;
+        }
+    }
+}
diff --git a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.h b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.h
new file mode 100644
index 000000000..2c3aa99e4
--- /dev/null
+++ b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.h
@@ -0,0 +1,45 @@
+#ifndef STORM_MODELCHECKER_SPARSE_MDP_PRCTL_MODELCHECKER_HELPER_H_
+#define STORM_MODELCHECKER_SPARSE_MDP_PRCTL_MODELCHECKER_HELPER_H_
+
+#include <vector>
+
+#include "src/models/sparse/StandardRewardModel.h"
+
+#include "src/storage/SparseMatrix.h"
+#include "src/storage/BitVector.h"
+#include "src/storage/MaximalEndComponent.h"
+
+#include "src/utility/solver.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            
+            template <typename ValueType, typename RewardModelType = storm::models::sparse::StandardRewardModel<ValueType>>
+            class SparseMdpPrctlHelper {
+            public:
+                static std::vector<ValueType> computeBoundedUntilProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+
+                static std::vector<ValueType> computeNextProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+
+                static std::vector<ValueType> computeUntilProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+
+                static std::vector<ValueType> computeUntilProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& MinMaxLinearEquationSolverFactory, bool qualitative);
+                
+                static std::vector<ValueType> computeInstantaneousRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, bool minimize, uint_fast64_t stepCount, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeCumulativeRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, bool minimize, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeReachabilityRewards(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                
+                static std::vector<ValueType> computeLongRunAverage(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                
+            private:
+                static ValueType computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::MaximalEndComponent const& mec);
+            };
+            
+        }
+    }
+}
+
+#endif /* STORM_MODELCHECKER_SPARSE_MDP_PRCTL_MODELCHECKER_HELPER_H_ */
\ No newline at end of file
diff --git a/src/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp b/src/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp
new file mode 100644
index 000000000..66c523331
--- /dev/null
+++ b/src/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp
@@ -0,0 +1,195 @@
+#include "src/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.h"
+
+#include "src/storage/dd/DdType.h"
+#include "src/storage/dd/CuddAdd.h"
+#include "src/storage/dd/CuddBdd.h"
+#include "src/storage/dd/CuddOdd.h"
+
+#include "src/utility/graph.h"
+
+#include "src/exceptions/InvalidPropertyException.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+     
+            template<storm::dd::DdType DdType, typename ValueType>
+            storm::dd::Add<DdType> SymbolicDtmcPrctlHelper<DdType, ValueType>::computeUntilProbabilities(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+                // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
+                // probability 0 and 1 of satisfying the until-formula.
+                std::pair<storm::dd::Bdd<DdType>, storm::dd::Bdd<DdType>> statesWithProbability01 = storm::utility::graph::performProb01(model, transitionMatrix, phiStates, psiStates);
+                storm::dd::Bdd<DdType> maybeStates = !statesWithProbability01.first && !statesWithProbability01.second && model.getReachableStates();
+                
+                // Perform some logging.
+                STORM_LOG_INFO("Found " << statesWithProbability01.first.getNonZeroCount() << " 'no' states.");
+                STORM_LOG_INFO("Found " << statesWithProbability01.second.getNonZeroCount() << " 'yes' states.");
+                STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
+                
+                // Check whether we need to compute exact probabilities for some states.
+                if (qualitative) {
+                    // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1.
+                    return statesWithProbability01.second.toAdd() + maybeStates.toAdd() * model.getManager().getConstant(0.5);
+                } else {
+                    // If there are maybe states, we need to solve an equation system.
+                    if (!maybeStates.isZero()) {
+                        // Create the ODD for the translation between symbolic and explicit storage.
+                        storm::dd::Odd<DdType> odd(maybeStates);
+                        
+                        // Create the matrix and the vector for the equation system.
+                        storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                        
+                        // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                        // non-maybe states in the matrix.
+                        storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                        
+                        // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
+                        // maybe states.
+                        storm::dd::Add<DdType> prob1StatesAsColumn = statesWithProbability01.second.toAdd();
+                        prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs());
+                        storm::dd::Add<DdType> subvector = submatrix * prob1StatesAsColumn;
+                        subvector = subvector.sumAbstract(model.getColumnVariables());
+                        
+                        // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
+                        // for solving the equation system (i.e. compute (I-A)).
+                        submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                        submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
+                        
+                        // Solve the equation system.
+                        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                        storm::dd::Add<DdType> result = solver->solveEquationSystem(model.getManager().getConstant(0.5) * maybeStatesAdd, subvector);
+                        
+                        return statesWithProbability01.second.toAdd() + result;
+                    } else {
+                        return statesWithProbability01.second.toAdd();
+                    }
+                }
+            }
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            storm::dd::Add<DdType> SymbolicDtmcPrctlHelper<DdType, ValueType>::computeNextProbabilities(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates) {
+                storm::dd::Add<DdType> result = transitionMatrix * nextStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd();
+                return result.sumAbstract(model.getColumnVariables());
+            }
+
+            template<storm::dd::DdType DdType, typename ValueType>
+            storm::dd::Add<DdType> SymbolicDtmcPrctlHelper<DdType, ValueType>::computeBoundedUntilProbabilities(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+                // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
+                // probability 0 or 1 of satisfying the until-formula.
+                storm::dd::Bdd<DdType> statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(model, transitionMatrix.notZero(), phiStates, psiStates, stepBound);
+                storm::dd::Bdd<DdType> maybeStates = statesWithProbabilityGreater0 && !psiStates && model.getReachableStates();
+                
+                // If there are maybe states, we need to perform matrix-vector multiplications.
+                if (!maybeStates.isZero()) {
+                    // Create the ODD for the translation between symbolic and explicit storage.
+                    storm::dd::Odd<DdType> odd(maybeStates);
+                    
+                    // Create the matrix and the vector for the equation system.
+                    storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                    
+                    // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                    // non-maybe states in the matrix.
+                    storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                    
+                    // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
+                    // maybe states.
+                    storm::dd::Add<DdType> prob1StatesAsColumn = psiStates.toAdd().swapVariables(model.getRowColumnMetaVariablePairs());
+                    storm::dd::Add<DdType> subvector = (submatrix * prob1StatesAsColumn).sumAbstract(model.getColumnVariables());
+                    
+                    // Finally cut away all columns targeting non-maybe states.
+                    submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                    
+                    // Perform the matrix-vector multiplication.
+                    std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                    storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getManager().getAddZero(), &subvector, stepBound);
+                    
+                    return psiStates.toAdd() + result;
+                } else {
+                    return psiStates.toAdd();
+                }
+            }
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            storm::dd::Add<DdType> SymbolicDtmcPrctlHelper<DdType, ValueType>::computeCumulativeRewards(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+                // Only compute the result if the model has at least one reward this->getModel().
+                STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // Compute the reward vector to add in each step based on the available reward models.
+                storm::dd::Add<DdType> totalRewardVector = model.hasStateRewards() ? model.getStateRewardVector() : model.getManager().getAddZero();
+                if (model.hasTransitionRewards()) {
+                    totalRewardVector += (transitionMatrix * model.getTransitionRewardMatrix()).sumAbstract(model.getColumnVariables());
+                }
+                
+                // Perform the matrix-vector multiplication.
+                std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix, model.getReachableStates(), model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                return solver->performMatrixVectorMultiplication(model.getManager().getAddZero(), &totalRewardVector, stepBound);
+            }
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            storm::dd::Add<DdType> SymbolicDtmcPrctlHelper<DdType, ValueType>::computeInstantaneousRewards(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+                // Only compute the result if the model has at least one reward this->getModel().
+                STORM_LOG_THROW(model.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // Perform the matrix-vector multiplication.
+                std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix, model.getReachableStates(), model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                return solver->performMatrixVectorMultiplication(model.getStateRewardVector(), nullptr, stepBound);
+            }
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            storm::dd::Add<DdType> SymbolicDtmcPrctlHelper<DdType, ValueType>::computeReachabilityRewards(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+                
+                // Only compute the result if there is at least one reward model.
+                STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // Determine which states have a reward of infinity by definition.
+                storm::dd::Bdd<DdType> infinityStates = storm::utility::graph::performProb1(model, transitionMatrix.notZero(), model.getReachableStates(), targetStates);
+                infinityStates = !infinityStates && model.getReachableStates();
+                storm::dd::Bdd<DdType> maybeStates = (!targetStates && !infinityStates) && model.getReachableStates();
+                STORM_LOG_INFO("Found " << infinityStates.getNonZeroCount() << " 'infinity' states.");
+                STORM_LOG_INFO("Found " << targetStates.getNonZeroCount() << " 'target' states.");
+                STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
+                
+                // Check whether we need to compute exact rewards for some states.
+                if (qualitative) {
+                    // Set the values for all maybe-states to 1 to indicate that their reward values
+                    // are neither 0 nor infinity.
+                    return infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + maybeStates.toAdd() * model.getManager().getConstant(storm::utility::one<ValueType>());
+                } else {
+                    // If there are maybe states, we need to solve an equation system.
+                    if (!maybeStates.isZero()) {
+                        // Create the ODD for the translation between symbolic and explicit storage.
+                        storm::dd::Odd<DdType> odd(maybeStates);
+                        
+                        // Create the matrix and the vector for the equation system.
+                        storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                        
+                        // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                        // non-maybe states in the matrix.
+                        storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                        
+                        // Then compute the state reward vector to use in the computation.
+                        storm::dd::Add<DdType> subvector = stateRewardVector ? maybeStatesAdd * stateRewardVector.get() : model.getManager().getAddZero();
+                        if (transitionRewardMatrix) {
+                            subvector += (submatrix * transitionRewardMatrix.get()).sumAbstract(model.getColumnVariables());
+                        }
+                        
+                        // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
+                        // for solving the equation system (i.e. compute (I-A)).
+                        submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                        submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
+                        
+                        // Solve the equation system.
+                        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                        storm::dd::Add<DdType> result = solver->solveEquationSystem(model.getManager().getConstant(0.5) * maybeStatesAdd, subvector);
+                        
+                        return infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + result;
+                    } else {
+                        return infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>());
+                    }
+                }
+            }
+            
+            template class SymbolicDtmcPrctlHelper<storm::dd::DdType::CUDD, double>;
+            
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.h b/src/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.h
new file mode 100644
index 000000000..5f1519ab0
--- /dev/null
+++ b/src/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.h
@@ -0,0 +1,35 @@
+#ifndef STORM_MODELCHECKER_SPARSE_DTMC_PRCTL_MODELCHECKER_HELPER_H_
+#define STORM_MODELCHECKER_SPARSE_DTMC_PRCTL_MODELCHECKER_HELPER_H_
+
+#include "src/models/symbolic/Model.h"
+
+#include "src/storage/dd/Add.h"
+#include "src/storage/dd/Bdd.h"
+
+#include "src/utility/solver.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            class SymbolicDtmcPrctlHelper {
+            public:
+                static storm::dd::Add<DdType> computeBoundedUntilProbabilities(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+                
+                static storm::dd::Add<DdType> computeNextProbabilities(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates);
+                
+                static storm::dd::Add<DdType> computeUntilProbabilities(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+                
+                static storm::dd::Add<DdType> computeCumulativeRewards(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+                
+                static storm::dd::Add<DdType> computeInstantaneousRewards(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+                
+                static storm::dd::Add<DdType> computeReachabilityRewards(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+            };
+            
+        }
+    }
+}
+
+#endif /* STORM_MODELCHECKER_SPARSE_DTMC_PRCTL_MODELCHECKER_HELPER_H_ */
\ No newline at end of file
diff --git a/src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h b/src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/modelchecker/propositional/SparsePropositionalModelChecker.h b/src/modelchecker/propositional/SparsePropositionalModelChecker.h
index 5270f4b37..f24ece9ab 100644
--- a/src/modelchecker/propositional/SparsePropositionalModelChecker.h
+++ b/src/modelchecker/propositional/SparsePropositionalModelChecker.h
@@ -11,6 +11,9 @@ namespace storm {
         template<typename SparseModelType>
         class SparsePropositionalModelChecker : public AbstractModelChecker {
         public:
+            typedef typename SparseModelType::ValueType ValueType;
+            typedef typename SparseModelType::RewardModelType RewardModelType;
+            
             explicit SparsePropositionalModelChecker(SparseModelType const& model);
             
             // The implemented methods of the AbstractModelChecker interface.
diff --git a/src/models/sparse/Ctmc.h b/src/models/sparse/Ctmc.h
index 0712cdb75..728c645ae 100644
--- a/src/models/sparse/Ctmc.h
+++ b/src/models/sparse/Ctmc.h
@@ -1,9 +1,6 @@
 #ifndef STORM_MODELS_SPARSE_CTMC_H_
 #define STORM_MODELS_SPARSE_CTMC_H_
 
-#include <memory>
-#include <vector>
-
 #include "src/models/sparse/DeterministicModel.h"
 #include "src/utility/OsDetection.h"
 
diff --git a/src/settings/ArgumentBuilder.h b/src/settings/ArgumentBuilder.h
index 18f03334d..8239286c5 100644
--- a/src/settings/ArgumentBuilder.h
+++ b/src/settings/ArgumentBuilder.h
@@ -190,7 +190,7 @@ return *this; \
 			ArgumentBuilder(ArgumentType type, std::string const& name, std::string const& description) : hasBeenBuilt(false), type(type), name(name), description(description), isOptional(false), hasDefaultValue(false), defaultValue_String(), defaultValue_Integer(), defaultValue_UnsignedInteger(), defaultValue_Double(), defaultValue_Boolean(), userValidationFunctions_String(), userValidationFunctions_Integer(), userValidationFunctions_UnsignedInteger(), userValidationFunctions_Double(), userValidationFunctions_Boolean() {
 				// Intentionally left empty.
 			}
-            
+
             // A flag that stores whether an argument has been built using this builder.
 			bool hasBeenBuilt;
             
diff --git a/src/storage/MaximalEndComponentDecomposition.cpp b/src/storage/MaximalEndComponentDecomposition.cpp
index f32db1b00..cc809e03f 100644
--- a/src/storage/MaximalEndComponentDecomposition.cpp
+++ b/src/storage/MaximalEndComponentDecomposition.cpp
@@ -14,12 +14,17 @@ namespace storm {
         
         template<typename ValueType>
         MaximalEndComponentDecomposition<ValueType>::MaximalEndComponentDecomposition(storm::models::sparse::NondeterministicModel<ValueType> const& model) {
-            performMaximalEndComponentDecomposition(model, storm::storage::BitVector(model.getNumberOfStates(), true));
+            performMaximalEndComponentDecomposition(model.getTransitionMatrix(), model.getBackwardTransitions(), storm::storage::BitVector(model.getNumberOfStates(), true));
+        }
+
+        template<typename ValueType>
+        MaximalEndComponentDecomposition<ValueType>::MaximalEndComponentDecomposition(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions) {
+            performMaximalEndComponentDecomposition(transitionMatrix, backwardTransitions, storm::storage::BitVector(transitionMatrix.getRowGroupCount(), true));
         }
         
         template<typename ValueType>
         MaximalEndComponentDecomposition<ValueType>::MaximalEndComponentDecomposition(storm::models::sparse::NondeterministicModel<ValueType> const& model, storm::storage::BitVector const& subsystem) {
-            performMaximalEndComponentDecomposition(model, subsystem);
+            performMaximalEndComponentDecomposition(model.getTransitionMatrix(), model.getBackwardTransitions(), subsystem);
         }
         
         template<typename ValueType>
@@ -45,30 +50,28 @@ namespace storm {
         }
         
         template <typename ValueType>
-        void MaximalEndComponentDecomposition<ValueType>::performMaximalEndComponentDecomposition(storm::models::sparse::NondeterministicModel<ValueType> const& model, storm::storage::BitVector const& subsystem) {
-            // Get some references for convenient access.
-            storm::storage::SparseMatrix<ValueType> backwardTransitions = model.getBackwardTransitions();
-            std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
-            storm::storage::SparseMatrix<ValueType> const& transitionMatrix = model.getTransitionMatrix();
+        void MaximalEndComponentDecomposition<ValueType>::performMaximalEndComponentDecomposition(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> backwardTransitions, storm::storage::BitVector const& subsystem) {
+            // Get some data for convenient access.
+            uint_fast64_t numberOfStates = transitionMatrix.getRowGroupCount();
+            std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = transitionMatrix.getRowGroupIndices();
             
             // Initialize the maximal end component list to be the full state space.
             std::list<StateBlock> endComponentStateSets;
             endComponentStateSets.emplace_back(subsystem.begin(), subsystem.end());
-            storm::storage::BitVector statesToCheck(model.getNumberOfStates());
+            storm::storage::BitVector statesToCheck(numberOfStates);
             
             // The iterator used here should really be a const_iterator.
             // However, gcc 4.8 (and assorted libraries) does not provide an erase(const_iterator) method for std::list
             // but only an erase(iterator). This is in compliance with the c++11 draft N3337, which specifies the change
             // from iterator to const_iterator only for "set, multiset, map [and] multimap".
-            // FIXME: As soon as gcc provides an erase(const_iterator) method, change this iterator back to a const_iterator.
-            for (std::list<StateBlock>::iterator mecIterator = endComponentStateSets.begin(); mecIterator != endComponentStateSets.end();) {
+            for (std::list<StateBlock>::const_iterator mecIterator = endComponentStateSets.begin(); mecIterator != endComponentStateSets.end();) {
                 StateBlock const& mec = *mecIterator;
 
                 // Keep track of whether the MEC changed during this iteration.
                 bool mecChanged = false;
                 
                 // Get an SCC decomposition of the current MEC candidate.
-                StronglyConnectedComponentDecomposition<ValueType> sccs(model, mec, true);
+                StronglyConnectedComponentDecomposition<ValueType> sccs(transitionMatrix, mec, true);
                 
                 // We need to do another iteration in case we have either more than once SCC or the SCC is smaller than
                 // the MEC canditate itself.
@@ -79,7 +82,7 @@ namespace storm {
                     statesToCheck.set(scc.begin(), scc.end());
                     
                     while (!statesToCheck.empty()) {
-                        storm::storage::BitVector statesToRemove(model.getNumberOfStates());
+                        storm::storage::BitVector statesToRemove(numberOfStates);
                         
                         for (auto state : statesToCheck) {
                             bool keepStateInMEC = false;
@@ -132,7 +135,7 @@ namespace storm {
                         }
                     }
                     
-                    std::list<StateBlock>::iterator eraseIterator(mecIterator);
+                    std::list<StateBlock>::const_iterator eraseIterator(mecIterator);
                     ++mecIterator;
                     endComponentStateSets.erase(eraseIterator);
                 } else {
diff --git a/src/storage/MaximalEndComponentDecomposition.h b/src/storage/MaximalEndComponentDecomposition.h
index d2dfa077a..672163e52 100644
--- a/src/storage/MaximalEndComponentDecomposition.h
+++ b/src/storage/MaximalEndComponentDecomposition.h
@@ -26,6 +26,14 @@ namespace storm  {
              */
             MaximalEndComponentDecomposition(storm::models::sparse::NondeterministicModel<ValueType> const& model);
             
+            /*
+             * Creates an MEC decomposition of the given model (represented by a row-grouped matrix).
+             *
+             * @param transitionMatrix The transition relation of model to decompose into MECs.
+             * @param backwardTransition The reversed transition relation.
+             */
+            MaximalEndComponentDecomposition(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions);
+            
             /*!
              * Creates an MEC decomposition of the given subsystem in the given model.
              *
@@ -67,10 +75,11 @@ namespace storm  {
              * Performs the actual decomposition of the given subsystem in the given model into MECs. As a side-effect
              * this stores the MECs found in the current decomposition.
              *
-             * @param model The model whose subsystem to decompose into MECs.
+             * @param transitionMatrix The transition matrix representing the system whose subsystem to decompose into MECs.
+             * @param backwardTransitions The reversed transition relation.
              * @param subsystem The subsystem to decompose.
              */
-            void performMaximalEndComponentDecomposition(storm::models::sparse::NondeterministicModel<ValueType> const& model, storm::storage::BitVector const& subsystem);
+            void performMaximalEndComponentDecomposition(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> backwardTransitions, storm::storage::BitVector const& subsystem);
         };
     }
 }
diff --git a/src/storage/StronglyConnectedComponentDecomposition.cpp b/src/storage/StronglyConnectedComponentDecomposition.cpp
index 9d032d709..a03a22b63 100644
--- a/src/storage/StronglyConnectedComponentDecomposition.cpp
+++ b/src/storage/StronglyConnectedComponentDecomposition.cpp
@@ -25,6 +25,18 @@ namespace storm {
             performSccDecomposition(model.getTransitionMatrix(), subsystem, dropNaiveSccs, onlyBottomSccs);
         }
         
+        template <typename ValueType, typename RewardModelType>
+        StronglyConnectedComponentDecomposition<ValueType, RewardModelType>::StronglyConnectedComponentDecomposition(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, StateBlock const& block, bool dropNaiveSccs, bool onlyBottomSccs) {
+            storm::storage::BitVector subsystem(transitionMatrix.getRowGroupCount(), block.begin(), block.end());
+            performSccDecomposition(transitionMatrix, subsystem, dropNaiveSccs, onlyBottomSccs);
+        }
+        
+        
+        template <typename ValueType, typename RewardModelType>
+        StronglyConnectedComponentDecomposition<ValueType, RewardModelType>::StronglyConnectedComponentDecomposition(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, bool dropNaiveSccs, bool onlyBottomSccs) {
+            performSccDecomposition(transitionMatrix, storm::storage::BitVector(transitionMatrix.getRowGroupCount(), true), dropNaiveSccs, onlyBottomSccs);
+        }
+
         template <typename ValueType, typename RewardModelType>
         StronglyConnectedComponentDecomposition<ValueType, RewardModelType>::StronglyConnectedComponentDecomposition(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& subsystem, bool dropNaiveSccs, bool onlyBottomSccs) {
             performSccDecomposition(transitionMatrix, subsystem, dropNaiveSccs, onlyBottomSccs);
diff --git a/src/storage/StronglyConnectedComponentDecomposition.h b/src/storage/StronglyConnectedComponentDecomposition.h
index 9f1e3d3b1..052c35cbe 100644
--- a/src/storage/StronglyConnectedComponentDecomposition.h
+++ b/src/storage/StronglyConnectedComponentDecomposition.h
@@ -62,6 +62,30 @@ namespace storm {
              */
             StronglyConnectedComponentDecomposition(storm::models::sparse::Model<ValueType, RewardModelType> const& model, storm::storage::BitVector const& subsystem, bool dropNaiveSccs = false, bool onlyBottomSccs = false);
 
+            /*
+             * Creates an SCC decomposition of the given subsystem in the given system (whose transition relation is
+             * given by a sparse matrix).
+             *
+             * @param transitionMatrix The transition matrix of the system to decompose.
+             * @param block The block to decompose into SCCs.
+             * @param dropNaiveSccs A flag that indicates whether trivial SCCs (i.e. SCCs consisting of just one state
+             * without a self-loop) are to be kept in the decomposition.
+             * @param onlyBottomSccs If set to true, only bottom SCCs, i.e. SCCs in which all states have no way of
+             * leaving the SCC), are kept.
+             */
+            StronglyConnectedComponentDecomposition(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, StateBlock const& block, bool dropNaiveSccs = false, bool onlyBottomSccs = false);
+            
+            /*
+             * Creates an SCC decomposition of the given system (whose transition relation is given by a sparse matrix).
+             *
+             * @param transitionMatrix The transition matrix of the system to decompose.
+             * @param dropNaiveSccs A flag that indicates whether trivial SCCs (i.e. SCCs consisting of just one state
+             * without a self-loop) are to be kept in the decomposition.
+             * @param onlyBottomSccs If set to true, only bottom SCCs, i.e. SCCs in which all states have no way of
+             * leaving the SCC), are kept.
+             */
+            StronglyConnectedComponentDecomposition(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, bool dropNaiveSccs = false, bool onlyBottomSccs = false);
+            
             /*
              * Creates an SCC decomposition of the given subsystem in the given system (whose transition relation is 
              * given by a sparse matrix).