diff --git a/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.cpp
new file mode 100644
index 000000000..022eca0af
--- /dev/null
+++ b/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.cpp
@@ -0,0 +1,263 @@
+#include "src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h"
+
+#include "src/storage/dd/CuddOdd.h"
+
+#include "src/utility/macros.h"
+#include "src/utility/graph.h"
+
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+#include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
+
+#include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/InvalidPropertyException.h"
+
+namespace storm {
+    namespace modelchecker {
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicDtmcPrctlModelChecker<DdType, ValueType>::SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory) : SymbolicPropositionalModelChecker<DdType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicDtmcPrctlModelChecker<DdType, ValueType>::SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model) : SymbolicPropositionalModelChecker<DdType>(model), linearEquationSolverFactory(new storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>()) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        bool SymbolicDtmcPrctlModelChecker<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> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+            // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
+            // probability 0 and 1 of satisfying the until-formula.
+            std::pair<storm::dd::Bdd<DdType>, storm::dd::Bdd<DdType>> statesWithProbability01 = storm::utility::graph::performProb01(model, transitionMatrix, phiStates, psiStates);
+            storm::dd::Bdd<DdType> maybeStates = !statesWithProbability01.first && !statesWithProbability01.second && model.getReachableStates();
+            
+            // Perform some logging.
+            STORM_LOG_INFO("Found " << statesWithProbability01.first.getNonZeroCount() << " 'no' states.");
+            STORM_LOG_INFO("Found " << statesWithProbability01.second.getNonZeroCount() << " 'yes' states.");
+            STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
+            
+            // Check whether we need to compute exact probabilities for some states.
+            if (qualitative) {
+                // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1.
+                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd() + maybeStates.toAdd() * model.getManager().getConstant(0.5)));
+            } else {
+                // If there are maybe states, we need to solve an equation system.
+                if (!maybeStates.isZero()) {
+                    // Create the ODD for the translation between symbolic and explicit storage.
+                    storm::dd::Odd<DdType> odd(maybeStates);
+                    
+                    // Create the matrix and the vector for the equation system.
+                    storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                    
+                    // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                    // non-maybe states in the matrix.
+                    storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                    
+                    // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
+                    // maybe states.
+                    storm::dd::Add<DdType> prob1StatesAsColumn = statesWithProbability01.second.toAdd();
+                    prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs());
+                    storm::dd::Add<DdType> subvector = submatrix * prob1StatesAsColumn;
+                    subvector = subvector.sumAbstract(model.getColumnVariables());
+                    
+                    // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
+                    // for solving the equation system (i.e. compute (I-A)).
+                    submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                    submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
+                    
+                    // Solve the equation system.
+                    std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                    storm::dd::Add<DdType> result = solver->solveEquationSystem(model.getManager().getConstant(0.5) * maybeStatesAdd, subvector);
+                    
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd() + result));
+                } else {
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd()));
+                }
+            }
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
+            std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return this->computeUntilProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->computeNextProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector())));
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeNextProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates) {
+            storm::dd::Add<DdType> result = transitionMatrix * nextStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd();
+            return result.sumAbstract(model.getColumnVariables());
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(pathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
+            std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
+            std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return this->computeBoundedUntilProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+            // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
+            // probability 0 or 1 of satisfying the until-formula.
+            storm::dd::Bdd<DdType> statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(model, transitionMatrix.notZero(), phiStates, psiStates, stepBound);
+            storm::dd::Bdd<DdType> maybeStates = statesWithProbabilityGreater0 && !psiStates && model.getReachableStates();
+            
+            // If there are maybe states, we need to perform matrix-vector multiplications.
+            if (!maybeStates.isZero()) {
+                // Create the ODD for the translation between symbolic and explicit storage.
+                storm::dd::Odd<DdType> odd(maybeStates);
+                
+                // Create the matrix and the vector for the equation system.
+                storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                
+                // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                // non-maybe states in the matrix.
+                storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                
+                // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
+                // maybe states.
+                storm::dd::Add<DdType> prob1StatesAsColumn = psiStates.toAdd().swapVariables(model.getRowColumnMetaVariablePairs());
+                storm::dd::Add<DdType> subvector = (submatrix * prob1StatesAsColumn).sumAbstract(model.getColumnVariables());
+                
+                // Finally cut away all columns targeting non-maybe states.
+                submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                
+                // Perform the matrix-vector multiplication.
+                std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getManager().getAddZero(), &subvector, stepBound);
+                
+                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd() + result));
+            } else {
+                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd()));
+            }
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
+            return this->computeCumulativeRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeCumulativeRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+            // Only compute the result if the model has at least one reward this->getModel().
+            STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+            
+            // Compute the reward vector to add in each step based on the available reward models.
+            storm::dd::Add<DdType> totalRewardVector = model.hasStateRewards() ? model.getStateRewardVector() : model.getManager().getAddZero();
+            if (model.hasTransitionRewards()) {
+                totalRewardVector += (transitionMatrix * model.getTransitionRewardMatrix()).sumAbstract(model.getColumnVariables());
+            }
+            
+            // Perform the matrix-vector multiplication.
+            std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix, model.getReachableStates(), model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+            storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getManager().getAddZero(), &totalRewardVector, stepBound);
+            
+            return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), result));
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
+            return this->computeInstantaneousRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory) {
+            // Only compute the result if the model has at least one reward this->getModel().
+            STORM_LOG_THROW(model.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+            
+            // Perform the matrix-vector multiplication.
+            std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix, model.getReachableStates(), model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+            storm::dd::Add<DdType> result = solver->performMatrixVectorMultiplication(model.getStateRewardVector(), nullptr, stepBound);
+            
+            return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), result));
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return this->computeReachabilityRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> SymbolicDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory, bool qualitative) {
+            
+            // Only compute the result if there is at least one reward model.
+            STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+            
+            // Determine which states have a reward of infinity by definition.
+            storm::dd::Bdd<DdType> infinityStates = storm::utility::graph::performProb1(model, transitionMatrix.notZero(), model.getReachableStates(), targetStates);
+            infinityStates = !infinityStates && model.getReachableStates();
+            storm::dd::Bdd<DdType> maybeStates = (!targetStates && !infinityStates) && model.getReachableStates();
+            STORM_LOG_INFO("Found " << infinityStates.getNonZeroCount() << " 'infinity' states.");
+            STORM_LOG_INFO("Found " << targetStates.getNonZeroCount() << " 'target' states.");
+            STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
+            
+            // Check whether we need to compute exact rewards for some states.
+            if (qualitative) {
+                // Set the values for all maybe-states to 1 to indicate that their reward values
+                // are neither 0 nor infinity.
+                return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + maybeStates.toAdd() * model.getManager().getConstant(storm::utility::one<ValueType>())));
+            } else {
+                // If there are maybe states, we need to solve an equation system.
+                if (!maybeStates.isZero()) {
+                    // Create the ODD for the translation between symbolic and explicit storage.
+                    storm::dd::Odd<DdType> odd(maybeStates);
+                    
+                    // Create the matrix and the vector for the equation system.
+                    storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                    
+                    // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                    // non-maybe states in the matrix.
+                    storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                    
+                    // Then compute the state reward vector to use in the computation.
+                    storm::dd::Add<DdType> subvector = stateRewardVector ? maybeStatesAdd * stateRewardVector.get() : model.getManager().getAddZero();
+                    if (transitionRewardMatrix) {
+                        subvector += (submatrix * transitionRewardMatrix.get()).sumAbstract(model.getColumnVariables());
+                    }
+                    
+                    // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
+                    // for solving the equation system (i.e. compute (I-A)).
+                    submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                    submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
+                    
+                    // Solve the equation system.
+                    std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                    storm::dd::Add<DdType> result = solver->solveEquationSystem(model.getManager().getConstant(0.5) * maybeStatesAdd, subvector);
+                    
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + result));
+                } else {
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>())));
+                }
+            }
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::models::symbolic::Dtmc<DdType> const& SymbolicDtmcPrctlModelChecker<DdType, ValueType>::getModel() const {
+            return this->template getModelAs<storm::models::symbolic::Dtmc<DdType>>();
+        }
+        
+        template class SymbolicDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double>;
+    }
+}
\ No newline at end of file
diff --git a/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h
new file mode 100644
index 000000000..9dca0335b
--- /dev/null
+++ b/src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h
@@ -0,0 +1,49 @@
+#ifndef STORM_MODELCHECKER_SYMBOLICDTMCPRCTLMODELCHECKER_H_
+#define STORM_MODELCHECKER_SYMBOLICDTMCPRCTLMODELCHECKER_H_
+
+#include "src/modelchecker/propositional/SymbolicPropositionalModelChecker.h"
+#include "src/models/symbolic/Dtmc.h"
+#include "src/utility/solver.h"
+
+namespace storm {
+    namespace modelchecker {
+        template<storm::dd::DdType DdType, typename ValueType>
+        class SymbolicCtmcCslModelChecker;
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        class SymbolicDtmcPrctlModelChecker : public SymbolicPropositionalModelChecker<DdType> {
+        public:
+            friend class SymbolicCtmcCslModelChecker<DdType, ValueType>;
+            
+            explicit SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model);
+            explicit SymbolicDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory);
+            
+            // 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, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            
+        protected:
+            storm::models::symbolic::Dtmc<DdType> const& getModel() const override;
+            
+        private:
+            // The methods that perform the actual checking.
+            static std::unique_ptr<CheckResult> computeBoundedUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+            static storm::dd::Add<DdType> computeNextProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates);
+            static std::unique_ptr<CheckResult> computeUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+            static std::unique_ptr<CheckResult> computeCumulativeRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+            static std::unique_ptr<CheckResult> computeInstantaneousRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory);
+            static std::unique_ptr<CheckResult> computeReachabilityRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType> const& linearEquationSolverFactory, bool qualitative);
+            
+            // An object that is used for retrieving linear equation solvers.
+            std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<DdType, ValueType>> linearEquationSolverFactory;
+        };
+        
+    } // namespace modelchecker
+} // namespace storm
+
+#endif /* STORM_MODELCHECKER_SYMBOLICDTMCPRCTLMODELCHECKER_H_ */
diff --git a/src/solver/FullySymbolicGameSolver.cpp b/src/solver/FullySymbolicGameSolver.cpp
deleted file mode 100644
index f02b50ae8..000000000
--- a/src/solver/FullySymbolicGameSolver.cpp
+++ /dev/null
@@ -1,70 +0,0 @@
-#include "src/solver/FullySymbolicGameSolver.h"
-
-#include "src/storage/dd/CuddBdd.h"
-#include "src/storage/dd/CuddAdd.h"
-
-namespace storm {
-    namespace solver {
-
-        template<storm::dd::DdType Type>
-        FullySymbolicGameSolver<Type>::FullySymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, 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, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : SymbolicGameSolver<Type>(gameMatrix, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, player1Variables, player2Variables), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), relative(relative) {
-            // Intentionally left empty.
-        }
-        
-        template<storm::dd::DdType Type>
-        FullySymbolicGameSolver<Type>::FullySymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, 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) : SymbolicGameSolver<Type>(gameMatrix, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, player1Variables, player2Variables) {
-            // Get the settings object to customize solving.
-            storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
-            storm::settings::modules::GeneralSettings const& generalSettings = storm::settings::generalSettings();
-            
-            // Get appropriate settings.
-            maximalNumberOfIterations = settings.getMaximalIterationCount();
-            precision = settings.getPrecision();
-            relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative;
-        }
-        
-        template<storm::dd::DdType Type>
-        storm::dd::Add<Type> FullySymbolicGameSolver<Type>::solveGame(bool player1Min, bool player2Min, storm::dd::Add<Type> const& x, storm::dd::Add<Type> const& b) const {
-            // Set up the environment.
-            storm::dd::Add<Type> xCopy = x;
-            uint_fast64_t iterations = 0;
-            bool converged = false;
-
-            while (!converged && iterations < maximalNumberOfIterations) {
-                // Compute tmp = A * x + b
-                storm::dd::Add<Type> xCopyAsColumn = xCopy.swapVariables(this->rowColumnMetaVariablePairs);
-                storm::dd::Add<Type> tmp = this->gameMatrix.multiplyMatrix(xCopyAsColumn, this->columnMetaVariables);
-                tmp += b;
-                
-                // Now abstract from player 2 and player 1 variables.
-                // TODO: can this order be reversed?
-                if (player2Min) {
-                    tmp = tmp.minAbstract(this->player2Variables);
-                } else {
-                    tmp = tmp.maxAbstract(this->player2Variables);
-                }
-                
-                if (player1Min) {
-                    tmp = tmp.minAbstract(this->player1Variables);
-                } else {
-                    tmp = tmp.maxAbstract(this->player1Variables);
-                }
-
-                // 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 class FullySymbolicGameSolver<storm::dd::DdType::CUDD>;
-        
-    }
-}
\ No newline at end of file
diff --git a/src/solver/FullySymbolicGameSolver.h b/src/solver/FullySymbolicGameSolver.h
deleted file mode 100644
index f3f5f014d..000000000
--- a/src/solver/FullySymbolicGameSolver.h
+++ /dev/null
@@ -1,67 +0,0 @@
-#ifndef STORM_SOLVER_FULLYSYMBOLICGAMESOLVER_H_
-#define STORM_SOLVER_FULLYSYMBOLICGAMESOLVER_H_
-
-#include "src/solver/SymbolicGameSolver.h"
-
-namespace storm {
-    namespace solver {
-        
-        /*!
-         * A interface that represents an abstract symbolic game solver.
-         */
-        template<storm::dd::DdType Type>
-        class FullySymbolicGameSolver : public SymbolicGameSolver<Type> {
-        public:
-            /*!
-             * Constructs a symbolic game solver with the given meta variable sets and pairs.
-             *
-             * @param gameMatrix The matrix defining the coefficients of the game.
-             * @param allRows A BDD characterizing all rows of the equation system.
-             * @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 rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta
-             * variables.
-             * @param player1Variables The meta variables used to encode the player 1 choices.
-             * @param player2Variables The meta variables used to encode the player 2 choices.
-             * @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.
-             */
-            FullySymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, 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);
-            
-            /*!
-             * Constructs a symbolic game solver with the given meta variable sets and pairs.
-             *
-             * @param gameMatrix The matrix defining the coefficients of the game.
-             * @param allRows A BDD characterizing all rows of the equation system.
-             * @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 rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta
-             * variables.
-             * @param player1Variables The meta variables used to encode the player 1 choices.
-             * @param player2Variables The meta variables used to encode the player 2 choices.
-             * @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.
-             */
-            FullySymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, 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, double precision, uint_fast64_t maximalNumberOfIterations, bool relative);
-            
-            virtual storm::dd::Add<Type> solveGame(bool player1Min, bool player2Min, storm::dd::Add<Type> const& x, storm::dd::Add<Type> const& b) const override;
-            
-        private:
-            // 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_FULLYSYMBOLICGAMESOLVER_H_ */
diff --git a/src/solver/SymbolicGameSolver.cpp b/src/solver/SymbolicGameSolver.cpp
index 8224e24c9..cafa4eb04 100644
--- a/src/solver/SymbolicGameSolver.cpp
+++ b/src/solver/SymbolicGameSolver.cpp
@@ -5,12 +5,64 @@
 
 namespace storm {
     namespace solver {
-
+        
         template<storm::dd::DdType Type>
         SymbolicGameSolver<Type>::SymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, 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) : gameMatrix(gameMatrix), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), player1Variables(player1Variables), player2Variables(player2Variables) {
+            // Get the settings object to customize solving.
+            storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
+            storm::settings::modules::GeneralSettings const& generalSettings = storm::settings::generalSettings();
+            
+            // Get appropriate settings.
+            maximalNumberOfIterations = settings.getMaximalIterationCount();
+            precision = settings.getPrecision();
+            relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative;
+        }
+        
+        template<storm::dd::DdType Type>
+        SymbolicGameSolver<Type>::SymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, 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, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : gameMatrix(gameMatrix), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), player1Variables(player1Variables), player2Variables(player2Variables), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), relative(relative) {
             // Intentionally left empty.
         }
-     
+        
+        template<storm::dd::DdType Type>
+        storm::dd::Add<Type> SymbolicGameSolver<Type>::solveGame(bool player1Min, bool player2Min, storm::dd::Add<Type> const& x, storm::dd::Add<Type> const& b) const {
+            // Set up the environment.
+            storm::dd::Add<Type> xCopy = x;
+            uint_fast64_t iterations = 0;
+            bool converged = false;
+            
+            while (!converged && iterations < maximalNumberOfIterations) {
+                // Compute tmp = A * x + b
+                storm::dd::Add<Type> xCopyAsColumn = xCopy.swapVariables(this->rowColumnMetaVariablePairs);
+                storm::dd::Add<Type> tmp = this->gameMatrix.multiplyMatrix(xCopyAsColumn, this->columnMetaVariables);
+                tmp += b;
+                
+                // Now abstract from player 2 and player 1 variables.
+                if (player2Min) {
+                    tmp = tmp.minAbstract(this->player2Variables);
+                } else {
+                    tmp = tmp.maxAbstract(this->player2Variables);
+                }
+                
+                if (player1Min) {
+                    tmp = tmp.minAbstract(this->player1Variables);
+                } else {
+                    tmp = tmp.maxAbstract(this->player1Variables);
+                }
+                
+                // 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 class SymbolicGameSolver<storm::dd::DdType::CUDD>;
         
     }
diff --git a/src/solver/SymbolicGameSolver.h b/src/solver/SymbolicGameSolver.h
index bbf5d539e..ff5d64056 100644
--- a/src/solver/SymbolicGameSolver.h
+++ b/src/solver/SymbolicGameSolver.h
@@ -10,7 +10,7 @@ namespace storm {
     namespace solver {
         
         /*!
-         * A interface that represents an abstract symbolic game solver.
+         * An interface that represents an abstract symbolic game solver.
          */
         template<storm::dd::DdType Type>
         class SymbolicGameSolver {
@@ -29,6 +29,24 @@ namespace storm {
              */
             SymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, 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);
             
+            /*!
+             * Constructs a symbolic game solver with the given meta variable sets and pairs.
+             *
+             * @param gameMatrix The matrix defining the coefficients of the game.
+             * @param allRows A BDD characterizing all rows of the equation system.
+             * @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 rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta
+             * variables.
+             * @param player1Variables The meta variables used to encode the player 1 choices.
+             * @param player2Variables The meta variables used to encode the player 2 choices.
+             * @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.
+             */
+            SymbolicGameSolver(storm::dd::Add<Type> const& gameMatrix, 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, double precision, uint_fast64_t maximalNumberOfIterations, bool relative);
+            
             /*!
              * Solves the equation system x = min/max(A*x + b) given by the parameters. Note that the matrix A has
              * to be given upon construction time of the solver object.
@@ -39,7 +57,7 @@ namespace storm {
              * @param b The vector to add after matrix-vector multiplication.
              * @return The solution vector.
              */
-            virtual storm::dd::Add<Type> solveGame(bool player1Min, bool player2Min, storm::dd::Add<Type> const& x, storm::dd::Add<Type> const& b) const = 0;
+            virtual storm::dd::Add<Type> solveGame(bool player1Min, bool player2Min, storm::dd::Add<Type> const& x, storm::dd::Add<Type> const& b) const;
         
         protected:
             // The matrix defining the coefficients of the linear equation system.
@@ -62,6 +80,15 @@ namespace storm {
             
             // The player 2 variables.
             std::set<storm::expressions::Variable> player2Variables;
+            
+            // 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
diff --git a/src/solver/SymbolicLinearEquationSolver.cpp b/src/solver/SymbolicLinearEquationSolver.cpp
new file mode 100644
index 000000000..24a9d12b0
--- /dev/null
+++ b/src/solver/SymbolicLinearEquationSolver.cpp
@@ -0,0 +1,85 @@
+#include "src/solver/SymbolicLinearEquationSolver.h"
+
+#include "src/storage/dd/CuddDdManager.h"
+#include "src/storage/dd/CuddAdd.h"
+
+namespace storm {
+    namespace solver {
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> 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, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : A(A), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), relative(relative) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> 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) : A(A), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs) {
+            // Get the settings object to customize solving.
+            storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
+            storm::settings::modules::GeneralSettings const& generalSettings = storm::settings::generalSettings();
+            
+            // 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>  SymbolicLinearEquationSolver<DdType, ValueType>::solveEquationSystem(storm::dd::Add<DdType> const& x, storm::dd::Add<DdType> const& b) const {
+            // Start by computing the Jacobi decomposition of the matrix A.
+            storm::dd::Add<DdType> diagonal = x.getDdManager()->getAddOne();
+            for (auto const& pair : rowColumnMetaVariablePairs) {
+                diagonal *= x.getDdManager()->getIdentity(pair.first).equals(x.getDdManager()->getIdentity(pair.second));
+                diagonal *= x.getDdManager()->getRange(pair.first).toAdd() * x.getDdManager()->getRange(pair.second).toAdd();
+            }
+            diagonal *= allRows.toAdd();
+            
+            storm::dd::Add<DdType> lu = diagonal.ite(this->A.getDdManager()->getAddZero(), this->A);
+            storm::dd::Add<DdType> dinv = diagonal / (diagonal * this->A);
+            
+            // Set up additional environment variables.
+            storm::dd::Add<DdType> xCopy = x;
+            uint_fast64_t iterationCount = 0;
+            bool converged = false;
+            
+            while (!converged && iterationCount < maximalNumberOfIterations) {
+                storm::dd::Add<DdType> xCopyAsColumn = xCopy.swapVariables(this->rowColumnMetaVariablePairs);
+                storm::dd::Add<DdType> tmp = lu.multiplyMatrix(xCopyAsColumn, this->columnMetaVariables);
+                tmp = b - tmp;
+                tmp = tmp.swapVariables(this->rowColumnMetaVariablePairs);
+                tmp = dinv.multiplyMatrix(tmp, this->columnMetaVariables);
+                
+                // 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;
+                }
+                
+                // Increase iteration count so we can abort if convergence is too slow.
+                ++iterationCount;
+            }
+            
+            return xCopy;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType> SymbolicLinearEquationSolver<DdType, ValueType>::performMatrixVectorMultiplication(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;
+                }
+            }
+            
+            return xCopy;
+        }
+        
+        template class SymbolicLinearEquationSolver<storm::dd::DdType::CUDD, double>;
+        
+    }
+}
diff --git a/src/solver/SymbolicLinearEquationSolver.h b/src/solver/SymbolicLinearEquationSolver.h
new file mode 100644
index 000000000..6cd437673
--- /dev/null
+++ b/src/solver/SymbolicLinearEquationSolver.h
@@ -0,0 +1,104 @@
+#ifndef STORM_SOLVER_SYMBOLICLINEAREQUATIONSOLVER_H_
+#define STORM_SOLVER_SYMBOLICLINEAREQUATIONSOLVER_H_
+
+#include <memory>
+#include <set>
+#include <boost/variant.hpp>
+
+#include "src/storage/expressions/Variable.h"
+#include "src/storage/dd/Bdd.h"
+#include "src/storage/dd/Add.h"
+#include "src/storage/dd/Odd.h"
+
+namespace storm {
+    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 SymbolicLinearEquationSolver {
+        public:
+            /*!
+             * 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 diagonal An ADD characterizing the elements on the diagonal of the matrix.
+             * @param allRows A BDD characterizing all rows of the equation system.
+             * @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 rowColumnMetaVariablePairs The pairs of row meta variables and the corresponding column meta
+             * variables.
+             */
+            SymbolicLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> 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);
+            
+            /*!
+             * 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 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 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.
+             */
+            SymbolicLinearEquationSolver(storm::dd::Add<DdType> const& A, storm::dd::Bdd<DdType> 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, 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 x The solution vector that has to be computed. 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.
+             */
+            virtual storm::dd::Add<DdType> solveEquationSystem(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 x The initial vector with which to perform matrix-vector multiplication. Its length must be equal
+             * to the number of rows 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 rows of A.
+             * @return The solution of the equation system.
+             */
+            virtual storm::dd::Add<DdType> performMatrixVectorMultiplication(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;
+            
+            // The row variables.
+            std::set<storm::expressions::Variable> rowMetaVariables;
+            
+            // The column variables.
+            std::set<storm::expressions::Variable> columnMetaVariables;
+            
+            // 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_SYMBOLICLINEAREQUATIONSOLVER_H_ */
diff --git a/src/utility/solver.cpp b/src/utility/solver.cpp
index 9b1983151..2be080beb 100644
--- a/src/utility/solver.cpp
+++ b/src/utility/solver.cpp
@@ -2,7 +2,7 @@
 
 #include "src/settings/SettingsManager.h"
 
-#include "src/solver/FullySymbolicGameSolver.h"
+#include "src/solver/SymbolicGameSolver.h"
 
 #include "src/solver/NativeLinearEquationSolver.h"
 #include "src/solver/GmmxxLinearEquationSolver.h"
@@ -17,9 +17,14 @@
 namespace storm {
     namespace utility {
         namespace solver {
+            template<storm::dd::DdType Type, typename ValueType>
+            std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>> SymbolicLinearEquationSolverFactory<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 {
+                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>
             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::FullySymbolicGameSolver<Type>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, player1Variables, player2Variables));
+                return std::unique_ptr<storm::solver::SymbolicGameSolver<Type>>(new storm::solver::SymbolicGameSolver<Type>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, player1Variables, player2Variables));
             }
             
             template<typename ValueType>
@@ -86,6 +91,7 @@ namespace storm {
                 return factory->create(name);
             }
             
+            template class SymbolicLinearEquationSolverFactory<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 9fe4389b6..38cc7fbc6 100644
--- a/src/utility/solver.h
+++ b/src/utility/solver.h
@@ -3,6 +3,7 @@
 
 #include "src/solver/SymbolicGameSolver.h"
 
+#include "src/solver/SymbolicLinearEquationSolver.h"
 #include "src/solver/LinearEquationSolver.h"
 #include "src/solver/MinMaxLinearEquationSolver.h"
 #include "src/solver/LpSolver.h"
@@ -14,6 +15,12 @@
 namespace storm {
     namespace utility {
         namespace solver {
+            template<storm::dd::DdType Type, typename ValueType>
+            class SymbolicLinearEquationSolverFactory {
+            public:
+                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>
             class SymbolicGameSolverFactory {
             public:
diff --git a/test/functional/modelchecker/SymbolicDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SymbolicDtmcPrctlModelCheckerTest.cpp
new file mode 100644
index 000000000..02d6ad5ba
--- /dev/null
+++ b/test/functional/modelchecker/SymbolicDtmcPrctlModelCheckerTest.cpp
@@ -0,0 +1,172 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/logic/Formulas.h"
+#include "src/utility/solver.h"
+#include "src/modelchecker/prctl/SymbolicDtmcPrctlModelChecker.h"
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+#include "src/parser/PrismParser.h"
+#include "src/builder/DdPrismModelBuilder.h"
+#include "src/models/symbolic/Dtmc.h"
+#include "src/settings/SettingsManager.h"
+
+TEST(SymbolicDtmcPrctlModelCheckerTest, Die) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/die.pm");
+    
+    // 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.buildRewards = true;
+    options.rewardModelName = "coin_flips";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    EXPECT_EQ(13, model->getNumberOfStates());
+    EXPECT_EQ(20, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Dtmc);
+    
+    std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::SymbolicDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>>(new storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("one");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*eventuallyFormula);
+    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.0/6.0, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0/6.0, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("two");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    result = checker.check(*eventuallyFormula);
+    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.0/6.0, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0/6.0, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("three");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    result = checker.check(*eventuallyFormula);
+    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(1.0/6.0, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0/6.0, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    auto done = std::make_shared<storm::logic::AtomicLabelFormula>("done");
+    auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(done);
+    
+    result = checker.check(*reachabilityRewardFormula);
+    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(3.6666622161865234, quantitativeResult4.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(3.6666622161865234, quantitativeResult4.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}
+
+TEST(SymbolicDtmcPrctlModelCheckerTest, Crowds) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/crowds-5-5.pm");
+    
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program);
+    EXPECT_EQ(8607, model->getNumberOfStates());
+    EXPECT_EQ(15113, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Dtmc);
+    
+    std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::SymbolicDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>>(new storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observe0Greater1");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*eventuallyFormula);
+    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.33288236360191303, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.33288236360191303, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observeIGreater1");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    result = checker.check(*eventuallyFormula);
+    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.15222081144084315, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.15222081144084315, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observeOnlyTrueSender");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    result = checker.check(*eventuallyFormula);
+    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.3215392962289586, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.3215392962289586, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}
+
+TEST(SymbolicDtmcPrctlModelCheckerTest, SynchronousLeader) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/leader-3-5.pm");
+    
+    // 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.buildRewards = true;
+    options.rewardModelName = "num_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(273, model->getNumberOfStates());
+    EXPECT_EQ(397, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Dtmc);
+    
+    std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::SymbolicDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>>(new storm::utility::solver::SymbolicLinearEquationSolverFactory<storm::dd::DdType::CUDD, double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*eventuallyFormula);
+    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.0, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto trueFormula = std::make_shared<storm::logic::BooleanLiteralFormula>(true);
+    auto boundedUntilFormula = std::make_shared<storm::logic::BoundedUntilFormula>(trueFormula, labelFormula, 20);
+    
+    result = checker.check(*boundedUntilFormula);
+    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.99999989760000074, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.99999989760000074, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
+    
+    result = checker.check(*reachabilityRewardFormula);
+    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(1.0416666666666643, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0416666666666643, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}
+