diff --git a/src/builder/DdPrismModelBuilder.cpp b/src/builder/DdPrismModelBuilder.cpp
index 461713e9c..408faf8eb 100644
--- a/src/builder/DdPrismModelBuilder.cpp
+++ b/src/builder/DdPrismModelBuilder.cpp
@@ -214,12 +214,18 @@ namespace storm {
             // Extract all the labels used in the formula.
             std::vector<std::shared_ptr<storm::logic::AtomicLabelFormula const>> atomicLabelFormulas = formula.getAtomicLabelFormulas();
             for (auto const& formula : atomicLabelFormulas) {
+                if (!labelsToBuild) {
+                    labelsToBuild = std::set<std::string>();
+                }
                 labelsToBuild.get().insert(formula.get()->getLabel());
             }
             
             // Extract all the expressions used in the formula.
             std::vector<std::shared_ptr<storm::logic::AtomicExpressionFormula const>> atomicExpressionFormulas = formula.getAtomicExpressionFormulas();
             for (auto const& formula : atomicExpressionFormulas) {
+                if (!expressionLabels) {
+                    expressionLabels = std::vector<storm::expressions::Expression>();
+                }
                 expressionLabels.get().push_back(formula.get()->getExpression());
             }
         }
@@ -240,6 +246,17 @@ namespace storm {
             }
         }
         
+        template <storm::dd::DdType Type>
+        struct DdPrismModelBuilder<Type>::SystemResult {
+            SystemResult(storm::dd::Add<Type> const& allTransitionsDd, DdPrismModelBuilder<Type>::ModuleDecisionDiagram const& globalModule, boost::optional<storm::dd::Add<Type>> const& stateActionDd) : allTransitionsDd(allTransitionsDd), globalModule(globalModule), stateActionDd(stateActionDd) {
+                // Intentionally left empty.
+            }
+            
+            storm::dd::Add<Type> allTransitionsDd;
+            typename DdPrismModelBuilder<Type>::ModuleDecisionDiagram globalModule;
+            boost::optional<storm::dd::Add<Type>> stateActionDd;
+        };
+        
         template <storm::dd::DdType Type>
         storm::dd::Add<Type> DdPrismModelBuilder<Type>::createUpdateDecisionDiagram(GenerationInformation& generationInfo, storm::prism::Module const& module, storm::dd::Add<Type> const& guard, storm::prism::Update const& update) {
             storm::dd::Add<Type> updateDd = generationInfo.manager->getAddOne();
@@ -645,7 +662,7 @@ namespace storm {
         }
         
         template <storm::dd::DdType Type>
-        std::tuple<storm::dd::Add<Type>, typename DdPrismModelBuilder<Type>::ModuleDecisionDiagram, boost::optional<storm::dd::Add<Type>>> DdPrismModelBuilder<Type>::createSystemDecisionDiagram(GenerationInformation& generationInfo) {
+        typename DdPrismModelBuilder<Type>::SystemResult DdPrismModelBuilder<Type>::createSystemDecisionDiagram(GenerationInformation& generationInfo) {
             // Create the initial offset mapping.
             std::map<uint_fast64_t, uint_fast64_t> synchronizingActionToOffsetMap;
             for (auto const& actionIndex : generationInfo.program.getSynchronizingActionIndices()) {
@@ -701,12 +718,17 @@ namespace storm {
             }
             
             storm::dd::Add<Type> result = createSystemFromModule(generationInfo, system);
-            
+
+            // For MDPs and DTMCs, we need a DD that maps states to the legal choices and states to the number of
+            // available choices, respectively. As it happens, the construction is the same in both cases.
+            boost::optional<storm::dd::Add<Type>> stateActionDd;
+            if (generationInfo.program.getModelType() == storm::prism::Program::ModelType::DTMC || generationInfo.program.getModelType() == storm::prism::Program::ModelType::MDP) {
+                stateActionDd = result.sumAbstract(generationInfo.columnMetaVariables);
+            }
+
             // For DTMCs, we normalize each row to 1 (to account for non-determinism).
-            boost::optional<storm::dd::Add<Type>> stateToChoiceCountMap;
             if (generationInfo.program.getModelType() == storm::prism::Program::ModelType::DTMC) {
-                stateToChoiceCountMap = result.sumAbstract(generationInfo.columnMetaVariables);
-                result = result / stateToChoiceCountMap.get();
+                result = result / stateActionDd.get();
             } else if (generationInfo.program.getModelType() == storm::prism::Program::ModelType::MDP) {
                 // For MDPs, we need to throw away the nondeterminism variables from the generation information that
                 // were never used.
@@ -716,11 +738,11 @@ namespace storm {
                 generationInfo.nondeterminismMetaVariables.resize(system.numberOfUsedNondeterminismVariables);
             }
             
-            return std::make_tuple(result, system, stateToChoiceCountMap);
+            return SystemResult(result, system, stateActionDd);
         }
         
         template <storm::dd::DdType Type>
-        storm::models::symbolic::StandardRewardModel<Type, double> DdPrismModelBuilder<Type>::createRewardModelDecisionDiagrams(GenerationInformation& generationInfo, storm::prism::RewardModel const& rewardModel, ModuleDecisionDiagram const& globalModule, storm::dd::Add<Type> const& transitionMatrix, storm::dd::Add<Type> const& reachableStatesAdd, boost::optional<storm::dd::Add<Type>> const& stateToChoiceCountMap) {
+        storm::models::symbolic::StandardRewardModel<Type, double> DdPrismModelBuilder<Type>::createRewardModelDecisionDiagrams(GenerationInformation& generationInfo, storm::prism::RewardModel const& rewardModel, ModuleDecisionDiagram const& globalModule, storm::dd::Add<Type> const& transitionMatrix, storm::dd::Add<Type> const& reachableStatesAdd, boost::optional<storm::dd::Add<Type>> const& stateActionDd) {
             
             // Start by creating the state reward vector.
             boost::optional<storm::dd::Add<Type>> stateRewards;
@@ -766,6 +788,14 @@ namespace storm {
                     }
                     
                     storm::dd::Add<Type> stateActionRewardDd = synchronization * states * rewards;
+
+                    // If we are building the state-action rewards for an MDP, we need to make sure that the encoding
+                    // of the nondeterminism is present in the reward vector, so we ne need to multiply it with the
+                    // legal state-actions.
+                    if (generationInfo.program.getModelType() == storm::prism::Program::ModelType::MDP) {
+                        STORM_LOG_THROW(static_cast<bool>(stateActionDd), storm::exceptions::InvalidStateException, "A DD that reflects the legal choices of each state was not build, but was expected to exist.");
+                        stateActionRewardDd *= stateActionDd.get();
+                    }
                     
                     // Perform some sanity checks.
                     STORM_LOG_WARN_COND(stateActionRewardDd.getMin() >= 0, "The reward model assigns negative rewards to some states.");
@@ -777,10 +807,9 @@ namespace storm {
                 
                 // Scale state-action rewards for DTMCs.
                 if (generationInfo.program.getModelType() == storm::prism::Program::ModelType::DTMC) {
-                    STORM_LOG_THROW(static_cast<bool>(stateToChoiceCountMap), storm::exceptions::InvalidStateException, "A DD that reflects the number of choices per state was not build, but was expected to exist.");
-                    stateActionRewards.get() /= stateToChoiceCountMap.get();
+                    STORM_LOG_THROW(static_cast<bool>(stateActionDd), storm::exceptions::InvalidStateException, "A DD that reflects the number of choices per state was not build, but was expected to exist.");
+                    stateActionRewards.get() /= stateActionDd.get();
                 }
-                stateActionRewards.get().exportToDot("rewards.dot");
             }
             
             // Then build the transition reward matrix.
@@ -827,8 +856,8 @@ namespace storm {
                 
                 // Scale transition rewards for DTMCs.
                 if (generationInfo.program.getModelType() == storm::prism::Program::ModelType::DTMC) {
-                    STORM_LOG_THROW(static_cast<bool>(stateToChoiceCountMap), storm::exceptions::InvalidStateException, "A DD that reflects the number of choices per state was not build, but was expected to exist.");
-                    stateActionRewards.get() /= stateToChoiceCountMap.get();
+                    STORM_LOG_THROW(static_cast<bool>(stateActionDd), storm::exceptions::InvalidStateException, "A DD that reflects the number of choices per state was not build, but was expected to exist.");
+                    transitionRewards.get() /= stateActionDd.get();
                 }
             }
             
@@ -869,10 +898,10 @@ namespace storm {
             // In particular, this creates the meta variables used to encode the model.
             GenerationInformation generationInfo(preparedProgram);
 
-            std::tuple<storm::dd::Add<Type>, ModuleDecisionDiagram, boost::optional<storm::dd::Add<Type>>> system = createSystemDecisionDiagram(generationInfo);
-            storm::dd::Add<Type> transitionMatrix = std::get<0>(system);
-            ModuleDecisionDiagram const& globalModule = std::get<1>(system);
-            boost::optional<storm::dd::Add<Type>> stateToChoiceCountMap = std::get<2>(system);
+            SystemResult system = createSystemDecisionDiagram(generationInfo);
+            storm::dd::Add<Type> transitionMatrix = system.allTransitionsDd;
+            ModuleDecisionDiagram const& globalModule = system.globalModule;
+            boost::optional<storm::dd::Add<Type>> stateActionDd = system.stateActionDd;
             
             // Cut the transitions and rewards to the reachable fragment of the state space.
             storm::dd::Bdd<Type> initialStates = createInitialStatesDecisionDiagram(generationInfo);
@@ -883,6 +912,9 @@ namespace storm {
             storm::dd::Bdd<Type> reachableStates = computeReachableStates(generationInfo, initialStates, transitionMatrixBdd);
             storm::dd::Add<Type> reachableStatesAdd = reachableStates.toAdd();
             transitionMatrix *= reachableStatesAdd;
+            if (stateActionDd) {
+                stateActionDd.get() *= reachableStatesAdd;
+            }
 
             // Detect deadlocks and 1) fix them if requested 2) throw an error otherwise.
             storm::dd::Bdd<Type> statesWithTransition = transitionMatrixBdd.existsAbstract(generationInfo.columnMetaVariables);
@@ -928,13 +960,8 @@ namespace storm {
             }
             
             std::unordered_map<std::string, storm::models::symbolic::StandardRewardModel<Type, double>> rewardModels;
-            if (options.buildAllRewardModels || !options.rewardModelsToBuild.empty()) {
-                for (auto const& rewardModel : preparedProgram.getRewardModels()) {
-                    if (options.buildAllRewardModels || options.rewardModelsToBuild.find(rewardModel.getName()) != options.rewardModelsToBuild.end()) {
-                        STORM_LOG_TRACE("Building reward structure.");
-                        rewardModels.emplace(rewardModel.getName(), createRewardModelDecisionDiagrams(generationInfo, rewardModel, globalModule, transitionMatrix, reachableStatesAdd, stateToChoiceCountMap));
-                    }
-                }
+            for (auto const& rewardModel : selectedRewardModels) {
+                rewardModels.emplace(rewardModel.get().getName(), createRewardModelDecisionDiagrams(generationInfo, rewardModel.get(), globalModule, transitionMatrix, reachableStatesAdd, stateActionDd));
             }
             
             // Build the labels that can be accessed as a shortcut.
diff --git a/src/builder/DdPrismModelBuilder.h b/src/builder/DdPrismModelBuilder.h
index 4380598c4..57ee60414 100644
--- a/src/builder/DdPrismModelBuilder.h
+++ b/src/builder/DdPrismModelBuilder.h
@@ -164,7 +164,12 @@ namespace storm {
             /*!
              * Structure to store all information required to generate the model from the program.
              */
-            class GenerationInformation; 
+            class GenerationInformation;
+            
+            /*!
+             * Structure to store the result of the system creation phase.
+             */
+            struct SystemResult;
         private:
             
             static storm::dd::Add<Type> encodeChoice(GenerationInformation& generationInfo, uint_fast64_t nondeterminismVariableOffset, uint_fast64_t numberOfBinaryVariables, int_fast64_t value);
@@ -189,9 +194,9 @@ namespace storm {
             
             static storm::dd::Add<Type> createSystemFromModule(GenerationInformation& generationInfo, ModuleDecisionDiagram const& module);
             
-            static storm::models::symbolic::StandardRewardModel<Type, double> createRewardModelDecisionDiagrams(GenerationInformation& generationInfo, storm::prism::RewardModel const& rewardModel, ModuleDecisionDiagram const& globalModule, storm::dd::Add<Type> const& transitionMatrix, storm::dd::Add<Type> const& reachableStatesAdd, boost::optional<storm::dd::Add<Type>> const& stateToChoiceCountMap);
+            static storm::models::symbolic::StandardRewardModel<Type, double> createRewardModelDecisionDiagrams(GenerationInformation& generationInfo, storm::prism::RewardModel const& rewardModel, ModuleDecisionDiagram const& globalModule, storm::dd::Add<Type> const& transitionMatrix, storm::dd::Add<Type> const& reachableStatesAdd, boost::optional<storm::dd::Add<Type>> const& stateActionDd);
             
-            static std::tuple<storm::dd::Add<Type>, ModuleDecisionDiagram, boost::optional<storm::dd::Add<Type>>> createSystemDecisionDiagram(GenerationInformation& generationInfo);
+            static SystemResult createSystemDecisionDiagram(GenerationInformation& generationInfo);
             
             static storm::dd::Bdd<Type> createInitialStatesDecisionDiagram(GenerationInformation& generationInfo);
 
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
index 653bf6ab7..aa4c620d6 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
@@ -45,7 +45,7 @@ namespace storm {
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();
             ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();
-            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper<ValueType>::computeBoundedUntilProbabilities(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *minMaxLinearEquationSolverFactory);
+            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper<ValueType>::computeBoundedUntilProbabilities(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *minMaxLinearEquationSolverFactory);
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
         }
         
diff --git a/src/modelchecker/prctl/SymbolicMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SymbolicMdpPrctlModelChecker.cpp
new file mode 100644
index 000000000..b5d7df0b8
--- /dev/null
+++ b/src/modelchecker/prctl/SymbolicMdpPrctlModelChecker.cpp
@@ -0,0 +1,91 @@
+#include "src/modelchecker/prctl/SymbolicMdpPrctlModelChecker.h"
+
+#include "src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h"
+
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+
+#include "src/utility/macros.h"
+#include "src/utility/graph.h"
+
+#include "src/settings/modules/GeneralSettings.h"
+
+#include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/InvalidPropertyException.h"
+
+namespace storm {
+    namespace modelchecker {
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicMdpPrctlModelChecker<DdType, ValueType>::SymbolicMdpPrctlModelChecker(storm::models::symbolic::Mdp<DdType> const& model, std::unique_ptr<storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory) : SymbolicPropositionalModelChecker<DdType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicMdpPrctlModelChecker<DdType, ValueType>::SymbolicMdpPrctlModelChecker(storm::models::symbolic::Mdp<DdType> const& model) : SymbolicPropositionalModelChecker<DdType>(model), linearEquationSolverFactory(new storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType>()) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        bool SymbolicMdpPrctlModelChecker<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> SymbolicMdpPrctlModelChecker<DdType, ValueType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidPropertyException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
+            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 storm::modelchecker::helper::SymbolicMdpPrctlHelper<DdType, ValueType>::computeUntilProbabilities(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicMdpPrctlModelChecker<DdType, ValueType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidPropertyException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
+            std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return storm::modelchecker::helper::SymbolicMdpPrctlHelper<DdType, ValueType>::computeNextProbabilities(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector());
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicMdpPrctlModelChecker<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.");
+            STORM_LOG_THROW(pathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidPropertyException, "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());
+            SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return storm::modelchecker::helper::SymbolicMdpPrctlHelper<DdType, ValueType>::computeBoundedUntilProbabilities(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicMdpPrctlModelChecker<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::InvalidPropertyException, "Formula needs to have a discrete time bound.");
+            return storm::modelchecker::helper::SymbolicMdpPrctlHelper<DdType, ValueType>::computeCumulativeRewards(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicMdpPrctlModelChecker<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::InvalidPropertyException, "Formula needs to have a discrete time bound.");
+            return storm::modelchecker::helper::SymbolicMdpPrctlHelper<DdType, ValueType>::computeInstantaneousRewards(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicMdpPrctlModelChecker<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::InvalidPropertyException, "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 storm::modelchecker::helper::SymbolicMdpPrctlHelper<DdType, ValueType>::computeReachabilityRewards(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), subResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::models::symbolic::Mdp<DdType> const& SymbolicMdpPrctlModelChecker<DdType, ValueType>::getModel() const {
+            return this->template getModelAs<storm::models::symbolic::Mdp<DdType>>();
+        }
+        
+        template class SymbolicMdpPrctlModelChecker<storm::dd::DdType::CUDD, double>;
+    }
+}
\ No newline at end of file
diff --git a/src/modelchecker/prctl/SymbolicMdpPrctlModelChecker.h b/src/modelchecker/prctl/SymbolicMdpPrctlModelChecker.h
new file mode 100644
index 000000000..bbef7ce0f
--- /dev/null
+++ b/src/modelchecker/prctl/SymbolicMdpPrctlModelChecker.h
@@ -0,0 +1,38 @@
+#ifndef STORM_MODELCHECKER_SYMBOLICMDPPRCTLMODELCHECKER_H_
+#define STORM_MODELCHECKER_SYMBOLICMDPPRCTLMODELCHECKER_H_
+
+#include "src/modelchecker/propositional/SymbolicPropositionalModelChecker.h"
+
+#include "src/models/symbolic/Mdp.h"
+
+#include "src/utility/solver.h"
+
+namespace storm {
+    namespace modelchecker {
+        template<storm::dd::DdType DdType, typename ValueType>
+        class SymbolicMdpPrctlModelChecker : public SymbolicPropositionalModelChecker<DdType> {
+            public:
+            explicit SymbolicMdpPrctlModelChecker(storm::models::symbolic::Mdp<DdType> const& model);
+            explicit SymbolicMdpPrctlModelChecker(storm::models::symbolic::Mdp<DdType> const& model, std::unique_ptr<storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory);
+            
+            // The implemented methods of the AbstractModelChecker interface.
+            virtual bool canHandle(storm::logic::Formula const& formula) const override;
+            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>()) 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:
+            // An object that is used for retrieving linear equation solvers.
+            std::unique_ptr<storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType>> linearEquationSolverFactory;
+        };
+        
+    } // namespace modelchecker
+} // namespace storm
+
+#endif /* STORM_MODELCHECKER_SYMBOLICMDPPRCTLMODELCHECKER_H_ */
diff --git a/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp b/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp
index 1cdc1f840..d62ca4dcc 100644
--- a/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp
+++ b/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp
@@ -210,8 +210,8 @@ namespace storm {
                         storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
                         
                         // Then compute the state reward vector to use in the computation.
-                        storm::dd::Add<DdType> subvector = rewardModel.getTotalRewardVector(submatrix, model.getColumnVariables());
-                        
+                        storm::dd::Add<DdType> subvector = rewardModel.getTotalRewardVector(maybeStatesAdd, submatrix, 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());
diff --git a/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp
index 46476d95b..044c5510b 100644
--- a/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp
+++ b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp
@@ -198,9 +198,9 @@ namespace storm {
                 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));
+                } else {
+                    infinityStates = storm::utility::graph::performProb1A(model, transitionMatrixBdd, model.getReachableStates(), targetStates, storm::utility::graph::performProbGreater0A(model, transitionMatrixBdd, model.getReachableStates(), targetStates));
                 }
                 infinityStates = !infinityStates && model.getReachableStates();
                 storm::dd::Bdd<DdType> maybeStates = (!targetStates && !infinityStates) && model.getReachableStates();
@@ -227,16 +227,17 @@ namespace storm {
                         storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
                         
                         // Then compute the state reward vector to use in the computation.
-                        storm::dd::Add<DdType> subvector = rewardModel.getTotalRewardVector(submatrix, model.getColumnVariables());
+                        storm::dd::Add<DdType> subvector = rewardModel.getTotalRewardVector(maybeStatesAdd, submatrix, 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);
+                        storm::dd::Add<DdType> stateActionAdd = (submatrix.notZero().existsAbstract(model.getColumnVariables()) || subvector.notZero()).toAdd();
+                        std::vector<uint_fast64_t> rowGroupSizes = stateActionAdd.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));
+                        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);
diff --git a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
index 461181343..c4f488581 100644
--- a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
+++ b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
@@ -179,27 +179,29 @@ namespace storm {
                     // 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 = totalStateRewardVectorGetter(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);
+                    if (!maybeStates.empty()) {
+                        // 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 = totalStateRewardVectorGetter(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.
diff --git a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
index 8203396aa..f65a834d1 100644
--- a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
+++ b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
@@ -26,6 +26,7 @@ namespace storm {
                 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);
@@ -159,9 +160,9 @@ namespace storm {
                 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));
+                    infinityStates = storm::utility::graph::performProb1E(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, trueStates, targetStates);
                 } else {
-                    infinityStates = std::move(storm::utility::graph::performProb1E(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, trueStates, targetStates));
+                    infinityStates = storm::utility::graph::performProb1A(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, trueStates, targetStates);
                 }
                 infinityStates.complement();
                 storm::storage::BitVector maybeStates = ~targetStates & ~infinityStates;
@@ -170,7 +171,7 @@ namespace storm {
                 LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
                 
                 // Create resulting vector.
-                std::vector<ValueType> result(transitionMatrix.getRowCount());
+                std::vector<ValueType> result(transitionMatrix.getRowCount(), storm::utility::zero<ValueType>());
                 
                 // Check whether we need to compute exact rewards for some states.
                 if (qualitative) {
@@ -179,29 +180,29 @@ namespace storm {
                     // 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);
+                    if (!maybeStates.empty()) {
+                        // 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;
@@ -210,7 +211,7 @@ namespace storm {
             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();
+                uint_fast64_t numberOfStates = transitionMatrix.getRowGroupCount();
                 if (psiStates.empty()) {
                     return std::vector<ValueType>(numberOfStates, storm::utility::zero<ValueType>());
                 }
diff --git a/src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp
index e69de29bb..6588e6f02 100644
--- a/src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp
+++ b/src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp
@@ -0,0 +1,201 @@
+#include "src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h"
+
+#include "src/storage/dd/CuddDdManager.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/models/symbolic/StandardRewardModel.h"
+
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+
+#include "src/exceptions/InvalidPropertyException.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace helper {
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            std::unique_ptr<CheckResult> SymbolicMdpPrctlHelper<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::SymbolicMinMaxLinearEquationSolverFactory<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;
+                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 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.
+                        submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                        
+                        // Now solve the resulting equation system.
+                        std::unique_ptr<storm::solver::SymbolicMinMaxLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getIllegalMask() && maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getNondeterminismVariables(), model.getRowColumnMetaVariablePairs());
+                        storm::dd::Add<DdType> result = solver->solveEquationSystem(minimize, model.getManager().getAddZero(), 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> SymbolicMdpPrctlHelper<DdType, ValueType>::computeNextProbabilities(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 std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), result.sumAbstract(model.getColumnVariables())));
+            }
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            std::unique_ptr<CheckResult> SymbolicMdpPrctlHelper<DdType, ValueType>::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::SymbolicMinMaxLinearEquationSolverFactory<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;
+                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 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());
+                    
+                    std::unique_ptr<storm::solver::SymbolicMinMaxLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getIllegalMask() && maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getNondeterminismVariables(), model.getRowColumnMetaVariablePairs());
+                    storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(minimize, 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> SymbolicMdpPrctlHelper<DdType, ValueType>::computeInstantaneousRewards(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+                // Only compute the result if the model has at least one reward this->getModel().
+                STORM_LOG_THROW(rewardModel.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+                
+                // Perform the matrix-vector multiplication.
+                std::unique_ptr<storm::solver::SymbolicMinMaxLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(model.getTransitionMatrix(), model.getReachableStates(), model.getIllegalMask(), model.getRowVariables(), model.getColumnVariables(), model.getNondeterminismVariables(), model.getRowColumnMetaVariablePairs());
+                storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(minimize, rewardModel.getStateRewardVector(), nullptr, stepBound);
+
+                return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), result));
+            }
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            std::unique_ptr<CheckResult> SymbolicMdpPrctlHelper<DdType, ValueType>::computeCumulativeRewards(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+                // 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.
+                storm::dd::Add<DdType> totalRewardVector = rewardModel.getTotalRewardVector(transitionMatrix, model.getColumnVariables());
+                
+                // Perform the matrix-vector multiplication.
+                std::unique_ptr<storm::solver::SymbolicMinMaxLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(model.getTransitionMatrix(), model.getReachableStates(), model.getIllegalMask(), model.getRowVariables(), model.getColumnVariables(), model.getNondeterminismVariables(), model.getRowColumnMetaVariablePairs());
+                storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(minimize, model.getManager().getAddZero(), &totalRewardVector, stepBound);
+                
+                return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), result));
+            }
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            std::unique_ptr<CheckResult> SymbolicMdpPrctlHelper<DdType, ValueType>::computeReachabilityRewards(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, RewardModelType const& rewardModel, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+                
+                // Only compute the result if there is at least one reward model.
+                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::dd::Bdd<DdType> infinityStates;
+                storm::dd::Bdd<DdType> transitionMatrixBdd = transitionMatrix.notZero();
+                if (minimize) {
+                    infinityStates = storm::utility::graph::performProb1E(model, transitionMatrixBdd, model.getReachableStates(), targetStates, storm::utility::graph::performProbGreater0E(model, transitionMatrixBdd, model.getReachableStates(), targetStates));
+                } else {
+                    infinityStates = storm::utility::graph::performProb1A(model, transitionMatrixBdd, model.getReachableStates(), targetStates, storm::utility::graph::performProbGreater0A(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 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 = rewardModel.getTotalRewardVector(maybeStatesAdd, submatrix, model.getColumnVariables());
+                        
+                        // Finally cut away all columns targeting non-maybe states.
+                        submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                        
+                        // Now solve the resulting equation system.
+                        std::unique_ptr<storm::solver::SymbolicMinMaxLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getIllegalMask() && maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getNondeterminismVariables(), model.getRowColumnMetaVariablePairs());
+                        storm::dd::Add<DdType> result = solver->solveEquationSystem(minimize, model.getManager().getAddZero(), 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>())));
+                    }
+                }
+            }
+            
+            template class SymbolicMdpPrctlHelper<storm::dd::DdType::CUDD, double>;
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h b/src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h
index e69de29bb..f76dd6915 100644
--- a/src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h
+++ b/src/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.h
@@ -0,0 +1,40 @@
+#ifndef STORM_MODELCHECKER_SYMBOLIC_MDP_PRCTL_MODELCHECKER_HELPER_H_
+#define STORM_MODELCHECKER_SYMBOLIC_MDP_PRCTL_MODELCHECKER_HELPER_H_
+
+#include "src/models/symbolic/NondeterministicModel.h"
+
+#include "src/storage/dd/Add.h"
+#include "src/storage/dd/Bdd.h"
+
+#include "src/utility/solver.h"
+
+namespace storm {
+    namespace modelchecker {
+        // Forward-declare result class.
+        class CheckResult;
+        
+        namespace helper {
+            
+            template<storm::dd::DdType DdType, typename ValueType>
+            class SymbolicMdpPrctlHelper {
+                public:
+                typedef typename storm::models::symbolic::NondeterministicModel<DdType>::RewardModelType RewardModelType;
+                
+                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::SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> 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::SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> computeCumulativeRewards(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> computeInstantaneousRewards(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+                
+                static std::unique_ptr<CheckResult> computeReachabilityRewards(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, RewardModelType const& rewardModel, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+            };
+            
+        }
+    }
+}
+
+#endif /* STORM_MODELCHECKER_SYMBOLIC_MDP_PRCTL_MODELCHECKER_HELPER_H_ */
diff --git a/src/models/symbolic/StandardRewardModel.cpp b/src/models/symbolic/StandardRewardModel.cpp
index 4762ff242..b98462fa1 100644
--- a/src/models/symbolic/StandardRewardModel.cpp
+++ b/src/models/symbolic/StandardRewardModel.cpp
@@ -82,6 +82,22 @@ namespace storm {
                 return this->optionalTransitionRewardMatrix;
             }
             
+            template <storm::dd::DdType Type, typename ValueType>
+            storm::dd::Add<Type> StandardRewardModel<Type, ValueType>::getTotalRewardVector(storm::dd::Add<Type> const& filterAdd, storm::dd::Add<Type> const& transitionMatrix, std::set<storm::expressions::Variable> const& columnVariables) const {
+                storm::dd::Add<Type> result = transitionMatrix.getDdManager()->getAddZero();
+                if (this->hasStateRewards()) {
+                    result += filterAdd * optionalStateRewardVector.get();
+                }
+                if (this->hasStateActionRewards()) {
+                    optionalStateActionRewardVector.get().exportToDot("statActRew.dot");
+                    result += filterAdd * optionalStateActionRewardVector.get();
+                }
+                if (this->hasTransitionRewards()) {
+                    result += (transitionMatrix * this->getTransitionRewardMatrix()).sumAbstract(columnVariables);
+                }
+                return result;
+            }
+            
             template <storm::dd::DdType Type, typename ValueType>
             storm::dd::Add<Type> StandardRewardModel<Type, ValueType>::getTotalRewardVector(storm::dd::Add<Type> const& transitionMatrix, std::set<storm::expressions::Variable> const& columnVariables) const {
                 storm::dd::Add<Type> result = transitionMatrix.getDdManager()->getAddZero();
diff --git a/src/models/symbolic/StandardRewardModel.h b/src/models/symbolic/StandardRewardModel.h
index 76abc9fdd..1eb1048d8 100644
--- a/src/models/symbolic/StandardRewardModel.h
+++ b/src/models/symbolic/StandardRewardModel.h
@@ -149,6 +149,18 @@ namespace storm {
                  */
                 storm::dd::Add<Type> getTotalRewardVector(storm::dd::Add<Type> const& transitionMatrix, std::set<storm::expressions::Variable> const& columnVariables) const;
                 
+                /*!
+                 * Creates a vector representing the complete reward vector based on the state-, state-action- and
+                 * transition-based rewards in the reward model. The state- and state-action rewards are filtered with
+                 * the given filter DD.
+                 *
+                 * @param filterAdd The DD used for filtering the rewards.
+                 * @param transitionMatrix The matrix that is used to weight the values of the transition reward matrix.
+                 * @param columnVariables The column variables of the model.
+                 * @return The full state-action reward vector.
+                 */
+                storm::dd::Add<Type> getTotalRewardVector(storm::dd::Add<Type> const& filterAdd, storm::dd::Add<Type> const& transitionMatrix, std::set<storm::expressions::Variable> const& columnVariables) const;
+                
                 /*!
                  * Creates a vector representing the complete reward vector based on the state-, state-action- and
                  * transition-based rewards in the reward model.
diff --git a/src/parser/DeterministicModelParser.cpp b/src/parser/DeterministicModelParser.cpp
index fd163c56a..1f2dc4ad9 100644
--- a/src/parser/DeterministicModelParser.cpp
+++ b/src/parser/DeterministicModelParser.cpp
@@ -38,11 +38,13 @@ namespace storm {
 			return result;
 		}
 
-        storm::models::sparse::Dtmc<double> DeterministicModelParser::parseDtmc(std::string const & transitionsFilename, std::string const & labelingFilename, std::string const & stateRewardFilename, std::string const & transitionRewardFilename) {
+        storm::models::sparse::Dtmc<double> DeterministicModelParser::parseDtmc(std::string const & transitionsFilename, std::string const & labelingFilename, std::string const& stateRewardFilename, std::string const& transitionRewardFilename) {
             DeterministicModelParser::Result parserResult(std::move(parseDeterministicModel(transitionsFilename, labelingFilename, stateRewardFilename, transitionRewardFilename)));
 
             std::unordered_map<std::string, storm::models::sparse::StandardRewardModel<double>> rewardModels;
-            rewardModels.insert(std::make_pair("", storm::models::sparse::StandardRewardModel<double>(parserResult.stateRewards, boost::optional<std::vector<double>>(), parserResult.transitionRewards)));
+            if (!stateRewardFilename.empty() || !transitionRewardFilename.empty()) {
+                rewardModels.insert(std::make_pair("", storm::models::sparse::StandardRewardModel<double>(parserResult.stateRewards, boost::optional<std::vector<double>>(), parserResult.transitionRewards)));
+            }
             return storm::models::sparse::Dtmc<double>(std::move(parserResult.transitionSystem), std::move(parserResult.labeling), std::move(rewardModels), boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
 		}
 
@@ -50,7 +52,9 @@ namespace storm {
             DeterministicModelParser::Result parserResult(std::move(parseDeterministicModel(transitionsFilename, labelingFilename, stateRewardFilename, transitionRewardFilename)));
 
             std::unordered_map<std::string, storm::models::sparse::StandardRewardModel<double>> rewardModels;
-            rewardModels.insert(std::make_pair("", storm::models::sparse::StandardRewardModel<double>(parserResult.stateRewards, boost::optional<std::vector<double>>(), parserResult.transitionRewards)));
+            if (!stateRewardFilename.empty() || !transitionRewardFilename.empty()) {
+                rewardModels.insert(std::make_pair("", storm::models::sparse::StandardRewardModel<double>(parserResult.stateRewards, boost::optional<std::vector<double>>(), parserResult.transitionRewards)));
+            }
             return storm::models::sparse::Ctmc<double>(std::move(parserResult.transitionSystem), std::move(parserResult.labeling), std::move(rewardModels), boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>());
 		}
 
diff --git a/src/settings/modules/GmmxxEquationSolverSettings.cpp b/src/settings/modules/GmmxxEquationSolverSettings.cpp
index dbaf8832a..525e369c7 100644
--- a/src/settings/modules/GmmxxEquationSolverSettings.cpp
+++ b/src/settings/modules/GmmxxEquationSolverSettings.cpp
@@ -105,7 +105,8 @@ namespace storm {
             }
             
             bool GmmxxEquationSolverSettings::check() const {
-                bool optionsSet = isLinearEquationSystemMethodSet() || isPreconditioningMethodSet() || isRestartIterationCountSet() | isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet();
+                // This list does not include the precision, because this option is shared with other modules.
+                bool optionsSet = isLinearEquationSystemMethodSet() || isPreconditioningMethodSet() || isRestartIterationCountSet() | isMaximalIterationCountSet() || isConvergenceCriterionSet();
                 
                 STORM_LOG_WARN_COND(storm::settings::generalSettings().getEquationSolver() == storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx || !optionsSet, "gmm++ is not selected as the equation solver, so setting options for gmm++ has no effect.");
                 
diff --git a/src/settings/modules/NativeEquationSolverSettings.cpp b/src/settings/modules/NativeEquationSolverSettings.cpp
index 5584308c5..a61291503 100644
--- a/src/settings/modules/NativeEquationSolverSettings.cpp
+++ b/src/settings/modules/NativeEquationSolverSettings.cpp
@@ -78,7 +78,8 @@ namespace storm {
             }
             
             bool NativeEquationSolverSettings::check() const {
-                bool optionSet = isLinearEquationSystemTechniqueSet() || isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet();
+                // This list does not include the precision, because this option is shared with other modules.
+                bool optionSet = isLinearEquationSystemTechniqueSet() || isMaximalIterationCountSet() || isConvergenceCriterionSet();
                 
                 STORM_LOG_WARN_COND(storm::settings::generalSettings().getEquationSolver() == storm::settings::modules::GeneralSettings::EquationSolver::Native || !optionSet, "Native is not selected as the equation solver, so setting options for native has no effect.");
                 
diff --git a/src/solver/SymbolicLinearEquationSolver.h b/src/solver/SymbolicLinearEquationSolver.h
index 007f28df6..7c2b24ba9 100644
--- a/src/solver/SymbolicLinearEquationSolver.h
+++ b/src/solver/SymbolicLinearEquationSolver.h
@@ -59,7 +59,7 @@ namespace storm {
              * The solution of the set of linear equations will be written to the vector x. Note that the matrix A has
              * to be given upon construction time of the solver object.
              *
-             * @param x The solution vector that has to be computed. Its length must be equal to the number of rows of A.
+             * @param x The initial guess for the solution vector. Its length must be equal to the number of rows of A.
              * @param b The right-hand side of the equation system. Its length must be equal to the number of rows of A.
              * @return The solution of the equation system.
              */
diff --git a/src/solver/SymbolicMinMaxLinearEquationSolver.cpp b/src/solver/SymbolicMinMaxLinearEquationSolver.cpp
new file mode 100644
index 000000000..e7df8cbaa
--- /dev/null
+++ b/src/solver/SymbolicMinMaxLinearEquationSolver.cpp
@@ -0,0 +1,96 @@
+#include "src/solver/SymbolicMinMaxLinearEquationSolver.h"
+
+#include "src/storage/dd/CuddDdManager.h"
+#include "src/storage/dd/CuddAdd.h"
+
+#include "src/storage/dd/Add.h"
+
+#include "src/utility/constants.h"
+
+#include "src/settings/SettingsManager.h"
+#include "src/settings/modules/NativeEquationSolverSettings.h"
+
+namespace storm {
+    namespace solver {
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::SymbolicMinMaxLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> const& allRows, storm::dd::Bdd<DdType> const& illegalMask, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& choiceVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : A(A), allRows(allRows), illegalMaskAdd(illegalMask.toAdd() * A.getDdManager()->getConstant(storm::utility::infinity<ValueType>())), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), choiceVariables(choiceVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), relative(relative) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::SymbolicMinMaxLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> const& allRows, storm::dd::Bdd<DdType> const& illegalMask, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& choiceVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) : A(A), allRows(allRows), illegalMaskAdd(illegalMask.toAdd() * A.getDdManager()->getConstant(storm::utility::infinity<ValueType>())), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), choiceVariables(choiceVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs) {
+            // Get the settings object to customize solving.
+            storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
+            
+            // Get appropriate settings.
+            maximalNumberOfIterations = settings.getMaximalIterationCount();
+            precision = settings.getPrecision();
+            relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType>  SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationSystem(bool minimize, storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const& b) const {
+            // Set up the environment.
+            storm::dd::Add<DdType> xCopy = x;
+            uint_fast64_t iterations = 0;
+            bool converged = false;
+            
+            while (!converged && iterations < maximalNumberOfIterations) {
+                // Compute tmp = A * x + b
+                storm::dd::Add<DdType> xCopyAsColumn = xCopy.swapVariables(this->rowColumnMetaVariablePairs);
+                storm::dd::Add<DdType> tmp = this->A.multiplyMatrix(xCopyAsColumn, this->columnMetaVariables);
+                tmp += b;
+                
+                if (minimize) {
+                    // This is a hack and only here because of the lack of a suitable minAbstract/maxAbstract function
+                    // that can properly deal with a restriction of the choices.
+                    tmp += illegalMaskAdd;
+                    tmp = tmp.minAbstract(this->choiceVariables);
+                } else {
+                    tmp = tmp.maxAbstract(this->choiceVariables);
+                }
+                
+                // Now check if the process already converged within our precision.
+                converged = xCopy.equalModuloPrecision(tmp, precision, relative);
+                
+                // If the method did not converge yet, we prepare the x vector for the next iteration.
+                if (!converged) {
+                    xCopy = tmp;
+                }
+                
+                ++iterations;
+            }
+                        
+            return xCopy;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType> SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::performMatrixVectorMultiplication(bool minimize, storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const* b, uint_fast64_t n) const {
+            storm::dd::Add<DdType> xCopy = x;
+            
+            // Perform matrix-vector multiplication while the bound is met.
+            for (uint_fast64_t i = 0; i < n; ++i) {
+                xCopy = xCopy.swapVariables(this->rowColumnMetaVariablePairs);
+                xCopy = this->A.multiplyMatrix(xCopy, this->columnMetaVariables);
+                if (b != nullptr) {
+                    xCopy += *b;
+                }
+                
+                if (minimize) {
+                    // This is a hack and only here because of the lack of a suitable minAbstract/maxAbstract function
+                    // that can properly deal with a restriction of the choices.
+                    xCopy += illegalMaskAdd;
+                    xCopy = xCopy.minAbstract(this->choiceVariables);
+                } else {
+                    xCopy = xCopy.maxAbstract(this->choiceVariables);
+                }
+            }
+            
+            return xCopy;
+        }
+        
+        template class SymbolicMinMaxLinearEquationSolver<storm::dd::DdType::CUDD, double>;
+        
+    }
+}
diff --git a/src/solver/SymbolicMinMaxLinearEquationSolver.h b/src/solver/SymbolicMinMaxLinearEquationSolver.h
new file mode 100644
index 000000000..ab3cd45e9
--- /dev/null
+++ b/src/solver/SymbolicMinMaxLinearEquationSolver.h
@@ -0,0 +1,126 @@
+#ifndef STORM_SOLVER_SYMBOLICMINMAXLINEAREQUATIONSOLVER_H_
+#define STORM_SOLVER_SYMBOLICMINMAXLINEAREQUATIONSOLVER_H_
+
+#include <memory>
+#include <set>
+#include <vector>
+#include <boost/variant.hpp>
+
+#include "src/storage/expressions/Variable.h"
+#include "src/storage/dd/DdType.h"
+
+#include "src/storage/dd/CuddAdd.h"
+
+namespace storm {
+    namespace dd {
+        template<storm::dd::DdType T> class Add;
+        template<storm::dd::DdType T> class Bdd;
+    }
+    
+    namespace solver {
+        /*!
+         * An interface that represents an abstract symbolic linear equation solver. In addition to solving a system of
+         * linear equations, the functionality to repeatedly multiply a matrix with a given vector is provided.
+         */
+        template<storm::dd::DdType DdType, typename ValueType>
+        class SymbolicMinMaxLinearEquationSolver {
+        public:
+            /*!
+             * Constructs a symbolic min/max linear equation solver with the given meta variable sets and pairs.
+             *
+             * @param A The matrix defining the coefficients of the linear equation system.
+             * @param diagonal An ADD characterizing the elements on the diagonal of the matrix.
+             * @param allRows A BDD characterizing all rows of the equation system.
+             * @param illegalMask A mask that characterizes all illegal choices (that are therefore not to be taken).
+             * @param rowMetaVariables The meta variables used to encode the rows of the matrix.
+             * @param columnMetaVariables The meta variables used to encode the columns of the matrix.
+             * @param choiceVariables The variables encoding the choices of each row group.
+             * @param rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta
+             * variables.
+             */
+            SymbolicMinMaxLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> const& allRows, storm::dd::Bdd<DdType> const& illegalMask, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& choiceVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
+            
+            /*!
+             * Constructs a symbolic linear equation solver with the given meta variable sets and pairs.
+             *
+             * @param A The matrix defining the coefficients of the linear equation system.
+             * @param allRows A BDD characterizing all rows of the equation system.
+             * @param illegalMask A mask that characterizes all illegal choices (that are therefore not to be taken).
+             * @param rowMetaVariables The meta variables used to encode the rows of the matrix.
+             * @param columnMetaVariables The meta variables used to encode the columns of the matrix.
+             * @param choiceVariables The variables encoding the choices of each row group.
+             * @param rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta
+             * variables.
+             * @param precision The precision to achieve.
+             * @param maximalNumberOfIterations The maximal number of iterations to perform when solving a linear
+             * equation system iteratively.
+             * @param relative Sets whether or not to use a relativ stopping criterion rather than an absolute one.
+             */
+            SymbolicMinMaxLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> const& allRows, storm::dd::Bdd<DdType> const& illegalMask, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& choiceVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, double precision, uint_fast64_t maximalNumberOfIterations, bool relative);
+            
+            /*!
+             * Solves the equation system A*x = b. The matrix A is required to be square and have a unique solution.
+             * The solution of the set of linear equations will be written to the vector x. Note that the matrix A has
+             * to be given upon construction time of the solver object.
+             *
+             * @param minimize If set, all the value of a group of rows is the taken as the minimum over all rows and as
+             * the maximum otherwise.
+             * @param x The initual guess for the solution vector. Its length must be equal to the number of row
+             * groups of A.
+             * @param b The right-hand side of the equation system. Its length must be equal to the number of row groups
+             * of A.
+             * @return The solution of the equation system.
+             */
+            virtual storm::dd::Add<DdType> solveEquationSystem(bool minimize, storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const& b) const;
+            
+            /*!
+             * Performs repeated matrix-vector multiplication, using x[0] = x and x[i + 1] = A*x[i] + b. After
+             * performing the necessary multiplications, the result is written to the input vector x. Note that the
+             * matrix A has to be given upon construction time of the solver object.
+             *
+             * @param minimize If set, all the value of a group of rows is the taken as the minimum over all rows and as
+             * the maximum otherwise.
+             * @param x The initial vector with which to perform matrix-vector multiplication. Its length must be equal
+             * to the number of row groups of A.
+             * @param b If non-null, this vector is added after each multiplication. If given, its length must be equal
+             * to the number of row groups of A.
+             * @return The solution of the equation system.
+             */
+            virtual storm::dd::Add<DdType> performMatrixVectorMultiplication(bool minimize, storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const* b = nullptr, uint_fast64_t n = 1) const;
+            
+        protected:
+            // The matrix defining the coefficients of the linear equation system.
+            storm::dd::Add<DdType> const& A;
+            
+            // A BDD characterizing all rows of the equation system.
+            storm::dd::Bdd<DdType> const& allRows;
+            
+            // An ADD characterizing the illegal choices.
+            storm::dd::Add<DdType> illegalMaskAdd;
+            
+            // The row variables.
+            std::set<storm::expressions::Variable> rowMetaVariables;
+            
+            // The column variables.
+            std::set<storm::expressions::Variable> columnMetaVariables;
+            
+            // The choice variables
+            std::set<storm::expressions::Variable> const& choiceVariables;
+            
+            // The pairs of meta variables used for renaming.
+            std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs;
+            
+            // The precision to achive.
+            double precision;
+            
+            // The maximal number of iterations to perform.
+            uint_fast64_t maximalNumberOfIterations;
+            
+            // A flag indicating whether a relative or an absolute stopping criterion is to be used.
+            bool relative;
+        };
+        
+    } // namespace solver
+} // namespace storm
+
+#endif /* STORM_SOLVER_SYMBOLICMINMAXLINEAREQUATIONSOLVER_H_ */
diff --git a/src/storage/dd/CuddAdd.cpp b/src/storage/dd/CuddAdd.cpp
index 633128dd0..af8eea9ca 100644
--- a/src/storage/dd/CuddAdd.cpp
+++ b/src/storage/dd/CuddAdd.cpp
@@ -753,8 +753,7 @@ namespace storm {
             // Create the explicit vector we need to fill later.
             std::vector<double> explicitVector(rowGroupIndices.back());
             
-            // Next, we split the matrix into one for each group. This only works if the group variables are at the very
-            // top.
+            // Next, we split the matrix into one for each group. This only works if the group variables are at the very top.
             std::vector<std::pair<Add<DdType::CUDD>, Add<DdType::CUDD>>> groups;
             splitGroupsRec(this->getCuddDdNode(), vector.getCuddDdNode(), groups, ddGroupVariableIndices, 0, ddGroupVariableIndices.size(), rowAndColumnMetaVariables, rowMetaVariables);
             
diff --git a/src/storage/dd/CuddOdd.cpp b/src/storage/dd/CuddOdd.cpp
index d851ffaa0..1b68dad48 100644
--- a/src/storage/dd/CuddOdd.cpp
+++ b/src/storage/dd/CuddOdd.cpp
@@ -1,7 +1,9 @@
 #include "src/storage/dd/CuddOdd.h"
 
+#include <fstream>
 #include <algorithm>
 #include <boost/functional/hash.hpp>
+#include <boost/algorithm/string/join.hpp>
 
 #include "src/exceptions/InvalidArgumentException.h"
 #include "src/utility/macros.h"
@@ -95,6 +97,14 @@ namespace storm {
             }
         }
         
+        uint_fast64_t Odd<DdType::CUDD>::getHeight() const {
+            if (this->isTerminalNode()) {
+                return 1;
+            } else {
+                return 1 + this->getElseSuccessor().getHeight();
+            }
+        }
+        
         bool Odd<DdType::CUDD>::isTerminalNode() const {
             return this->elseNode == nullptr && this->thenNode == nullptr;
         }
@@ -163,6 +173,63 @@ namespace storm {
             }
         }
         
+        void Odd<DdType::CUDD>::exportToDot(std::string const& filename) const {
+            std::ofstream dotFile;
+            dotFile.open (filename);
+            
+            // Print header.
+            dotFile << "digraph \"ODD\" {" << std::endl << "center=true;" << std::endl << "edge [dir = none];" << std::endl;
+
+            // Print levels as ranks.
+            dotFile << "{ node [shape = plaintext];" << std::endl << "edge [style = invis];" << std::endl;
+            std::vector<std::string> levelNames;
+            for (uint_fast64_t level = 0; level < this->getHeight(); ++level) {
+                levelNames.push_back("\"" + std::to_string(level) + "\"");
+            }
+            dotFile << boost::join(levelNames, " -> ") << ";";
+            dotFile << "}" << std::endl;
+            
+            std::map<uint_fast64_t, std::vector<std::reference_wrapper<storm::dd::Odd<DdType::CUDD> const>>> levelToOddNodesMap;
+            this->addToLevelToOddNodesMap(levelToOddNodesMap);
+            
+            for (auto const& levelNodes : levelToOddNodesMap) {
+                dotFile << "{ rank = same; \"" << levelNodes.first << "\"" << std::endl;;
+                for (auto const& node : levelNodes.second) {
+                    dotFile << "\"" << &node.get() << "\";" << std::endl;
+                }
+                dotFile << "}" << std::endl;
+            }
+            
+            std::set<storm::dd::Odd<DdType::CUDD> const*> printedNodes;
+            for (auto const& levelNodes : levelToOddNodesMap) {
+                for (auto const& node : levelNodes.second) {
+                    if (printedNodes.find(&node.get()) != printedNodes.end()) {
+                        continue;
+                    } else {
+                        printedNodes.insert(&node.get());
+                    }
+                    
+                    dotFile << "\"" << &node.get() << "\" [label=\"" << levelNodes.first << "\"];" << std::endl;
+                    if (!node.get().isTerminalNode()) {
+                        dotFile << "\"" << &node.get() << "\" -> \"" << &node.get().getElseSuccessor() << "\" [style=dashed, label=\"0\"];" << std::endl;
+                        dotFile << "\"" << &node.get() << "\" -> \"" << &node.get().getThenSuccessor() << "\" [style=solid, label=\"" << node.get().getElseOffset() << "\"];" << std::endl;
+                    }
+                }
+            }
+            
+            dotFile << "}" << std::endl;
+            
+            dotFile.close();
+        }
+        
+        void Odd<DdType::CUDD>::addToLevelToOddNodesMap(std::map<uint_fast64_t, std::vector<std::reference_wrapper<storm::dd::Odd<DdType::CUDD> const>>>& levelToOddNodesMap, uint_fast64_t level) const {
+            levelToOddNodesMap[level].push_back(*this);
+            if (!this->isTerminalNode()) {
+                this->getElseSuccessor().addToLevelToOddNodesMap(levelToOddNodesMap, level + 1);
+                this->getThenSuccessor().addToLevelToOddNodesMap(levelToOddNodesMap, level + 1);
+            }
+        }
+        
         std::shared_ptr<Odd<DdType::CUDD>> Odd<DdType::CUDD>::buildOddFromAddRec(DdNode* dd, Cudd const& manager, uint_fast64_t currentLevel, uint_fast64_t maxLevel, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<std::unordered_map<DdNode*, std::shared_ptr<Odd<DdType::CUDD>>>>& uniqueTableForLevels) {
             // Check whether the ODD for this node has already been computed (for this level) and if so, return this instead.
             auto const& iterator = uniqueTableForLevels[currentLevel].find(dd);
diff --git a/src/storage/dd/CuddOdd.h b/src/storage/dd/CuddOdd.h
index e00bdd2c7..518213c81 100644
--- a/src/storage/dd/CuddOdd.h
+++ b/src/storage/dd/CuddOdd.h
@@ -96,6 +96,13 @@ namespace storm {
              */
             uint_fast64_t getNodeCount() const;
             
+            /*!
+             * Retrieves the height of the ODD.
+             *
+             * @return The height of the ODD.
+             */
+            uint_fast64_t getHeight() const;
+            
             /*!
              * Checks whether the given ODD node is a terminal node, i.e. has no successors.
              *
@@ -122,6 +129,13 @@ namespace storm {
              */
             void expandExplicitVector(storm::dd::Odd<DdType::CUDD> const& newOdd, std::vector<double> const& oldValues, std::vector<double>& newValues) const;
             
+            /*!
+             * Exports the ODD in the dot format to the given file.
+             *
+             * @param filename The name of the file to which to write the dot output.
+             */
+            void exportToDot(std::string const& filename) const;
+            
         private:
             // Declare a hash functor that is used for the unique tables in the construction process.
             class HashFunctor {
@@ -129,6 +143,14 @@ namespace storm {
                 std::size_t operator()(std::pair<DdNode*, bool> const& key) const;
             };
             
+            /*!
+             * Adds all nodes below the current one to the given mapping.
+             *
+             * @param levelToOddNodesMap A mapping of the level to the ODD node.
+             * @param The level of the current node.
+             */
+            void addToLevelToOddNodesMap(std::map<uint_fast64_t, std::vector<std::reference_wrapper<storm::dd::Odd<DdType::CUDD> const>>>& levelToOddNodesMap, uint_fast64_t level = 0) const;
+            
             /*!
              * Constructs an offset-labeled DD with the given topmost DD node, else- and then-successor.
              *
diff --git a/src/storage/prism/RewardModel.cpp b/src/storage/prism/RewardModel.cpp
index 9b1213f2b..3f38f7bd8 100644
--- a/src/storage/prism/RewardModel.cpp
+++ b/src/storage/prism/RewardModel.cpp
@@ -86,7 +86,7 @@ namespace storm {
             if (rewardModel.getName() != "") {
                 std::cout << " \"" << rewardModel.getName() << "\"";
             }
-            std::cout << std::endl;
+            stream << std::endl;
             for (auto const& reward : rewardModel.getStateRewards()) {
                 stream << reward << std::endl;
             }
diff --git a/src/storage/prism/Update.cpp b/src/storage/prism/Update.cpp
index b97fb6cc8..147e9b42e 100644
--- a/src/storage/prism/Update.cpp
+++ b/src/storage/prism/Update.cpp
@@ -56,7 +56,7 @@ namespace storm {
         std::ostream& operator<<(std::ostream& stream, Update const& update) {
             stream << update.getLikelihoodExpression() << " : ";
             if (update.getAssignments().empty()) {
-                stream << " true ";
+                stream << "true";
             } else {
                 std::vector<std::string> assignmentsAsStrings;
                 std::for_each(update.getAssignments().cbegin(), update.getAssignments().cend(), [&assignmentsAsStrings] (Assignment const& assignment) { std::stringstream stream; stream << assignment; assignmentsAsStrings.push_back(stream.str()); });
diff --git a/src/utility/cli.h b/src/utility/cli.h
index a47b12a50..73f18315e 100644
--- a/src/utility/cli.h
+++ b/src/utility/cli.h
@@ -48,12 +48,14 @@
 
 // Headers for model checking.
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
-#include "src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h"
 #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
-#include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
 #include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h"
-#include "src/modelchecker/csl/HybridCtmcCslModelChecker.h"
 #include "src/modelchecker/prctl/HybridMdpPrctlModelChecker.h"
+#include "src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h"
+#include "src/modelchecker/prctl/SymbolicMdpPrctlModelChecker.h"
+#include "src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h"
+#include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
+#include "src/modelchecker/csl/HybridCtmcCslModelChecker.h"
 #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
 #include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
 
@@ -289,7 +291,7 @@ namespace storm {
 #endif
             
             template<storm::dd::DdType DdType>
-            void verifySymbolicModel(boost::optional<storm::prism::Program> const& program, std::shared_ptr<storm::models::symbolic::Model<DdType>> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
+            void verifySymbolicModelWithHybridEngine(boost::optional<storm::prism::Program> const& program, std::shared_ptr<storm::models::symbolic::Model<DdType>> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
                 for (auto const& formula : formulas) {
                     std::cout << std::endl << "Model checking property: " << *formula << " ...";
                     std::unique_ptr<storm::modelchecker::CheckResult> result;
@@ -326,6 +328,38 @@ namespace storm {
                 }
             }
             
+            template<storm::dd::DdType DdType>
+            void verifySymbolicModelWithSymbolicEngine(boost::optional<storm::prism::Program> const& program, std::shared_ptr<storm::models::symbolic::Model<DdType>> model, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
+                for (auto const& formula : formulas) {
+                    std::cout << std::endl << "Model checking property: " << *formula << " ...";
+                    std::unique_ptr<storm::modelchecker::CheckResult> result;
+                    if (model->getType() == storm::models::ModelType::Dtmc) {
+                        std::shared_ptr<storm::models::symbolic::Dtmc<DdType>> dtmc = model->template as<storm::models::symbolic::Dtmc<DdType>>();
+                        storm::modelchecker::SymbolicDtmcPrctlModelChecker<DdType, double> modelchecker(*dtmc);
+                        if (modelchecker.canHandle(*formula)) {
+                            result = modelchecker.check(*formula);
+                        }
+                    } else if (model->getType() == storm::models::ModelType::Mdp) {
+                        std::shared_ptr<storm::models::symbolic::Mdp<DdType>> mdp = model->template as<storm::models::symbolic::Mdp<DdType>>();
+                        storm::modelchecker::SymbolicMdpPrctlModelChecker<DdType, double> modelchecker(*mdp);
+                        if (modelchecker.canHandle(*formula)) {
+                            result = modelchecker.check(*formula);
+                        }
+                    } else {
+                        STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This functionality is not yet implemented.");
+                    }
+                    
+                    if (result) {
+                        std::cout << " done." << std::endl;
+                        std::cout << "Result (initial states): ";
+                        result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<DdType>(model->getReachableStates(), model->getInitialStates()));
+                        std::cout << *result << std::endl;
+                    } else {
+                        std::cout << " skipped, because the modelling formalism is currently unsupported." << std::endl;
+                    }
+                }
+            }
+            
             template<typename ValueType>
             void buildAndCheckSymbolicModel(boost::optional<storm::prism::Program> const& program, std::vector<std::shared_ptr<storm::logic::Formula>> const& formulas) {
                 // Now we are ready to actually build the model.
@@ -346,10 +380,9 @@ namespace storm {
                         verifySparseModel<ValueType>(program, model->as<storm::models::sparse::Model<ValueType>>(), formulas);
                     } else if (model->isSymbolicModel()) {
                         if (storm::settings::generalSettings().getEngine() == storm::settings::modules::GeneralSettings::Engine::Hybrid) {
-                            verifySymbolicModel(program, model->as<storm::models::symbolic::Model<storm::dd::DdType::CUDD>>(), formulas);
+                            verifySymbolicModelWithHybridEngine(program, model->as<storm::models::symbolic::Model<storm::dd::DdType::CUDD>>(), formulas);
                         } else {
-                            // Not handled yet.
-                            STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This functionality is not yet implemented.");
+                            verifySymbolicModelWithSymbolicEngine(program, model->as<storm::models::symbolic::Model<storm::dd::DdType::CUDD>>(), formulas);
                         }
                     } else {
                         STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Invalid input model type.");
diff --git a/src/utility/graph.h b/src/utility/graph.h
index 45d90dffb..1f424c98d 100644
--- a/src/utility/graph.h
+++ b/src/utility/graph.h
@@ -723,12 +723,24 @@ namespace storm {
                         for(typename storm::storage::SparseMatrix<T>::const_iterator predecessorEntryIt = backwardTransitions.begin(currentState), predecessorEntryIte = backwardTransitions.end(currentState); predecessorEntryIt != predecessorEntryIte; ++predecessorEntryIt) {
                             if (phiStates.get(predecessorEntryIt->getColumn()) && !nextStates.get(predecessorEntryIt->getColumn())) {
                                 // Check whether the predecessor has only successors in the current state set for all of the
-                                // nondeterminstic choices.
-                                bool allSuccessorsInCurrentStatesForAllChoices = true;
-                                for (typename storm::storage::SparseMatrix<T>::const_iterator successorEntryIt = transitionMatrix.begin(nondeterministicChoiceIndices[predecessorEntryIt->getColumn()]), successorEntryIte = transitionMatrix.begin(nondeterministicChoiceIndices[predecessorEntryIt->getColumn() + 1]); successorEntryIt != successorEntryIte; ++successorEntryIt) {
-                                    if (!currentStates.get(successorEntryIt->getColumn())) {
-                                        allSuccessorsInCurrentStatesForAllChoices = false;
-                                        goto afterCheckLoop;
+                                // nondeterminstic choices and that for each choice there exists a successor that is already
+                                // in the next states.
+                                bool addToStatesWithProbability1 = true;
+                                for (uint_fast64_t row = nondeterministicChoiceIndices[predecessorEntryIt->getColumn()]; row < nondeterministicChoiceIndices[predecessorEntryIt->getColumn() + 1]; ++row) {
+                                    bool hasAtLeastOneSuccessorWithProbability1 = false;
+                                    for (typename storm::storage::SparseMatrix<T>::const_iterator successorEntryIt = transitionMatrix.begin(row), successorEntryIte = transitionMatrix.end(row); successorEntryIt != successorEntryIte; ++successorEntryIt) {
+                                        if (!currentStates.get(successorEntryIt->getColumn())) {
+                                            addToStatesWithProbability1 = false;
+                                            goto afterCheckLoop;
+                                        }
+                                        if (nextStates.get(successorEntryIt->getColumn())) {
+                                            hasAtLeastOneSuccessorWithProbability1 = true;
+                                        }
+                                    }
+                                    
+                                    if (!hasAtLeastOneSuccessorWithProbability1) {
+                                        addToStatesWithProbability1 = false;
+                                        break;
                                     }
                                 }
                                 
@@ -736,7 +748,7 @@ namespace storm {
                                 // If all successors for all nondeterministic choices are in the current state set, we
                                 // add it to the set of states for the next iteration and perform a backward search from
                                 // that state.
-                                if (allSuccessorsInCurrentStatesForAllChoices) {
+                                if (addToStatesWithProbability1) {
                                     nextStates.set(predecessorEntryIt->getColumn(), true);
                                     stack.push_back(predecessorEntryIt->getColumn());
                                 }
diff --git a/src/utility/solver.cpp b/src/utility/solver.cpp
index 14cd65396..0e818456c 100644
--- a/src/utility/solver.cpp
+++ b/src/utility/solver.cpp
@@ -25,6 +25,11 @@ namespace storm {
                 return std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>>(new storm::solver::SymbolicLinearEquationSolver<Type, ValueType>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs));
             }
             
+            template<storm::dd::DdType Type, typename ValueType>
+            std::unique_ptr<storm::solver::SymbolicMinMaxLinearEquationSolver<Type, ValueType>> SymbolicMinMaxLinearEquationSolverFactory<Type, ValueType>::create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, storm::dd::Bdd<Type> const& illegalMask, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& choiceVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const {
+                return std::unique_ptr<storm::solver::SymbolicMinMaxLinearEquationSolver<Type, ValueType>>(new storm::solver::SymbolicMinMaxLinearEquationSolver<Type, ValueType>(A, allRows, illegalMask, rowMetaVariables, columnMetaVariables, choiceVariables, rowColumnMetaVariablePairs));
+            }
+            
             template<storm::dd::DdType Type>
             std::unique_ptr<storm::solver::SymbolicGameSolver<Type>> SymbolicGameSolverFactory<Type>::create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables) const {
                 return std::unique_ptr<storm::solver::SymbolicGameSolver<Type>>(new storm::solver::SymbolicGameSolver<Type>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, player1Variables, player2Variables));
@@ -117,6 +122,7 @@ namespace storm {
             }
             
             template class SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>;
+            template class SymbolicMinMaxLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>;
             template class SymbolicGameSolverFactory<storm::dd::DdType::CUDD>;
             template class LinearEquationSolverFactory<double>;
             template class GmmxxLinearEquationSolverFactory<double>;
diff --git a/src/utility/solver.h b/src/utility/solver.h
index cd56695f0..715ad3972 100644
--- a/src/utility/solver.h
+++ b/src/utility/solver.h
@@ -2,7 +2,7 @@
 #define STORM_UTILITY_SOLVER_H_
 
 #include "src/solver/SymbolicGameSolver.h"
-
+#include "src/solver/SymbolicMinMaxLinearEquationSolver.h"
 #include "src/solver/SymbolicLinearEquationSolver.h"
 #include "src/solver/LinearEquationSolver.h"
 #include "src/solver/NativeLinearEquationSolver.h"
@@ -23,6 +23,12 @@ namespace storm {
                 virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>> create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const;
             };
             
+            template<storm::dd::DdType Type, typename ValueType>
+            class SymbolicMinMaxLinearEquationSolverFactory {
+                public:
+                virtual std::unique_ptr<storm::solver::SymbolicMinMaxLinearEquationSolver<Type, ValueType>> create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, storm::dd::Bdd<Type> const& illegalMask, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& choiceVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const;
+            };
+            
             template<storm::dd::DdType Type>
             class SymbolicGameSolverFactory {
             public:
diff --git a/test/functional/modelchecker/GmmxxHybridMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxHybridMdpPrctlModelCheckerTest.cpp
index 8162d0a5c..31ab04a6b 100644
--- a/test/functional/modelchecker/GmmxxHybridMdpPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/GmmxxHybridMdpPrctlModelCheckerTest.cpp
@@ -102,8 +102,8 @@ TEST(GmmxxHybridMdpPrctlModelCheckerTest, Dice) {
     result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
     storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult7 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
     
-    EXPECT_NEAR(7.3333283960819244, quantitativeResult7.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
-    EXPECT_NEAR(7.3333283960819244, quantitativeResult7.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333294987678528, quantitativeResult7.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333294987678528, quantitativeResult7.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
     
     formula = parser.parseFromString("Rmax=? [F \"done\"]");
     
@@ -111,8 +111,8 @@ TEST(GmmxxHybridMdpPrctlModelCheckerTest, Dice) {
     result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
     storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult8 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
     
-    EXPECT_NEAR(7.3333283960819244, quantitativeResult8.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
-    EXPECT_NEAR(7.3333283960819244, quantitativeResult8.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333294987678528, quantitativeResult8.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333294987678528, quantitativeResult8.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
 }
 
 TEST(GmmxxHybridMdpPrctlModelCheckerTest, AsynchronousLeader) {
@@ -181,8 +181,8 @@ TEST(GmmxxHybridMdpPrctlModelCheckerTest, AsynchronousLeader) {
     result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
     storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult5 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
     
-    EXPECT_NEAR(4.2856925589077264, quantitativeResult5.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
-    EXPECT_NEAR(4.2856925589077264, quantitativeResult5.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856896106114934, quantitativeResult5.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856896106114934, quantitativeResult5.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
     
     formula = parser.parseFromString("Rmax=? [F \"elected\"]");
     
@@ -190,6 +190,6 @@ TEST(GmmxxHybridMdpPrctlModelCheckerTest, AsynchronousLeader) {
     result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
     storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult6 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
     
-    EXPECT_NEAR(4.2856953906798676, quantitativeResult6.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
-    EXPECT_NEAR(4.2856953906798676, quantitativeResult6.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856896106114934, quantitativeResult6.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856896106114934, quantitativeResult6.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
 }
diff --git a/test/functional/modelchecker/NativeHybridMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/NativeHybridMdpPrctlModelCheckerTest.cpp
index 5b67f31e2..e7bcd7e65 100644
--- a/test/functional/modelchecker/NativeHybridMdpPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/NativeHybridMdpPrctlModelCheckerTest.cpp
@@ -100,8 +100,8 @@ TEST(NativeHybridMdpPrctlModelCheckerTest, Dice) {
     result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
     storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult7 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
     
-    EXPECT_NEAR(7.3333283960819244, quantitativeResult7.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
-    EXPECT_NEAR(7.3333283960819244, quantitativeResult7.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333294987678528, quantitativeResult7.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333294987678528, quantitativeResult7.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
     
     formula = parser.parseFromString("Rmax=? [F \"done\"]");
     
@@ -109,8 +109,8 @@ TEST(NativeHybridMdpPrctlModelCheckerTest, Dice) {
     result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
     storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult8 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
     
-    EXPECT_NEAR(7.3333283960819244, quantitativeResult8.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
-    EXPECT_NEAR(7.3333283960819244, quantitativeResult8.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333294987678528, quantitativeResult8.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333294987678528, quantitativeResult8.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
 }
 
 TEST(NativeHybridMdpPrctlModelCheckerTest, AsynchronousLeader) {
@@ -164,7 +164,7 @@ TEST(NativeHybridMdpPrctlModelCheckerTest, AsynchronousLeader) {
     EXPECT_NEAR(0.0625, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
     EXPECT_NEAR(0.0625, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
     
-    formula = parser.parseFromString("Pmax=? [F \"elected\"]");
+    formula = parser.parseFromString("Pmax=? [F<=25 \"elected\"]");
     
     result = checker.check(*formula);
     result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
@@ -179,8 +179,8 @@ TEST(NativeHybridMdpPrctlModelCheckerTest, AsynchronousLeader) {
     result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
     storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult5 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
     
-    EXPECT_NEAR(4.2856925589077264, quantitativeResult5.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
-    EXPECT_NEAR(4.2856925589077264, quantitativeResult5.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856896106114934, quantitativeResult5.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856896106114934, quantitativeResult5.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
     
     formula = parser.parseFromString("Rmax=? [F \"elected\"]");
     
@@ -188,6 +188,6 @@ TEST(NativeHybridMdpPrctlModelCheckerTest, AsynchronousLeader) {
     result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
     storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult6 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
     
-    EXPECT_NEAR(4.2856953906798676, quantitativeResult6.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
-    EXPECT_NEAR(4.2856953906798676, quantitativeResult6.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856896106114934, quantitativeResult6.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856896106114934, quantitativeResult6.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
 }
diff --git a/test/functional/modelchecker/SymbolicMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SymbolicMdpPrctlModelCheckerTest.cpp
new file mode 100644
index 000000000..0642e19d2
--- /dev/null
+++ b/test/functional/modelchecker/SymbolicMdpPrctlModelCheckerTest.cpp
@@ -0,0 +1,192 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/logic/Formulas.h"
+#include "src/utility/solver.h"
+#include "src/modelchecker/prctl/SymbolicMdpPrctlModelChecker.h"
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+#include "src/parser/FormulaParser.h"
+#include "src/parser/PrismParser.h"
+#include "src/builder/DdPrismModelBuilder.h"
+#include "src/models/symbolic/Dtmc.h"
+#include "src/models/symbolic/StandardRewardModel.h"
+#include "src/settings/SettingsManager.h"
+
+#include "src/settings/modules/GeneralSettings.h"
+
+TEST(SymbolicMdpPrctlModelCheckerTest, Dice) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/two_dice.nm");
+    
+    // A parser that we use for conveniently constructing the formulas.
+    storm::parser::FormulaParser parser;
+    
+    // Build the die model with its reward model.
+#ifdef WINDOWS
+    storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#else
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#endif
+    options.buildAllRewardModels = false;
+    options.rewardModelsToBuild.insert("coinflips");
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    EXPECT_EQ(169ul, model->getNumberOfStates());
+    EXPECT_EQ(436ul, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Mdp);
+    
+    std::shared_ptr<storm::models::symbolic::Mdp<storm::dd::DdType::CUDD>> mdp = model->as<storm::models::symbolic::Mdp<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::SymbolicMdpPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*mdp, std::unique_ptr<storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>>(new storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>()));
+    
+    std::shared_ptr<storm::logic::Formula> formula = parser.parseFromString("Pmin=? [F \"two\"]");
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0277777612209320068, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0277777612209320068, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    formula = parser.parseFromString("Pmax=? [F \"two\"]");
+    
+    result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0277777612209320068, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0277777612209320068, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    formula = parser.parseFromString("Pmin=? [F \"three\"]");
+    
+    result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0555555224418640136, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0555555224418640136, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    formula = parser.parseFromString("Pmax=? [F \"three\"]");
+    
+    result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult4 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0555555224418640136, quantitativeResult4.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0555555224418640136, quantitativeResult4.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    formula = parser.parseFromString("Pmin=? [F \"four\"]");
+    
+    result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult5 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.083333283662796020508, quantitativeResult5.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.083333283662796020508, quantitativeResult5.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    formula = parser.parseFromString("Pmax=? [F \"four\"]");
+    
+    result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult6 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.083333283662796020508, quantitativeResult6.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.083333283662796020508, quantitativeResult6.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    formula = parser.parseFromString("Rmin=? [F \"done\"]");
+    
+    result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult7 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(7.3333272933959961, quantitativeResult7.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333272933959961, quantitativeResult7.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    formula = parser.parseFromString("Rmax=? [F \"done\"]");
+    
+    result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult8 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(7.3333272933959961, quantitativeResult8.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333272933959961, quantitativeResult8.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}
+
+TEST(SymbolicMdpPrctlModelCheckerTest, AsynchronousLeader) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/leader4.nm");
+    
+    // A parser that we use for conveniently constructing the formulas.
+    storm::parser::FormulaParser parser;
+    
+    // Build the die model with its reward model.
+#ifdef WINDOWS
+    storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#else
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#endif
+    options.buildAllRewardModels = false;
+    options.rewardModelsToBuild.insert("rounds");
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    EXPECT_EQ(3172ul, model->getNumberOfStates());
+    EXPECT_EQ(7144ul, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Mdp);
+    
+    std::shared_ptr<storm::models::symbolic::Mdp<storm::dd::DdType::CUDD>> mdp = model->as<storm::models::symbolic::Mdp<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::SymbolicMdpPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*mdp, std::unique_ptr<storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>>(new storm::utility::solver::SymbolicMinMaxLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>()));
+    
+    std::shared_ptr<storm::logic::Formula> formula = parser.parseFromString("Pmin=? [F \"elected\"]");
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    formula = parser.parseFromString("Pmax=? [F \"elected\"]");
+    
+    result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    formula = parser.parseFromString("Pmin=? [F<=25 \"elected\"]");
+    
+    result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0625, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0625, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    formula = parser.parseFromString("Pmax=? [F<=25 \"elected\"]");
+    
+    result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult4 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0625, quantitativeResult4.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0625, quantitativeResult4.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    formula = parser.parseFromString("Rmin=? [F \"elected\"]");
+    
+    result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult5 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(4.2856890848060498, quantitativeResult5.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856890848060498, quantitativeResult5.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    formula = parser.parseFromString("Rmax=? [F \"elected\"]");
+    
+    result = checker.check(*formula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult6 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(4.2856890848060498, quantitativeResult6.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856890848060498, quantitativeResult6.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}