diff --git a/src/builder/ExplicitPrismModelBuilder.cpp b/src/builder/ExplicitPrismModelBuilder.cpp
index 8ba671a94..bb6b0ede7 100644
--- a/src/builder/ExplicitPrismModelBuilder.cpp
+++ b/src/builder/ExplicitPrismModelBuilder.cpp
@@ -131,6 +131,7 @@ namespace storm {
             // all expressions in the program so we can then evaluate them without having to store the values of the
             // constants in the state (i.e., valuation).
             preparedProgram = preparedProgram.substituteConstants();
+                
             storm::prism::RewardModel rewardModel = storm::prism::RewardModel();
             
             // Select the appropriate reward model.
@@ -565,7 +566,11 @@ namespace storm {
                             }
                             
                             for (auto const& stateProbabilityPair : choice) {
-                                globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second / totalNumberOfChoices;
+                                if (discreteTimeModel) {
+                                    globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second / totalNumberOfChoices;
+                                } else {
+                                    globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second;
+                                }
                                 
                                 // Now add all rewards that match this choice.
                                 for (auto const& transitionReward : transitionRewards) {
@@ -581,7 +586,11 @@ namespace storm {
                             }
                             
                             for (auto const& stateProbabilityPair : choice) {
-                                globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second / totalNumberOfChoices;
+                                if (discreteTimeModel) {
+                                    globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second / totalNumberOfChoices;
+                                } else {
+                                    globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second;
+                                }
                                 
                                 // Now add all rewards that match this choice.
                                 for (auto const& transitionReward : transitionRewards) {
diff --git a/src/counterexamples/MILPMinimalLabelSetGenerator.h b/src/counterexamples/MILPMinimalLabelSetGenerator.h
index af6d91dd8..7bec64e63 100644
--- a/src/counterexamples/MILPMinimalLabelSetGenerator.h
+++ b/src/counterexamples/MILPMinimalLabelSetGenerator.h
@@ -629,7 +629,7 @@ namespace storm {
                 uint_fast64_t numberOfConstraintsCreated = 0;
 
                 for (auto label : choiceInformation.knownLabels) {
-                    storm::expressions::Expression constraint = variableInformation.labelToVariableMap.at(label) == solver.getConstant(0);
+                    storm::expressions::Expression constraint = variableInformation.labelToVariableMap.at(label) == solver.getConstant(1);
                     solver.addConstraint("KnownLabels" + std::to_string(numberOfConstraintsCreated), constraint);
                     ++numberOfConstraintsCreated;
                 }
@@ -929,10 +929,9 @@ namespace storm {
                 double maximalReachabilityProbability = 0;
                 if (checkThresholdFeasible) {
                     storm::modelchecker::SparseMdpPrctlModelChecker<T> modelchecker(labeledMdp);
-                    std::unique_ptr<storm::modelchecker::CheckResult> result = modelchecker.check(pathFormula);
-                    storm::modelchecker::ExplicitQuantitativeCheckResult<double> const& quantitativeResult = result->asExplicitQuantitativeCheckResult<double>();
+                    std::vector<T> result = modelchecker.computeUntilProbabilitiesHelper(false, phiStates, psiStates, false);
                     for (auto state : labeledMdp.getInitialStates()) {
-                        maximalReachabilityProbability = std::max(maximalReachabilityProbability, quantitativeResult[state]);
+                        maximalReachabilityProbability = std::max(maximalReachabilityProbability, result[state]);
                     }
                     STORM_LOG_THROW((strictBound && maximalReachabilityProbability >= probabilityThreshold) || (!strictBound && maximalReachabilityProbability > probabilityThreshold), storm::exceptions::InvalidArgumentException, "Given probability threshold " << probabilityThreshold << " can not be " << (strictBound ? "achieved" : "exceeded") << " in model with maximal reachability probability of " << maximalReachabilityProbability << ".");
                     std::cout << std::endl << "Maximal reachability in model is " << maximalReachabilityProbability << "." << std::endl << std::endl;
@@ -953,6 +952,8 @@ namespace storm {
                 //  (4.2) Construct constraint system.
                 buildConstraintSystem(*solver, labeledMdp, psiStates, stateInformation, choiceInformation, variableInformation, probabilityThreshold, strictBound, includeSchedulerCuts);
                 
+                solver->writeModelToFile("model.lp");
+                
                 // (4.3) Optimize the model.
                 solver->optimize();
                 
diff --git a/src/counterexamples/SMTMinimalCommandSetGenerator.h b/src/counterexamples/SMTMinimalCommandSetGenerator.h
index f099674dc..cdd04e8fd 100644
--- a/src/counterexamples/SMTMinimalCommandSetGenerator.h
+++ b/src/counterexamples/SMTMinimalCommandSetGenerator.h
@@ -1482,7 +1482,7 @@ namespace storm {
              * @param commandSet The currently chosen set of commands.
              * @param variableInformation A structure with information about the variables of the solver.
              */
-            static void analyzeInsufficientProbabilitySolution(storm::solver::SmtSolver& solver, storm::models::Mdp<T> const& subMdp, storm::models::Mdp<T> const& originalMdp, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, boost::container::flat_set<uint_fast64_t> const& commandSet, VariableInformation& variableInformation, RelevancyInformation const& relevancyInformation) {
+            static void analyzeInsufficientProbabilitySolution(storm::solver::SmtSolver& solver, storm::models::sparse::Mdp<T> const& subMdp, storm::models::sparse::Mdp<T> const& originalMdp, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, boost::container::flat_set<uint_fast64_t> const& commandSet, VariableInformation& variableInformation, RelevancyInformation const& relevancyInformation) {
 
                 LOG4CPLUS_DEBUG(logger, "Analyzing solution with insufficient probability.");
 
@@ -1627,10 +1627,9 @@ namespace storm {
                 double maximalReachabilityProbability = 0;
                 if (checkThresholdFeasible) {
                     storm::modelchecker::SparseMdpPrctlModelChecker<T> modelchecker(labeledMdp);
-                    std::unique_ptr<storm::modelchecker::CheckResult> result = modelchecker.check(pathFormula);
-                    storm::modelchecker::ExplicitQuantitativeCheckResult<double> const& quantitativeResult = result->asExplicitQuantitativeCheckResult<double>();
+                    std::vector<T> result = modelchecker.computeUntilProbabilitiesHelper(false, phiStates, psiStates, false);
                     for (auto state : labeledMdp.getInitialStates()) {
-                        maximalReachabilityProbability = std::max(maximalReachabilityProbability, quantitativeResult[state]);
+                        maximalReachabilityProbability = std::max(maximalReachabilityProbability, result[state]);
                     }
                     STORM_LOG_THROW((strictBound && maximalReachabilityProbability >= probabilityThreshold) || (!strictBound && maximalReachabilityProbability > probabilityThreshold), storm::exceptions::InvalidArgumentException, "Given probability threshold " << probabilityThreshold << " can not be " << (strictBound ? "achieved" : "exceeded") << " in model with maximal reachability probability of " << maximalReachabilityProbability << ".");
                     std::cout << std::endl << "Maximal reachability in model is " << maximalReachabilityProbability << "." << std::endl << std::endl;
diff --git a/src/modelchecker/AbstractModelChecker.cpp b/src/modelchecker/AbstractModelChecker.cpp
index 7004245a2..c7ac1b18b 100644
--- a/src/modelchecker/AbstractModelChecker.cpp
+++ b/src/modelchecker/AbstractModelChecker.cpp
@@ -87,7 +87,7 @@ namespace storm {
         std::unique_ptr<CheckResult> AbstractModelChecker::computeLongRunAverage(storm::logic::StateFormula const& eventuallyFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This model checker does not support the computation of long-run averages.");
         }
-
+        
         std::unique_ptr<CheckResult> AbstractModelChecker::computeExpectedTimes(storm::logic::EventuallyFormula const& eventuallyFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This model checker does not support the computation of expected times.");
         }
@@ -164,6 +164,12 @@ namespace storm {
             std::unique_ptr<CheckResult> result;
             if (stateFormula.hasOptimalityType()) {
                 result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), qualitative, stateFormula.getOptimalityType());
+            } else if (stateFormula.hasBound()) {
+                if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || stateFormula.getComparisonType() == storm::logic::ComparisonType::LessEqual) {
+                    result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), qualitative, storm::logic::OptimalityType::Maximize);
+                } else {
+                    result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), qualitative, storm::logic::OptimalityType::Minimize);
+                }
             } else {
                 result = this->computeProbabilities(stateFormula.getSubformula().asPathFormula(), qualitative);
             }
@@ -190,6 +196,12 @@ namespace storm {
             std::unique_ptr<CheckResult> result;
             if (stateFormula.hasOptimalityType()) {
                 result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), qualitative, stateFormula.getOptimalityType());
+            } else if (stateFormula.hasBound()) {
+                if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || stateFormula.getComparisonType() == storm::logic::ComparisonType::LessEqual) {
+                    result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), qualitative, storm::logic::OptimalityType::Maximize);
+                } else {
+                    result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), qualitative, storm::logic::OptimalityType::Minimize);
+                }
             } else {
                 result = this->computeRewards(stateFormula.getSubformula().asRewardPathFormula(), qualitative);
             }
@@ -216,6 +228,12 @@ namespace storm {
             std::unique_ptr<CheckResult> result;
             if (stateFormula.hasOptimalityType()) {
                 result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), qualitative, stateFormula.getOptimalityType());
+            } else if (stateFormula.hasBound()) {
+                if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || stateFormula.getComparisonType() == storm::logic::ComparisonType::LessEqual) {
+                    result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), qualitative, storm::logic::OptimalityType::Maximize);
+                } else {
+                    result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), qualitative, storm::logic::OptimalityType::Minimize);
+                }
             } else {
                 result = this->computeExpectedTimes(stateFormula.getSubformula().asEventuallyFormula(), qualitative);
             }
@@ -231,13 +249,19 @@ namespace storm {
         std::unique_ptr<CheckResult> AbstractModelChecker::checkLongRunAverageOperatorFormula(storm::logic::LongRunAverageOperatorFormula const& stateFormula) {
             STORM_LOG_THROW(stateFormula.getSubformula().isEventuallyFormula(), storm::exceptions::InvalidArgumentException, "The given formula is invalid.");
             
+            std::unique_ptr<CheckResult> result;
             if (stateFormula.hasOptimalityType()) {
-                return this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false, stateFormula.getOptimalityType());
+                result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false, stateFormula.getOptimalityType());
+            } else if (stateFormula.hasBound()) {
+                if (stateFormula.getComparisonType() == storm::logic::ComparisonType::Less || stateFormula.getComparisonType() == storm::logic::ComparisonType::LessEqual) {
+                    result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false, storm::logic::OptimalityType::Maximize);
+                } else {
+                    result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false, storm::logic::OptimalityType::Minimize);
+                }
             } else {
-                return this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false);
+                result = this->computeLongRunAverage(stateFormula.getSubformula().asStateFormula(), false);
             }
             
-            std::unique_ptr<CheckResult> result;
             if (stateFormula.hasBound()) {
                 return result->asQuantitativeCheckResult().compareAgainstBound(stateFormula.getComparisonType(), stateFormula.getBound());
             } else {
diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
new file mode 100644
index 000000000..480afd57d
--- /dev/null
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
@@ -0,0 +1,160 @@
+#include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
+
+#include <vector>
+
+#include "src/utility/macros.h"
+#include "src/utility/vector.h"
+#include "src/utility/graph.h"
+#include "src/utility/solver.h"
+
+#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
+#include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
+#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
+
+#include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/InvalidPropertyException.h"
+#include "src/exceptions/NotImplementedException.h"
+
+namespace storm {
+    namespace modelchecker {
+        template<class ValueType>
+        SparseCtmcCslModelChecker<ValueType>::SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model) : SparsePropositionalModelChecker<ValueType>(model), linearEquationSolver(storm::utility::solver::getLinearEquationSolver<ValueType>()) {
+            // Intentionally left empty.
+        }
+        
+        template<class ValueType>
+        SparseCtmcCslModelChecker<ValueType>::SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model, std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>&& linearEquationSolver) : SparsePropositionalModelChecker<ValueType>(model), linearEquationSolver(std::move(linearEquationSolver)) {
+            // Intentionally left empty.
+        }
+        
+        template<class ValueType>
+        bool SparseCtmcCslModelChecker<ValueType>::canHandle(storm::logic::Formula const& formula) const {
+            // FIXME: refine.
+            return formula.isCslStateFormula() || formula.isCslPathFormula() || formula.isRewardPathFormula();
+        }
+        
+        template<class ValueType>
+        std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<ValueType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(pathFormula.isIntervalBounded(), storm::exceptions::InvalidPropertyException, "Cannot treat non-interval bounded until.");
+            
+            std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
+            std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
+            ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();;
+            ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();
+            std::pair<double, double> const& intervalBounds =  pathFormula.getIntervalBounds();
+            std::unique_ptr<CheckResult> result = std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeBoundedUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, intervalBounds.first, intervalBounds.second)));
+            
+            return result;
+        }
+        
+        template<class ValueType>
+        std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<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());
+            ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
+            std::vector<ValueType> result = SparseDtmcPrctlModelChecker<ValueType>::computeNextProbabilitiesHelper(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolver);
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(result)));
+        }
+        
+        template<class ValueType>
+        std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<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());
+            ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();
+            ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolver)));
+        }
+        
+        template<class ValueType>
+        storm::models::sparse::Ctmc<ValueType> const& SparseCtmcCslModelChecker<ValueType>::getModel() const {
+            return this->template getModelAs<storm::models::sparse::Ctmc<ValueType>>();
+        }
+        
+        template<class ValueType>
+        std::vector<ValueType> SparseCtmcCslModelChecker<ValueType>::computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& exitRates, bool qualitative, double lowerBound, double upperBound) const {
+            // If the time bounds are [0, inf], we rather call untimed reachability.
+            storm::utility::ConstantsComparator<ValueType> comparator;
+            if (comparator.isZero(lowerBound) && comparator.isInfinity(upperBound)) {
+                return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), phiStates, psiStates, qualitative, *this->linearEquationSolver)));
+            }
+
+            // From this point on, we know that we have to solve a more complicated problem [t, t'] with either t != 0
+            // or t' != inf.
+            
+            // Create the result vector.
+            std::vector<ValueType> result(this->getModel().getNumberOfStates(), storm::utility::zero<ValueType>());
+            
+            // If we identify the states that have probability 0 of reaching the target states, we can exclude them from the
+            // further computations.
+            storm::storage::BitVector statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(this->getModel(), this->getModel().getBackwardTransitions(), phiStates, psiStates);
+            STORM_LOG_INFO("Found " << statesWithProbabilityGreater0.getNumberOfSetBits() << " states with probability greater 0.");
+            storm::storage::BitVector statesWithProbabilityGreater0NonPsi = statesWithProbabilityGreater0 & ~psiStates;
+            STORM_LOG_INFO("Found " << statesWithProbabilityGreater0NonPsi.getNumberOfSetBits() << " 'maybe' states.");
+            
+            if (!statesWithProbabilityGreater0NonPsi.empty()) {
+                if (comparator.isZero(upperBound)) {
+                    // In this case, the interval is of the form [0, 0].
+                    storm::utility::vector::setVectorValues<ValueType>(result, psiStates, storm::utility::one<ValueType>());
+                } else {
+                    // Find the maximal rate of all 'maybe' states to take it as the uniformization rate.
+                    ValueType uniformizationRate = 0;
+                    for (auto const& state : statesWithProbabilityGreater0NonPsi) {
+                        uniformizationRate = std::max(uniformizationRate, exitRates[state]);
+                    }
+                    STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                    
+                    if (comparator.isZero(lowerBound)) {
+                        // In this case, the interval is of the form [0, t].
+                        // Note that this excludes [0, inf] since this is untimed reachability and we considered this case earlier.
+
+                        // Compute the uniformized matrix.
+                        storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), statesWithProbabilityGreater0, psiStates, uniformizationRate, exitRates);
+                        
+                    } else if (comparator.isInfinity(upperBound)) {
+                        // In this case, the interval is of the form [t, inf] with t != 0.
+                        
+                    } else {
+                        // In this case, the interval is of the form [t, t'] with t != 0 and t' != inf.
+                        
+                    }
+                }
+            }
+            
+            return result;
+        }
+        
+        template<class ValueType>
+        storm::storage::SparseMatrix<ValueType> SparseCtmcCslModelChecker<ValueType>::computeUniformizedMatrix(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& maybeStates, storm::storage::BitVector const& absorbingStates, ValueType uniformizationRate, std::vector<ValueType> const& exitRates) {
+            // Create the submatrix that only contains the states with a positive probability (including the
+            // psi states) and reserve space for elements on the diagonal.
+            storm::storage::SparseMatrix<ValueType> uniformizedMatrix = transitionMatrix.getSubmatrix(false, maybeStates, maybeStates, true);
+            
+            // Make the appropriate states absorbing.
+            uniformizedMatrix.makeRowsAbsorbing(absorbingStates % maybeStates);
+            
+            // Now we need to perform the actual uniformization. That is, all entries need to be divided by
+            // the uniformization rate, and the diagonal needs to be set to the negative exit rate of the
+            // state plus the self-loop rate and then increased by one.
+            uint_fast64_t currentRow = 0;
+            for (auto const& state : maybeStates) {
+                for (auto& element : uniformizedMatrix.getRow(currentRow)) {
+                    if (element.getColumn() == currentRow) {
+                        if (absorbingStates.get(state)) {
+                            // Nothing to do here, since the state has already been made absorbing.
+                        } else {
+                            element.setValue(-(exitRates[state] + element.getValue()) / uniformizationRate + storm::utility::one<ValueType>());
+                        }
+                    } else {
+                        element.setValue(element.getValue() / uniformizationRate);
+                    }
+                }
+                ++currentRow;
+            }
+        }
+        
+        template<class ValueType>
+        std::vector<ValueType> SparseCtmcCslModelChecker<ValueType>::computeUntilProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver) {
+            return SparseDtmcPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(transitionMatrix, backwardTransitions, phiStates, psiStates, qualitative, linearEquationSolver);
+        }
+        
+    } // namespace modelchecker
+} // namespace storm
\ No newline at end of file
diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/modelchecker/csl/SparseCtmcCslModelChecker.h
new file mode 100644
index 000000000..a88c36173
--- /dev/null
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.h
@@ -0,0 +1,49 @@
+#ifndef STORM_MODELCHECKER_SPARSECTMCCSLMODELCHECKER_H_
+#define STORM_MODELCHECKER_SPARSECTMCCSLMODELCHECKER_H_
+
+#include "src/modelchecker/propositional/SparsePropositionalModelChecker.h"
+#include "src/models/sparse/Ctmc.h"
+#include "src/solver/LinearEquationSolver.h"
+
+namespace storm {
+    namespace modelchecker {
+        
+        template<class ValueType>
+        class SparseCtmcCslModelChecker : public SparsePropositionalModelChecker<ValueType> {
+        public:
+            explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model);
+            explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model, std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>&& linearEquationSolver);
+            
+            // 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;
+            
+        protected:
+            storm::models::sparse::Ctmc<ValueType> const& getModel() const override;
+            
+        private:
+            // The methods that perform the actual checking.
+            std::vector<ValueType> computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& exitRates, bool qualitative, double lowerBound, double upperBound) const;
+            static std::vector<ValueType> computeUntilProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver);
+            
+            /*!
+             * Computes the matrix representing the transitions of the uniformized CTMC.
+             *
+             * @param transitionMatrix The matrix to uniformize.
+             * @param maybeStates The states that need to be considered.
+             * @param absorbingStates The states that need to be made absorbing.
+             * @param uniformizationRate The rate to be used for uniformization.
+             * @param exitRates The exit rates of all states.
+             */
+            storm::storage::SparseMatrix<ValueType> computeUniformizedMatrix(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& maybeStates, storm::storage::BitVector const& absorbingStates, ValueType uniformizationRate, std::vector<ValueType> const& exitRates);
+            
+            // An object that is used for solving linear equations and performing matrix-vector multiplication.
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> linearEquationSolver;
+        };
+        
+    } // namespace modelchecker
+} // namespace storm
+
+#endif /* STORM_MODELCHECKER_SPARSECTMCCSLMODELCHECKER_H_ */
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
index 7dbc8c9b5..a4e9d7ae9 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
@@ -5,6 +5,7 @@
 #include "src/utility/macros.h"
 #include "src/utility/vector.h"
 #include "src/utility/graph.h"
+#include "src/utility/solver.h"
 
 #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
@@ -75,14 +76,13 @@ namespace storm {
         }
         
         template<typename ValueType>
-        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeNextProbabilitiesHelper(storm::storage::BitVector const& nextStates) {
+        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeNextProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver) {
             // Create the vector with which to multiply and initialize it correctly.
-            std::vector<ValueType> result(this->getModel().getNumberOfStates());
+            std::vector<ValueType> result(transitionMatrix.getRowCount());
             storm::utility::vector::setVectorValues(result, nextStates, storm::utility::one<ValueType>());
             
             // Perform one single matrix-vector multiplication.
-            STORM_LOG_THROW(linearEquationSolver != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
-            this->linearEquationSolver->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), result);
+            linearEquationSolver.performMatrixVectorMultiplication(transitionMatrix, result);
             return result;
         }
         
@@ -90,14 +90,14 @@ namespace storm {
         std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<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());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeNextProbabilitiesHelper(subResult.getTruthValuesVector())));
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeNextProbabilitiesHelper(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolver)));
         }
         
         template<typename ValueType>
-        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const {
+        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver) {
             // We need to identify the states which have to be taken out of the matrix, i.e.
             // all states that have probability 0 and 1 of satisfying the until-formula.
-            std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(this->getModel(), phiStates, psiStates);
+            std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(backwardTransitions, phiStates, psiStates);
             storm::storage::BitVector statesWithProbability0 = std::move(statesWithProbability01.first);
             storm::storage::BitVector statesWithProbability1 = std::move(statesWithProbability01.second);
             
@@ -108,7 +108,7 @@ namespace storm {
             STORM_LOG_INFO("Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
             
             // Create resulting vector.
-            std::vector<ValueType> result(this->getModel().getNumberOfStates());
+            std::vector<ValueType> result(transitionMatrix.getRowCount());
             
             // Check whether we need to compute exact probabilities for some states.
             if (qualitative) {
@@ -119,7 +119,7 @@ namespace storm {
                     // In this case we have have to compute the probabilities.
                     
                     // We can eliminate the rows and columns from the original transition probability matrix.
-                    storm::storage::SparseMatrix<ValueType> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(true, maybeStates, maybeStates, true);
+                    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.
@@ -132,11 +132,10 @@ namespace storm {
                     
                     // Prepare the right-hand side of the equation system. For entry i this corresponds to
                     // the accumulated probability of going from state i to some 'yes' state.
-                    std::vector<ValueType> b = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1);
+                    std::vector<ValueType> b = transitionMatrix.getConstrainedRowSumVector(maybeStates, statesWithProbability1);
                     
                     // Now solve the created system of linear equations.
-                    STORM_LOG_THROW(linearEquationSolver != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
-                    this->linearEquationSolver->solveEquationSystem(submatrix, x, b);
+                    linearEquationSolver.solveEquationSystem(submatrix, x, b);
                     
                     // Set values of resulting vector according to result.
                     storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
@@ -156,7 +155,7 @@ namespace storm {
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();;
             ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();;
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative)));
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolver)));
         }
         
         template<typename ValueType>
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
index 967ac5d8c..7133c12c8 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
@@ -8,6 +8,8 @@
 
 namespace storm {
     namespace modelchecker {
+        // Forward-declare CTMC model checker so we can make it a friend.
+        template<typename ValueType> class SparseCtmcCslModelChecker;
         
         template<class ValueType>
         class SparseDtmcPrctlModelChecker : public SparsePropositionalModelChecker<ValueType> {
@@ -30,8 +32,8 @@ namespace storm {
         private:
             // The methods that perform the actual checking.
             std::vector<ValueType> computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound) const;
-            std::vector<ValueType> computeNextProbabilitiesHelper(storm::storage::BitVector const& nextStates);
-            std::vector<ValueType> computeUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const;
+            static std::vector<ValueType> computeNextProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver);
+            static std::vector<ValueType> computeUntilProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver);
             std::vector<ValueType> computeInstantaneousRewardsHelper(uint_fast64_t stepCount) const;
             std::vector<ValueType> computeCumulativeRewardsHelper(uint_fast64_t stepBound) const;
             std::vector<ValueType> computeReachabilityRewardsHelper(storm::storage::BitVector const& targetStates, bool qualitative) const;
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
index 5078cd8d8..8f30f22b2 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
@@ -70,7 +70,7 @@ namespace storm {
         
         template<typename ValueType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<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 this->getModel().");
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
             std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();
@@ -93,7 +93,7 @@ namespace storm {
         
         template<typename ValueType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<ValueType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
-            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic this->getModel().");
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
             std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeNextProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector())));
@@ -162,7 +162,7 @@ namespace storm {
         
         template<typename ValueType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<ValueType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
-            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic this->getModel().");
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
             std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();
@@ -201,7 +201,7 @@ namespace storm {
         
         template<typename ValueType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, 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 this->getModel().");
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeCumulativeRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, rewardPathFormula.getStepBound())));
         }
         
@@ -221,7 +221,7 @@ namespace storm {
         
         template<typename ValueType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, 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 this->getModel().");
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeInstantaneousRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, rewardPathFormula.getStepCount())));
         }
         
@@ -308,7 +308,7 @@ namespace storm {
         
         template<typename ValueType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, 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 this->getModel().");
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector(), qualitative)));
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
index 8f51d9003..cfe6188e2 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
@@ -10,6 +10,9 @@ namespace storm {
     namespace counterexamples {
         template<typename ValueType>
         class SMTMinimalCommandSetGenerator;
+        
+        template<typename ValueType>
+        class MILPMinimalLabelSetGenerator;
     }
     
     namespace modelchecker {
@@ -22,7 +25,8 @@ namespace storm {
         class SparseMdpPrctlModelChecker : public SparsePropositionalModelChecker<ValueType> {
         public:
             friend class SparseMarkovAutomatonCslModelChecker<ValueType>;
-            friend class counterexamples::SMTMinimalCommandSetGenerator<ValueType>;
+            friend class storm::counterexamples::SMTMinimalCommandSetGenerator<ValueType>;
+            friend class storm::counterexamples::MILPMinimalLabelSetGenerator<ValueType>;
             
             explicit SparseMdpPrctlModelChecker(storm::models::sparse::Mdp<ValueType> const& model);
             explicit SparseMdpPrctlModelChecker(storm::models::sparse::Mdp<ValueType> const& model, std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> nondeterministicLinearEquationSolver);
diff --git a/src/modelchecker/results/ExplicitQuantitativeCheckResult.cpp b/src/modelchecker/results/ExplicitQuantitativeCheckResult.cpp
index 6c5b4520f..f54984ac3 100644
--- a/src/modelchecker/results/ExplicitQuantitativeCheckResult.cpp
+++ b/src/modelchecker/results/ExplicitQuantitativeCheckResult.cpp
@@ -116,28 +116,28 @@ namespace storm {
                 switch (comparisonType) {
                     case logic::Less:
                         for (uint_fast64_t index = 0; index < valuesAsVector.size(); ++index) {
-                            if (result[index] < bound) {
+                            if (valuesAsVector[index] < bound) {
                                 result.set(index);
                             }
                         }
                         break;
                     case logic::LessEqual:
                         for (uint_fast64_t index = 0; index < valuesAsVector.size(); ++index) {
-                            if (result[index] <= bound) {
+                            if (valuesAsVector[index] <= bound) {
                                 result.set(index);
                             }
                         }
                         break;
                     case logic::Greater:
                         for (uint_fast64_t index = 0; index < valuesAsVector.size(); ++index) {
-                            if (result[index] > bound) {
+                            if (valuesAsVector[index] > bound) {
                                 result.set(index);
                             }
                         }
                         break;
                     case logic::GreaterEqual:
                         for (uint_fast64_t index = 0; index < valuesAsVector.size(); ++index) {
-                            if (result[index] >= bound) {
+                            if (valuesAsVector[index] >= bound) {
                                 result.set(index);
                             }
                         }
diff --git a/src/storage/DeterministicModelBisimulationDecomposition.cpp b/src/storage/DeterministicModelBisimulationDecomposition.cpp
index 3563031af..85a66d02b 100644
--- a/src/storage/DeterministicModelBisimulationDecomposition.cpp
+++ b/src/storage/DeterministicModelBisimulationDecomposition.cpp
@@ -1493,4 +1493,4 @@ namespace storm {
         template class DeterministicModelBisimulationDecomposition<storm::RationalFunction>;
 #endif
     }
-}
\ No newline at end of file
+}
diff --git a/src/storage/expressions/BinaryNumericalFunctionExpression.cpp b/src/storage/expressions/BinaryNumericalFunctionExpression.cpp
index be0f30dc0..4e95ea551 100644
--- a/src/storage/expressions/BinaryNumericalFunctionExpression.cpp
+++ b/src/storage/expressions/BinaryNumericalFunctionExpression.cpp
@@ -6,6 +6,7 @@
 #include "src/storage/expressions/DoubleLiteralExpression.h"
 #include "src/utility/macros.h"
 #include "src/exceptions/InvalidTypeException.h"
+#include "src/exceptions/InvalidStateException.h"
 
 namespace storm {
     namespace expressions {
@@ -65,7 +66,7 @@ namespace storm {
             std::shared_ptr<BaseExpression const> firstOperandSimplified = this->getFirstOperand()->simplify();
             std::shared_ptr<BaseExpression const> secondOperandSimplified = this->getSecondOperand()->simplify();
             
-            if (firstOperandSimplified->isLiteral() && secondOperandSimplified->isLiteral()) {
+            if (firstOperandSimplified->isLiteral() && secondOperandSimplified->isLiteral() && this->getOperatorType() != OperatorType::Divide) {
                 if (this->hasIntegerType()) {
                     int_fast64_t firstOperandEvaluation = firstOperandSimplified->evaluateAsInt();
                     int_fast64_t secondOperandEvaluation = secondOperandSimplified->evaluateAsInt();
@@ -74,10 +75,10 @@ namespace storm {
                         case OperatorType::Plus: newValue = firstOperandEvaluation + secondOperandEvaluation; break;
                         case OperatorType::Minus: newValue = firstOperandEvaluation - secondOperandEvaluation; break;
                         case OperatorType::Times: newValue = firstOperandEvaluation * secondOperandEvaluation; break;
-                        case OperatorType::Divide: newValue = firstOperandEvaluation / secondOperandEvaluation; break;
                         case OperatorType::Min: newValue = std::min(firstOperandEvaluation, secondOperandEvaluation); break;
                         case OperatorType::Max: newValue = std::max(firstOperandEvaluation, secondOperandEvaluation); break;
                         case OperatorType::Power: newValue = static_cast<int_fast64_t>(std::pow(firstOperandEvaluation, secondOperandEvaluation)); break;
+                        case OperatorType::Divide: STORM_LOG_THROW(false, storm::exceptions::InvalidStateException, "Unable to simplify division."); break;
                     }
                     return std::shared_ptr<BaseExpression>(new IntegerLiteralExpression(this->getManager(), newValue));
                 } else if (this->hasRationalType()) {
@@ -88,10 +89,10 @@ namespace storm {
                         case OperatorType::Plus: newValue = firstOperandEvaluation + secondOperandEvaluation; break;
                         case OperatorType::Minus: newValue = firstOperandEvaluation - secondOperandEvaluation; break;
                         case OperatorType::Times: newValue = firstOperandEvaluation * secondOperandEvaluation; break;
-                        case OperatorType::Divide: newValue = firstOperandEvaluation / secondOperandEvaluation; break;
                         case OperatorType::Min: newValue = std::min(firstOperandEvaluation, secondOperandEvaluation); break;
                         case OperatorType::Max: newValue = std::max(firstOperandEvaluation, secondOperandEvaluation); break;
                         case OperatorType::Power: newValue = static_cast<int_fast64_t>(std::pow(firstOperandEvaluation, secondOperandEvaluation)); break;
+                        case OperatorType::Divide: STORM_LOG_THROW(false, storm::exceptions::InvalidStateException, "Unable to simplify division."); break;
                     }
                     return std::shared_ptr<BaseExpression>(new DoubleLiteralExpression(this->getManager(), newValue));
                 }
diff --git a/src/storage/expressions/UnaryNumericalFunctionExpression.cpp b/src/storage/expressions/UnaryNumericalFunctionExpression.cpp
index 6e2a13ae0..edfb22702 100644
--- a/src/storage/expressions/UnaryNumericalFunctionExpression.cpp
+++ b/src/storage/expressions/UnaryNumericalFunctionExpression.cpp
@@ -76,7 +76,7 @@ namespace storm {
                         case OperatorType::Floor: value = std::floor(boost::get<double>(operandEvaluation)); break;
                         case OperatorType::Ceil: value = std::ceil(boost::get<double>(operandEvaluation)); break;
                     }
-                    return std::shared_ptr<BaseExpression>(new DoubleLiteralExpression(this->getManager(), boost::get<double>(result)));
+                    return std::shared_ptr<BaseExpression>(new DoubleLiteralExpression(this->getManager(), value));
                 }
             }
             
diff --git a/src/utility/cli.h b/src/utility/cli.h
index b24e57934..f86acabf8 100644
--- a/src/utility/cli.h
+++ b/src/utility/cli.h
@@ -325,6 +325,11 @@ namespace storm {
                     }
                     options.addConstantDefinitionsFromString(program, settings.getConstantDefinitionString());
                     
+                    // Generate command labels if we are going to build a counterexample later.
+                    if (storm::settings::counterexampleGeneratorSettings().isMinimalCommandSetGenerationSet()) {
+                        options.buildCommandLabels = true;
+                    }
+                    
                     result = storm::builder::ExplicitPrismModelBuilder<ValueType>::translateProgram(program, options);
                 } else if (settings.getEngine() == storm::settings::modules::GeneralSettings::Engine::Dd) {
                     typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
@@ -336,6 +341,7 @@ namespace storm {
                     result = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
                 }
                 
+                // Then, build the model from the symbolic description.
                 return result;
             }
             
diff --git a/src/utility/numerical.h b/src/utility/numerical.h
new file mode 100644
index 000000000..d14036859
--- /dev/null
+++ b/src/utility/numerical.h
@@ -0,0 +1,131 @@
+#include <vector>
+#include <tuple>
+#include <cmath>
+
+#include "src/utility/macros.h"
+#include "src/utility/constants.h"
+#include "src/exceptions/InvalidArgumentException.h"
+
+namespace storm {
+    namespace utility {
+        namespace numerical {
+            template<typename ValueType>
+            std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> getFoxGlynnCutoff(ValueType lambda, ValueType underflow, ValueType overflow, ValueType accuracy) {
+                storm::utility::ConstantsComparator<ValueType> comparator;
+                STORM_LOG_THROW(!comparator.isZero(lambda), storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: lambda must not be zero.");
+                
+                // This code is a modified version of the one in PRISM. According to their implementation, for lambda
+                // smaller than 400, we compute the result using the naive method.
+                if (lambda < 25) {
+                    ValueType eToPowerMinusLambda = std::exp(-lambda);
+                    ValueType targetValue = (1 - accuracy) / eToPowerMinusLambda;
+                    std::vector<ValueType> weights;
+                    
+                    ValueType exactlyKEvents = 1;
+                    ValueType atMostKEvents = exactlyKEvents;
+                    weights.push_back(exactlyKEvents * eToPowerMinusLambda);
+                    
+                    uint_fast64_t k = 1;
+                    do {
+                        exactlyKEvents *= lambda / k;
+                        atMostKEvents += exactlyKEvents;
+                        weights.push_back(exactlyKEvents * eToPowerMinusLambda);
+                        ++k;
+                    } while (atMostKEvents < targetValue);
+                    
+                    return std::make_tuple(0, k - 1, 1.0, weights);
+                } else {
+                    STORM_LOG_THROW(accuracy >= 1e-10, storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: the accuracy must not be below 1e-10.");
+                    
+                    // Factor from Fox&Glynn's paper. The paper does not explain where it comes from.
+                    ValueType factor = 1e+10;
+                    
+                    // Now start the Finder algorithm to find the truncation points.
+                    ValueType m = std::floor(lambda);
+                    uint_fast64_t leftTruncationPoint = 0, rightTruncationPoint = 0;
+                    {
+                        // Factors used by the corollaries explained in Fox & Glynns paper.
+                        // Square root of pi.
+                        ValueType sqrtpi = 1.77245385090551602729;
+                        
+                        // Square root of 2.
+                        ValueType sqrt2 = 1.41421356237309504880;
+                        
+                        // Set up a_\lambda, b_\lambda, and the square root of lambda.
+                        ValueType aLambda = 0, bLambda = 0, sqrtLambda = 0;
+                        if (m < 400) {
+                            sqrtLambda = std::sqrt(400.0);
+                            aLambda = (1.0 + 1.0 / 400.0) * exp(0.0625) * sqrt2;
+                            bLambda = (1.0 + 1.0 / 400.0) * exp(0.125 / 400.0);
+                        } else {
+                            sqrtLambda = std::sqrt(lambda);
+                            aLambda = (1.0 + 1.0 / lambda) * exp(0.0625) * sqrt2;
+                            bLambda = (1.0 + 1.0 / lambda) * exp(0.125 / lambda);
+                        }
+                        
+                        // Use Corollary 1 from the paper to find the right truncation point.
+                        uint_fast64_t k = 3;
+                        
+                        ValueType dkl = 1.0 / (1 - std::exp(-(2.0 / 9.0) * (k * sqrt2 * sqrtLambda + 1.5)));
+                        
+                        // According to David Jansen the upper bound can be ignored to achieve more accurate results.
+                        // Right hand side of the equation in Corollary 1.
+                        while ((accuracy / 2.0) < (aLambda * dkl * std::exp(-(k*k / 2.0)) / (k * sqrt2 * sqrtpi))) {
+                            ++k;
+                            
+                            // d(k,Lambda) from the paper.
+                            dkl = 1.0 / (1 - std::exp(-(2.0 / 9.0)*(k * sqrt2 * sqrtLambda + 1.5)));
+                        }
+                        
+                        // Left hand side of the equation in Corollary 1.
+                        rightTruncationPoint = static_cast<uint_fast64_t>(std::ceil((m + k * sqrt2 * sqrtLambda + 1.5)));
+                        
+                        // Use Corollary 2 to find left truncation point.
+                        k = 3;
+                        
+                        // Right hand side of the equation in Corollary 2.
+                        while ((accuracy / 2.0) < ((bLambda * exp(-(k*k / 2.0))) / (k * sqrt2 * sqrtpi))) {
+                            ++k;
+                        }
+                        
+                        // Left hand side of the equation in Corollary 2.
+                        leftTruncationPoint = std::max(static_cast<uint_fast64_t>(std::trunc(m - k * sqrtLambda - 1.5)), uint_fast64_t(0));
+                        
+                        // Check for underflow.
+                        STORM_LOG_THROW(static_cast<uint_fast64_t>(std::trunc((m - k * sqrtLambda - 1.5))) >= 0, storm::exceptions::OutOfRangeException, "Error in Fox-Glynn algorithm: Underflow of left truncation point.");
+                    }
+                    
+                    std::vector<ValueType> weights(rightTruncationPoint - leftTruncationPoint + 1);
+                    weights[m - leftTruncationPoint] = overflow / (factor * (rightTruncationPoint - leftTruncationPoint));
+                    for (uint_fast64_t j = m; j > leftTruncationPoint; --j) {
+                        weights[j - 1 - leftTruncationPoint] = (j / lambda) * weights[j - leftTruncationPoint];
+                    }
+                    
+                    for (uint_fast64_t j = m; j < rightTruncationPoint; ++j) {
+                        weights[j + 1 - leftTruncationPoint] = (lambda / (j + 1)) * weights[j - leftTruncationPoint];
+                    }
+                    
+                    
+                    // Compute the total weight and start with smallest to avoid roundoff errors.
+                    ValueType totalWeight = storm::utility::zero<ValueType>();
+                    
+                    uint_fast64_t s = leftTruncationPoint;
+                    uint_fast64_t t = rightTruncationPoint;
+                    while (s < t) {
+                        if (weights[s - leftTruncationPoint] <= weights[t - leftTruncationPoint]) {
+                            totalWeight += weights[s - leftTruncationPoint];
+                            ++s;
+                        } else {
+                            totalWeight += weights[t - leftTruncationPoint];
+                            --t;
+                        }
+                    }
+                    
+                    totalWeight += weights[s - leftTruncationPoint];
+                    
+                    return std::make_tuple(leftTruncationPoint, rightTruncationPoint, totalWeight, weights);
+                }
+            }
+        }
+    }
+}
\ No newline at end of file