diff --git a/src/modelchecker/reachability/DirectEncoding.h b/src/modelchecker/reachability/DirectEncoding.h
index ef5b8db3f..69ddeeb28 100644
--- a/src/modelchecker/reachability/DirectEncoding.h
+++ b/src/modelchecker/reachability/DirectEncoding.h
@@ -20,7 +20,7 @@ namespace storm
             {
             public:
                 template<typename T>
-                std::string encodeAsSmt2(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<T> const& oneStepProbabilities, std::vector<carl::Variable> const& parameters, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& finalStates, typename T::CoeffType const& threshold, bool lessequal = false) {
+                std::string encodeAsSmt2(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<T> const& oneStepProbabilities, std::set<carl::Variable> const& parameters, storm::storage::BitVector const& initialStates, typename T::CoeffType const& threshold, bool lessequal = false) {
                     
                     carl::io::WriteTosmt2Stream smt2;
                     uint_fast64_t numberOfStates = transitionMatrix.getRowCount();
@@ -37,25 +37,23 @@ namespace storm
                     for (uint_fast64_t state = 0; state < numberOfStates; ++state) {
                         carl::Variable stateVar = vpool.getFreshVariable("s_" + std::to_string(state));
                         stateVars.push_back(stateVar);
-                        if (!finalStates[state]) {
-                            smt2 << ("state_bound_" + std::to_string(state));
-                            smt2 << carl::io::smt2node::AND;
-                            smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(stateVar), carl::CompareRelation::GT);
-                            smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(stateVar) - Polynomial::PolyType(1), carl::CompareRelation::LT);
-                            smt2 << carl::io::smt2node::CLOSENODE;
-                        }
+                        smt2 << ("state_bound_" + std::to_string(state));
+                        smt2 << carl::io::smt2node::AND;
+                        smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(stateVar), carl::CompareRelation::GT);
+                        smt2 << carl::Constraint<Polynomial::PolyType>(Polynomial::PolyType(stateVar) - Polynomial::PolyType(1), carl::CompareRelation::LT);
+                        smt2 << carl::io::smt2node::CLOSENODE;
                     }
                     
                     smt2.setAutomaticLineBreaks(true);
                     Polynomial::PolyType initStateReachSum;
                     for (uint_fast64_t state = 0; state < numberOfStates; ++state) {
-                        T reachpropPol(0);
+                        T reachpropPol;
                         for (auto const& transition : transitionMatrix.getRow(state)) {
-                            reachpropPol += transition.getValue() * stateVars[transition.getColumn()];
+//                            reachpropPol += transition.getValue() * stateVars[transition.getColumn()];
                         }
                         reachpropPol += oneStepProbabilities[state];
                         smt2 << ("transition_" + std::to_string(state));
-                        smt2 << carl::Constraint<T>(reachpropPol - stateVars[state], carl::CompareRelation::EQ);
+//                        smt2 << carl::Constraint<Polynomial::PolyType>(reachpropPol - stateVars[state], carl::CompareRelation::EQ);
                     }
                     
                     smt2 << ("reachability");
diff --git a/src/modelchecker/reachability/SparseSccModelChecker.cpp b/src/modelchecker/reachability/SparseSccModelChecker.cpp
index 51ba38c64..1284af269 100644
--- a/src/modelchecker/reachability/SparseSccModelChecker.cpp
+++ b/src/modelchecker/reachability/SparseSccModelChecker.cpp
@@ -10,37 +10,60 @@
 #include "src/utility/graph.h"
 #include "src/utility/vector.h"
 #include "src/utility/macros.h"
+#include "src/utility/ConstantsComparator.h"
 
 namespace storm {
     namespace modelchecker {
         namespace reachability {
             
             template<typename ValueType>
-            static ValueType&& simplify(ValueType&& value) {
-                // In the general case, we don't to anything here, but merely return the value. If something else is
-                // supposed to happen here, the templated function can be specialized for this particular type.
-                return std::forward<ValueType>(value);
-            }
-            
-            static RationalFunction&& simplify(RationalFunction&& value) {
-                // In the general case, we don't to anything here, but merely return the value. If something else is
-                // supposed to happen here, the templated function can be specialized for this particular type.
-                STORM_LOG_TRACE("Simplifying " << value << " ... ");
-                value.simplify();
-                STORM_LOG_TRACE("done.");
-                return std::forward<RationalFunction>(value);
-            }
-            
-            template<typename IndexType, typename ValueType>
-            static storm::storage::MatrixEntry<IndexType, ValueType>&& simplify(storm::storage::MatrixEntry<IndexType, ValueType>&& matrixEntry) {
-                simplify(matrixEntry.getValue());
-                return std::move(matrixEntry);
-            }
-            
-            template<typename IndexType, typename ValueType>
-            static storm::storage::MatrixEntry<IndexType, ValueType>& simplify(storm::storage::MatrixEntry<IndexType, ValueType>& matrixEntry) {
-                matrixEntry.setValue(simplify(matrixEntry.getValue()));
-                return matrixEntry;
+            ValueType SparseSccModelChecker<ValueType>::computeReachabilityProbability(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType>& oneStepProbabilities, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, boost::optional<std::vector<std::size_t>> const& distances) {
+                auto totalTimeStart = std::chrono::high_resolution_clock::now();
+                auto conversionStart = std::chrono::high_resolution_clock::now();
+
+                // Create a bit vector that represents the subsystem of states we still have to eliminate.
+                storm::storage::BitVector subsystem = storm::storage::BitVector(transitionMatrix.getRowCount(), true);
+                
+                // Then, we convert the reduced matrix to a more flexible format to be able to perform state elimination more easily.
+                FlexibleSparseMatrix<ValueType> flexibleMatrix = getFlexibleSparseMatrix(transitionMatrix);
+                FlexibleSparseMatrix<ValueType> flexibleBackwardTransitions = getFlexibleSparseMatrix(backwardTransitions, true);
+                auto conversionEnd = std::chrono::high_resolution_clock::now();
+                
+                // Then, we recursively treat all SCCs.
+                auto modelCheckingStart = std::chrono::high_resolution_clock::now();
+                std::vector<storm::storage::sparse::state_type> entryStateQueue;
+                uint_fast64_t maximalDepth = treatScc(flexibleMatrix, oneStepProbabilities, initialStates, subsystem, transitionMatrix, flexibleBackwardTransitions, false, 0, storm::settings::parametricSettings().getMaximalSccSize(), entryStateQueue, distances);
+                
+                // If the entry states were to be eliminated last, we need to do so now.
+                STORM_LOG_DEBUG("Eliminating " << entryStateQueue.size() << " entry states as a last step.");
+                if (storm::settings::parametricSettings().isEliminateEntryStatesLastSet()) {
+                    for (auto const& state : entryStateQueue) {
+                        eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions);
+                    }
+                }
+                
+                auto modelCheckingEnd = std::chrono::high_resolution_clock::now();
+                auto totalTimeEnd = std::chrono::high_resolution_clock::now();
+                
+                if (storm::settings::generalSettings().isShowStatisticsSet()) {
+                    auto conversionTime = conversionEnd - conversionStart;
+                    auto modelCheckingTime = modelCheckingEnd - modelCheckingStart;
+                    auto totalTime = totalTimeEnd - totalTimeStart;
+                    
+                    STORM_PRINT_AND_LOG(std::endl);
+                    STORM_PRINT_AND_LOG("Time breakdown:" << std::endl);
+                    STORM_PRINT_AND_LOG("    * time for conversion: " << std::chrono::duration_cast<std::chrono::milliseconds>(conversionTime).count() << "ms" << std::endl);
+                    STORM_PRINT_AND_LOG("    * time for checking: " << std::chrono::duration_cast<std::chrono::milliseconds>(modelCheckingTime).count() << "ms" << std::endl);
+                    STORM_PRINT_AND_LOG("------------------------------------------" << std::endl);
+                    STORM_PRINT_AND_LOG("    * total time: " << std::chrono::duration_cast<std::chrono::milliseconds>(totalTime).count() << "ms" << std::endl);
+                    STORM_PRINT_AND_LOG(std::endl);
+                    STORM_PRINT_AND_LOG("Other:" << std::endl);
+                    STORM_PRINT_AND_LOG("    * number of states eliminated: " << transitionMatrix.getRowCount() << std::endl);
+                    STORM_PRINT_AND_LOG("    * maximal depth of SCC decomposition: " << maximalDepth << std::endl);
+                }
+                
+                // Now, we return the value for the only initial state.
+                return storm::utility::simplify(oneStepProbabilities[*initialStates.begin()]);
             }
             
             template<typename ValueType>
@@ -74,30 +97,26 @@ namespace storm {
                 
                 // Do some sanity checks to establish some required properties.
                 STORM_LOG_THROW(dtmc.getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::IllegalArgumentException, "Input model is required to have exactly one initial state.");
-                storm::storage::sparse::state_type initialStateIndex = *dtmc.getInitialStates().begin();
                 
                 // Then, compute the subset of states that has a probability of 0 or 1, respectively.
-                auto totalTimeStart = std::chrono::high_resolution_clock::now();
-                auto preprocessingStart = std::chrono::high_resolution_clock::now();
                 std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(dtmc, phiStates, psiStates);
                 storm::storage::BitVector statesWithProbability0 = statesWithProbability01.first;
                 storm::storage::BitVector statesWithProbability1 = statesWithProbability01.second;
                 storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
                 
                 // If the initial state is known to have either probability 0 or 1, we can directly return the result.
-                if (!maybeStates.get(initialStateIndex)) {
-                    return statesWithProbability0.get(initialStateIndex) ? storm::utility::constantZero<ValueType>() : storm::utility::constantOne<ValueType>();
+                if (dtmc.getInitialStates().isDisjointFrom(maybeStates)) {
+                    STORM_LOG_DEBUG("The probability of all initial states was found in a preprocessing step.");
+                    return statesWithProbability0.get(*dtmc.getInitialStates().begin()) ? storm::utility::constantZero<ValueType>() : storm::utility::constantOne<ValueType>();
                 }
-                
+
                 // Determine the set of states that is reachable from the initial state without jumping over a target state.
                 storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(dtmc.getTransitionMatrix(), dtmc.getInitialStates(), maybeStates, statesWithProbability1);
                 
                 // Subtract from the maybe states the set of states that is not reachable (on a path from the initial to a target state).
                 maybeStates &= reachableStates;
-                auto preprocessingEnd = std::chrono::high_resolution_clock::now();
-                
+
                 // Create a vector for the probabilities to go to a state with probability 1 in one step.
-                auto conversionStart = std::chrono::high_resolution_clock::now();
                 std::vector<ValueType> oneStepProbabilities = dtmc.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1);
                 
                 // Determine the set of initial states of the sub-DTMC.
@@ -108,9 +127,10 @@ namespace storm {
                 
                 // To be able to apply heuristics later, we now determine the distance of each state to the initial state.
                 std::vector<std::pair<storm::storage::sparse::state_type, std::size_t>> stateQueue;
-                stateQueue.reserve(maybeStates.getNumberOfSetBits());
-                std::vector<std::size_t> distances(maybeStates.getNumberOfSetBits());
-                storm::storage::BitVector statesInQueue(maybeStates.getNumberOfSetBits());
+                stateQueue.reserve(submatrix.getRowCount());
+                storm::storage::BitVector statesInQueue(submatrix.getRowCount());
+                std::vector<std::size_t> distances(submatrix.getRowCount());
+                
                 storm::storage::sparse::state_type currentPosition = 0;
                 for (auto const& initialState : newInitialStates) {
                     stateQueue.emplace_back(initialState, 0);
@@ -131,51 +151,7 @@ namespace storm {
                     ++currentPosition;
                 }
                 
-                // Create a bit vector that represents the subsystem of states we still have to eliminate.
-                storm::storage::BitVector subsystem = storm::storage::BitVector(maybeStates.getNumberOfSetBits(), true);
-                
-                // Then, we convert the reduced matrix to a more flexible format to be able to perform state elimination more easily.
-                FlexibleSparseMatrix<ValueType> flexibleMatrix = getFlexibleSparseMatrix(submatrix);
-                FlexibleSparseMatrix<ValueType> flexibleBackwardTransitions = getFlexibleSparseMatrix(submatrix.transpose(), true);
-                auto conversionEnd = std::chrono::high_resolution_clock::now();
-                
-                // Then, we recursively treat all SCCs.
-                auto modelCheckingStart = std::chrono::high_resolution_clock::now();
-                std::vector<storm::storage::sparse::state_type> entryStateQueue;
-                uint_fast64_t maximalDepth = treatScc(dtmc, flexibleMatrix, oneStepProbabilities, newInitialStates, subsystem, submatrix, flexibleBackwardTransitions, false, 0, storm::settings::parametricSettings().getMaximalSccSize(), entryStateQueue, distances);
-
-                // If the entry states were to be eliminated last, we need to do so now.
-                STORM_LOG_DEBUG("Eliminating " << entryStateQueue.size() << " entry states as a last step.");
-                if (storm::settings::parametricSettings().isEliminateEntryStatesLastSet()) {
-                    for (auto const& state : entryStateQueue) {
-                        eliminateState(flexibleMatrix, oneStepProbabilities, state, flexibleBackwardTransitions);
-                    }
-                }
-                
-                auto modelCheckingEnd = std::chrono::high_resolution_clock::now();
-                auto totalTimeEnd = std::chrono::high_resolution_clock::now();
-                
-                if (storm::settings::generalSettings().isShowStatisticsSet()) {
-                    auto preprocessingTime = preprocessingEnd - preprocessingStart;
-                    auto conversionTime = conversionEnd - conversionStart;
-                    auto modelCheckingTime = modelCheckingEnd - modelCheckingStart;
-                    auto totalTime = totalTimeEnd - totalTimeStart;
-                    
-                    STORM_PRINT_AND_LOG(std::endl);
-                    STORM_PRINT_AND_LOG("Time breakdown:" << std::endl);
-                    STORM_PRINT_AND_LOG("    * time for preprocessing: " << std::chrono::duration_cast<std::chrono::milliseconds>(preprocessingTime).count() << "ms" << std::endl);
-                    STORM_PRINT_AND_LOG("    * time for conversion: " << std::chrono::duration_cast<std::chrono::milliseconds>(conversionTime).count() << "ms" << std::endl);
-                    STORM_PRINT_AND_LOG("    * time for checking: " << std::chrono::duration_cast<std::chrono::milliseconds>(modelCheckingTime).count() << "ms" << std::endl);
-                    STORM_PRINT_AND_LOG("------------------------------------------" << std::endl);
-                    STORM_PRINT_AND_LOG("    * total time: " << std::chrono::duration_cast<std::chrono::milliseconds>(totalTime).count() << "ms" << std::endl);
-                    STORM_PRINT_AND_LOG(std::endl);
-                    STORM_PRINT_AND_LOG("Other:" << std::endl);
-                    STORM_PRINT_AND_LOG("    * number of states eliminated: " << maybeStates.getNumberOfSetBits() << std::endl);
-                    STORM_PRINT_AND_LOG("    * maximal depth of SCC decomposition: " << maximalDepth << std::endl);
-                }
-                
-                // Now, we return the value for the only initial state.
-                return simplify(oneStepProbabilities[initialStateIndex]);
+                return computeReachabilityProbability(submatrix, oneStepProbabilities, submatrix.transpose(), newInitialStates, phiStates, psiStates, distances);
             }
             
             template<typename ValueType>
@@ -193,7 +169,7 @@ namespace storm {
             }
             
             template<typename ValueType>
-            uint_fast64_t SparseSccModelChecker<ValueType>::treatScc(storm::models::Dtmc<ValueType> const& dtmc, FlexibleSparseMatrix<ValueType>& matrix, std::vector<ValueType>& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix<ValueType> const& forwardTransitions, FlexibleSparseMatrix<ValueType>& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level, uint_fast64_t maximalSccSize, std::vector<storm::storage::sparse::state_type>& entryStateQueue, std::vector<std::size_t> const& distances) {
+            uint_fast64_t SparseSccModelChecker<ValueType>::treatScc(FlexibleSparseMatrix<ValueType>& matrix, std::vector<ValueType>& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix<ValueType> const& forwardTransitions, FlexibleSparseMatrix<ValueType>& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level, uint_fast64_t maximalSccSize, std::vector<storm::storage::sparse::state_type>& entryStateQueue, boost::optional<std::vector<std::size_t>> const& distances) {
                 uint_fast64_t maximalDepth = level;
                 
                 // If the SCCs are large enough, we try to split them further.
@@ -220,7 +196,9 @@ namespace storm {
                     
                     STORM_LOG_DEBUG("Eliminating " << trivialSccs.size() << " trivial SCCs.");
                     if (storm::settings::parametricSettings().isSortTrivialSccsSet()) {
-                        std::sort(trivialSccs.begin(), trivialSccs.end(), [&distances] (std::pair<storm::storage::sparse::state_type, uint_fast64_t> const& state1, std::pair<storm::storage::sparse::state_type, uint_fast64_t> const& state2) -> bool { return distances[state1.first] > distances[state2.first]; } );
+                        STORM_LOG_THROW(distances, storm::exceptions::IllegalFunctionCallException, "Cannot sort according to distances because none were provided.");
+                        std::vector<std::size_t> const& actualDistances = distances.get();
+                        std::sort(trivialSccs.begin(), trivialSccs.end(), [&actualDistances] (std::pair<storm::storage::sparse::state_type, uint_fast64_t> const& state1, std::pair<storm::storage::sparse::state_type, uint_fast64_t> const& state2) -> bool { return actualDistances[state1.first] > actualDistances[state2.first]; } );
                     }
                     for (auto const& stateIndexPair : trivialSccs) {
                         eliminateState(matrix, oneStepProbabilities, stateIndexPair.first, backwardTransitions);
@@ -237,7 +215,7 @@ namespace storm {
                         storm::storage::BitVector newSccAsBitVector(forwardTransitions.getRowCount(), newScc.begin(), newScc.end());
                         
                         // Determine the set of entry states of the SCC.
-                        storm::storage::BitVector entryStates(dtmc.getNumberOfStates());
+                        storm::storage::BitVector entryStates(forwardTransitions.getRowCount());
                         for (auto const& state : newScc) {
                             for (auto const& predecessor : backwardTransitions.getRow(state)) {
                                 if (predecessor.getValue() != storm::utility::constantZero<ValueType>() && !newSccAsBitVector.get(predecessor.getColumn())) {
@@ -247,7 +225,7 @@ namespace storm {
                         }
                         
                         // Recursively descend in SCC-hierarchy.
-                        uint_fast64_t depth = treatScc(dtmc, matrix, oneStepProbabilities, entryStates, newSccAsBitVector, forwardTransitions, backwardTransitions, !storm::settings::parametricSettings().isEliminateEntryStatesLastSet(), level + 1, maximalSccSize, entryStateQueue, distances);
+                        uint_fast64_t depth = treatScc(matrix, oneStepProbabilities, entryStates, newSccAsBitVector, forwardTransitions, backwardTransitions, !storm::settings::parametricSettings().isEliminateEntryStatesLastSet(), level + 1, maximalSccSize, entryStateQueue, distances);
                         maximalDepth = std::max(maximalDepth, depth);
                     }
                     
@@ -335,7 +313,7 @@ namespace storm {
                 std::size_t scaledSuccessors = 0;
                 if (hasSelfLoop) {
                     loopProbability = storm::utility::constantOne<ValueType>() / (storm::utility::constantOne<ValueType>() - loopProbability);
-                    simplify(loopProbability);
+                    storm::utility::simplify(loopProbability);
                     for (auto& entry : matrix.getRow(state)) {
                         // Only scale the non-diagonal entries.
                         if (entry.getColumn() != state) {
@@ -409,7 +387,7 @@ namespace storm {
                             auto tmpResult = *first2 * multiplyFactor;
                             multiplicationTime += std::chrono::high_resolution_clock::now() - multiplicationClock;
                             simplifyClock = std::chrono::high_resolution_clock::now();
-                            *result = simplify(std::move(tmpResult));
+                            *result = storm::utility::simplify(tmpResult);
                             simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock;
                             ++first2;
                         } else if (first1->getColumn() < first2->getColumn()) {
@@ -420,14 +398,14 @@ namespace storm {
                             auto tmp1 = multiplyFactor * first2->getValue();
                             multiplicationTime += std::chrono::high_resolution_clock::now() - multiplicationClock;
                             simplifyClock = std::chrono::high_resolution_clock::now();
-                            tmp1 = simplify(std::move(tmp1));
+                            tmp1 = storm::utility::simplify(tmp1);
                             multiplicationClock = std::chrono::high_resolution_clock::now();
                             simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock;
                             additionClock = std::chrono::high_resolution_clock::now();
                             auto tmp2 = first1->getValue() + tmp1;
                             additionTime += std::chrono::high_resolution_clock::now() - additionClock;
                             simplifyClock = std::chrono::high_resolution_clock::now();
-                            tmp2 = simplify(std::move(tmp2));
+                            tmp2 = storm::utility::simplify(tmp2);
                             simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock;
                             *result = storm::storage::MatrixEntry<typename FlexibleSparseMatrix<ValueType>::index_type, typename FlexibleSparseMatrix<ValueType>::value_type>(first1->getColumn(), tmp2);
                             ++first1;
@@ -440,7 +418,7 @@ namespace storm {
                             auto tmpResult = *first2 * multiplyFactor;
                             multiplicationTime += std::chrono::high_resolution_clock::now() - multiplicationClock;
                             simplifyClock = std::chrono::high_resolution_clock::now();
-                            *result = simplify(std::move(tmpResult));
+                            *result = storm::utility::simplify(tmpResult);
                             simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock;
                         }
                     }
@@ -453,13 +431,13 @@ namespace storm {
                     auto tmp1 = multiplyFactor * oneStepProbabilities[state];
                     multiplicationTime += std::chrono::high_resolution_clock::now() - multiplicationClock;
                     simplifyClock = std::chrono::high_resolution_clock::now();
-                    tmp1 = simplify(std::move(tmp1));
+                    tmp1 = storm::utility::simplify(tmp1);
                     simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock;
                     additionClock2 = std::chrono::high_resolution_clock::now();
                     auto tmp2 = oneStepProbabilities[predecessor] + tmp1;
                     additionTime2 += std::chrono::high_resolution_clock::now() - additionClock2;
                     simplifyClock = std::chrono::high_resolution_clock::now();
-                    tmp2 = simplify(std::move(tmp2));
+                    tmp2 = storm::utility::simplify(tmp2);
                     simplifyTime += std::chrono::high_resolution_clock::now() - simplifyClock;
                     oneStepProbabilities[predecessor] = tmp2;
                     
diff --git a/src/modelchecker/reachability/SparseSccModelChecker.h b/src/modelchecker/reachability/SparseSccModelChecker.h
index a12cf9196..6c530448e 100644
--- a/src/modelchecker/reachability/SparseSccModelChecker.h
+++ b/src/modelchecker/reachability/SparseSccModelChecker.h
@@ -36,9 +36,11 @@ namespace storm {
             class SparseSccModelChecker {
             public:
                 static ValueType computeReachabilityProbability(storm::models::Dtmc<ValueType> const& dtmc, std::shared_ptr<storm::properties::prctl::PrctlFilter<double>> const& filterFormula);
-                
+
+                static ValueType computeReachabilityProbability(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType>& oneStepProbabilities, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, boost::optional<std::vector<std::size_t>> const& distances = {});
+
             private:
-                static uint_fast64_t treatScc(storm::models::Dtmc<ValueType> const& dtmc, FlexibleSparseMatrix<ValueType>& matrix, std::vector<ValueType>& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix<ValueType> const& forwardTransitions, FlexibleSparseMatrix<ValueType>& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level, uint_fast64_t maximalSccSize, std::vector<storm::storage::sparse::state_type>& entryStateQueue, std::vector<std::size_t> const& distances);
+                static uint_fast64_t treatScc(FlexibleSparseMatrix<ValueType>& matrix, std::vector<ValueType>& oneStepProbabilities, storm::storage::BitVector const& entryStates, storm::storage::BitVector const& scc, storm::storage::SparseMatrix<ValueType> const& forwardTransitions, FlexibleSparseMatrix<ValueType>& backwardTransitions, bool eliminateEntryStates, uint_fast64_t level, uint_fast64_t maximalSccSize, std::vector<storm::storage::sparse::state_type>& entryStateQueue, boost::optional<std::vector<std::size_t>> const& distances);
                 static FlexibleSparseMatrix<ValueType> getFlexibleSparseMatrix(storm::storage::SparseMatrix<ValueType> const& matrix, bool setAllValuesToOne = false);
                 static void eliminateState(FlexibleSparseMatrix<ValueType>& matrix, std::vector<ValueType>& oneStepProbabilities, uint_fast64_t state, FlexibleSparseMatrix<ValueType>& backwardTransitions);
                 static bool eliminateStateInPlace(storm::storage::SparseMatrix<ValueType>& matrix, std::vector<ValueType>& oneStepProbabilities, uint_fast64_t state, storm::storage::SparseMatrix<ValueType>& backwardTransitions);
diff --git a/src/properties/prctl/PrctlFilter.h b/src/properties/prctl/PrctlFilter.h
index dc7c2c382..e72d700c8 100644
--- a/src/properties/prctl/PrctlFilter.h
+++ b/src/properties/prctl/PrctlFilter.h
@@ -194,7 +194,6 @@ namespace storm {
 					} else {
 
 						if(this->actions.empty()) {
-
 							// There are no filter actions but only the raw state formula. So just print that.
 							return child->toString();
 						}
diff --git a/src/storage/DeterministicModelStrongBisimulationDecomposition.cpp b/src/storage/DeterministicModelStrongBisimulationDecomposition.cpp
index 27cf123ff..adf78306c 100644
--- a/src/storage/DeterministicModelStrongBisimulationDecomposition.cpp
+++ b/src/storage/DeterministicModelStrongBisimulationDecomposition.cpp
@@ -604,10 +604,9 @@ namespace storm {
             // Now build (a) and (b) by traversing all blocks.
             for (uint_fast64_t blockIndex = 0; blockIndex < this->blocks.size(); ++blockIndex) {
                 auto const& block = this->blocks[blockIndex];
-                
+
                 // Pick one representative state. It doesn't matter which state it is, because they all behave equally.
                 storm::storage::sparse::state_type representativeState = *block.begin();
-                
                 Block const& oldBlock = partition.getBlock(representativeState);
                 
                 // If the block is absorbing, we simply add a self-loop.
@@ -655,7 +654,7 @@ namespace storm {
                 Block const& initialBlock = partition.getBlock(initialState);
                 newLabeling.addAtomicPropositionToState("init", initialBlock.getId());
             }
-            
+                        
             // FIXME:
             // If reward structures are allowed, the quotient structures need to be built here.
             
@@ -701,6 +700,12 @@ namespace storm {
         
         template<typename ValueType>
         void DeterministicModelStrongBisimulationDecomposition<ValueType>::refineBlockProbabilities(Block& block, Partition& partition, std::deque<Block*>& splitterQueue) {
+            // First, we simplify all the values. For most types this does not change anything, but for rational
+            // functions, this achieves a canonicity that is used for sorting the rational functions later.
+            for (auto probabilityIt = partition.getBeginOfValues(block), probabilityIte = partition.getEndOfValues(block); probabilityIt != probabilityIte; ++probabilityIt) {
+                storm::utility::simplify(*probabilityIt);
+            }
+            
             // Sort the states in the block based on their probabilities.
             std::sort(partition.getBeginOfStates(block), partition.getEndOfStates(block), [&partition] (storm::storage::sparse::state_type const& a, storm::storage::sparse::state_type const& b) { return partition.getValue(a) < partition.getValue(b); } );
             
diff --git a/src/stormParametric.cpp b/src/stormParametric.cpp
index f5649a70f..f01fe0936 100644
--- a/src/stormParametric.cpp
+++ b/src/stormParametric.cpp
@@ -1,6 +1,7 @@
 // Include generated headers.
 #include "storm-config.h"
-#include "storm-version.h"
+
+#include <tuple>
 
 // Include other headers.
 #include "src/exceptions/BaseException.h"
@@ -13,11 +14,120 @@
 #include "src/storage/DeterministicModelStrongBisimulationDecomposition.h"
 #include "src/modelchecker/reachability/SparseSccModelChecker.h"
 #include "src/storage/parameters.h"
+#include "src/models/Dtmc.h"
+#include "src/properties/prctl/PrctlFilter.h"
+
+std::tuple<storm::RationalFunction, boost::optional<storm::storage::SparseMatrix<storm::RationalFunction>>, boost::optional<std::vector<storm::RationalFunction>>, boost::optional<storm::storage::BitVector>, boost::optional<double>, boost::optional<bool>> computeReachabilityProbability(storm::models::Dtmc<storm::RationalFunction> const& dtmc, std::shared_ptr<storm::properties::prctl::PrctlFilter<double>> const& filterFormula) {
+    // The first thing we need to do is to make sure the formula is of the correct form and - if so - extract
+    // the bitvector representation of the atomic propositions.
+    
+    std::shared_ptr<storm::properties::prctl::AbstractStateFormula<double>> stateFormula = std::dynamic_pointer_cast<storm::properties::prctl::AbstractStateFormula<double>>(filterFormula->getChild());
+    std::shared_ptr<storm::properties::prctl::AbstractPathFormula<double>> pathFormula;
+    boost::optional<double> threshold;
+    boost::optional<bool> strict;
+    if (stateFormula != nullptr) {
+        std::shared_ptr<storm::properties::prctl::ProbabilisticBoundOperator<double>> probabilisticBoundFormula = std::dynamic_pointer_cast<storm::properties::prctl::ProbabilisticBoundOperator<double>>(stateFormula);
+        STORM_LOG_THROW(probabilisticBoundFormula != nullptr, storm::exceptions::InvalidPropertyException, "Illegal formula " << *filterFormula << " for parametric model checking. Note that only unbounded reachability properties are permitted.");
+        STORM_LOG_THROW(probabilisticBoundFormula->getComparisonOperator() == storm::properties::ComparisonType::LESS_EQUAL || probabilisticBoundFormula->getComparisonOperator() == storm::properties::ComparisonType::LESS, storm::exceptions::InvalidPropertyException, "Illegal formula " << *filterFormula << " for parametric model checking. Note that only unbounded reachability properties with upper probability bounds are permitted.");
+        
+        threshold = probabilisticBoundFormula->getBound();
+        strict = probabilisticBoundFormula->getComparisonOperator() == storm::properties::ComparisonType::LESS;
+        pathFormula = probabilisticBoundFormula->getChild();
+    } else {
+        pathFormula = std::dynamic_pointer_cast<storm::properties::prctl::AbstractPathFormula<double>>(filterFormula->getChild());
+    }
+    
+    STORM_LOG_THROW(pathFormula != nullptr, storm::exceptions::InvalidPropertyException, "Illegal formula " << *filterFormula << " for parametric model checking. Note that only unbounded reachability properties are permitted.");
+    
+    std::shared_ptr<storm::properties::prctl::Until<double>> untilFormula = std::dynamic_pointer_cast<storm::properties::prctl::Until<double>>(pathFormula);
+    std::shared_ptr<storm::properties::prctl::AbstractStateFormula<double>> phiStateFormula;
+    std::shared_ptr<storm::properties::prctl::AbstractStateFormula<double>> psiStateFormula;
+    if (untilFormula != nullptr) {
+        phiStateFormula = untilFormula->getLeft();
+        psiStateFormula = untilFormula->getRight();
+    } else {
+        std::shared_ptr<storm::properties::prctl::Eventually<double>> eventuallyFormula = std::dynamic_pointer_cast<storm::properties::prctl::Eventually<double>>(pathFormula);
+        STORM_LOG_THROW(eventuallyFormula != nullptr, storm::exceptions::InvalidPropertyException, "Illegal formula " << *filterFormula << " for parametric model checking. Note that only unbounded reachability properties are permitted.");
+        phiStateFormula = std::shared_ptr<storm::properties::prctl::Ap<double>>(new storm::properties::prctl::Ap<double>("true"));
+        psiStateFormula = eventuallyFormula->getChild();
+    }
+    
+    // Now we need to make sure the formulas defining the phi and psi states are just labels.
+    std::shared_ptr<storm::properties::prctl::Ap<double>> phiStateFormulaApFormula = std::dynamic_pointer_cast<storm::properties::prctl::Ap<double>>(phiStateFormula);
+    std::shared_ptr<storm::properties::prctl::Ap<double>> psiStateFormulaApFormula = std::dynamic_pointer_cast<storm::properties::prctl::Ap<double>>(psiStateFormula);
+    STORM_LOG_THROW(phiStateFormulaApFormula.get() != nullptr, storm::exceptions::InvalidPropertyException, "Illegal formula " << *phiStateFormula << " for parametric model checking. Note that only atomic propositions are admitted in that position.");
+    STORM_LOG_THROW(psiStateFormulaApFormula.get() != nullptr, storm::exceptions::InvalidPropertyException, "Illegal formula " << *psiStateFormula << " for parametric model checking. Note that only atomic propositions are admitted in that position.");
+    
+    // Now retrieve the appropriate bitvectors from the atomic propositions.
+    storm::storage::BitVector phiStates = phiStateFormulaApFormula->getAp() != "true" ? dtmc.getLabeledStates(phiStateFormulaApFormula->getAp()) : storm::storage::BitVector(dtmc.getNumberOfStates(), true);
+    storm::storage::BitVector psiStates = dtmc.getLabeledStates(psiStateFormulaApFormula->getAp());
+    
+    // Do some sanity checks to establish some required properties.
+    STORM_LOG_THROW(dtmc.getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::IllegalArgumentException, "Input model is required to have exactly one initial state.");
+    
+    // Then, compute the subset of states that has a probability of 0 or 1, respectively.
+    std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(dtmc, phiStates, psiStates);
+    storm::storage::BitVector statesWithProbability0 = statesWithProbability01.first;
+    storm::storage::BitVector statesWithProbability1 = statesWithProbability01.second;
+    storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
+    
+    // If the initial state is known to have either probability 0 or 1, we can directly return the result.
+    if (dtmc.getInitialStates().isDisjointFrom(maybeStates)) {
+        STORM_LOG_DEBUG("The probability of all initial states was found in a preprocessing step.");
+        return statesWithProbability0.get(*dtmc.getInitialStates().begin()) ? storm::utility::constantZero<storm::RationalFunction>() : storm::utility::constantOne<storm::RationalFunction>();
+    }
+    
+    // Determine the set of states that is reachable from the initial state without jumping over a target state.
+    storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(dtmc.getTransitionMatrix(), dtmc.getInitialStates(), maybeStates, statesWithProbability1);
+    
+    // Subtract from the maybe states the set of states that is not reachable (on a path from the initial to a target state).
+    maybeStates &= reachableStates;
+    
+    // Create a vector for the probabilities to go to a state with probability 1 in one step.
+    std::vector<storm::RationalFunction> oneStepProbabilities = dtmc.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1);
+    
+    // Determine the set of initial states of the sub-DTMC.
+    storm::storage::BitVector newInitialStates = dtmc.getInitialStates() % maybeStates;
+    
+    // We then build the submatrix that only has the transitions of the maybe states.
+    storm::storage::SparseMatrix<storm::RationalFunction> submatrix = dtmc.getTransitionMatrix().getSubmatrix(false, maybeStates, maybeStates);
+    
+    // To be able to apply heuristics later, we now determine the distance of each state to the initial state.
+    std::vector<std::pair<storm::storage::sparse::state_type, std::size_t>> stateQueue;
+    stateQueue.reserve(submatrix.getRowCount());
+    storm::storage::BitVector statesInQueue(submatrix.getRowCount());
+    std::vector<std::size_t> distances(submatrix.getRowCount());
+    
+    storm::storage::sparse::state_type currentPosition = 0;
+    for (auto const& initialState : newInitialStates) {
+        stateQueue.emplace_back(initialState, 0);
+        statesInQueue.set(initialState);
+    }
+    
+    // Perform a BFS.
+    while (currentPosition < stateQueue.size()) {
+        std::pair<storm::storage::sparse::state_type, std::size_t> const& stateDistancePair = stateQueue[currentPosition];
+        distances[stateDistancePair.first] = stateDistancePair.second;
+        
+        for (auto const& successorEntry : submatrix.getRow(stateDistancePair.first)) {
+            if (!statesInQueue.get(successorEntry.getColumn())) {
+                stateQueue.emplace_back(successorEntry.getColumn(), stateDistancePair.second + 1);
+                statesInQueue.set(successorEntry.getColumn());
+            }
+        }
+        ++currentPosition;
+    }
+    
+    storm::modelchecker::reachability::SparseSccModelChecker<storm::RationalFunction> modelchecker;
+    
+    return std::make_tuple(modelchecker.computeReachabilityProbability(submatrix, oneStepProbabilities, submatrix.transpose(), newInitialStates, phiStates, psiStates, distances),submatrix, oneStepProbabilities, newInitialStates, threshold, strict);
+}
+
 /*!
  * Main entry point of the executable storm.
  */
 int main(const int argc, const char** argv) {
-    try {
+//    try {
         storm::utility::cli::setUp();
         storm::utility::cli::printHeader(argc, argv);
         bool optionsCorrect = storm::utility::cli::parseOptions(argc, argv);
@@ -28,7 +138,7 @@ int main(const int argc, const char** argv) {
         // From this point on we are ready to carry out the actual computations.
         // Program Translation Time Measurement, Start
         std::chrono::high_resolution_clock::time_point programTranslationStart = std::chrono::high_resolution_clock::now();
-
+        
         // First, we build the model using the given symbolic model description and constant definitions.
         std::string const& programFile = storm::settings::generalSettings().getSymbolicModelFilename();
         std::string const& constants = storm::settings::generalSettings().getConstantDefinitionString();
@@ -36,31 +146,31 @@ int main(const int argc, const char** argv) {
         std::shared_ptr<storm::models::AbstractModel<storm::RationalFunction>> model = storm::adapters::ExplicitModelAdapter<storm::RationalFunction>::translateProgram(program, constants);
         
         model->printModelInformationToStream(std::cout);
-
+        
         // Program Translation Time Measurement, End
         std::chrono::high_resolution_clock::time_point programTranslationEnd = std::chrono::high_resolution_clock::now();
         std::cout << "Parsing and translating the model took " << std::chrono::duration_cast<std::chrono::milliseconds>(programTranslationEnd - programTranslationStart).count() << "ms." << std::endl << std::endl;
-
+        
         std::shared_ptr<storm::models::Dtmc<storm::RationalFunction>> dtmc = model->as<storm::models::Dtmc<storm::RationalFunction>>();
         
         // Perform bisimulation minimization if requested.
         if (storm::settings::generalSettings().isBisimulationSet()) {
             storm::storage::DeterministicModelStrongBisimulationDecomposition<storm::RationalFunction> bisimulationDecomposition(*dtmc, true);
             dtmc = bisimulationDecomposition.getQuotient()->as<storm::models::Dtmc<storm::RationalFunction>>();
+            
+            dtmc->printModelInformationToStream(std::cout);
         }
         
         assert(dtmc);
         storm::modelchecker::reachability::CollectConstraints<storm::RationalFunction> constraintCollector;
         constraintCollector(*dtmc);
-
-        
-        storm::modelchecker::reachability::SparseSccModelChecker<storm::RationalFunction> modelChecker;
         
         STORM_LOG_THROW(storm::settings::generalSettings().isPctlPropertySet(), storm::exceptions::InvalidSettingsException, "Unable to perform model checking without a property.");
         std::shared_ptr<storm::properties::prctl::PrctlFilter<double>> filterFormula = storm::parser::PrctlParser::parsePrctlFormula(storm::settings::generalSettings().getPctlProperty());
         
-        storm::RationalFunction value = modelChecker.computeReachabilityProbability(*dtmc, filterFormula);
-        STORM_PRINT_AND_LOG(std::endl << "computed value " << value << std::endl);
+        std::tuple<storm::RationalFunction, boost::optional<storm::storage::SparseMatrix<storm::RationalFunction>>, boost::optional<std::vector<storm::RationalFunction>>, boost::optional<storm::storage::BitVector>, boost::optional<double>, boost::optional<bool>> result = computeReachabilityProbability(*dtmc, filterFormula);
+        storm::RationalFunction valueFunction = std::get<0>(result);
+        STORM_PRINT_AND_LOG(std::endl << "computed value " << valueFunction << std::endl);
         
         // Get variables from parameter definitions in prism program.
         std::set<storm::Variable> parameters;
@@ -73,26 +183,26 @@ int main(const int argc, const char** argv) {
                 parameters.insert(p);
             }
         }
-        // 
-        STORM_LOG_ASSERT(parameters == value.gatherVariables(), "Parameters in result and program definition do not coincide.");
+        //
+        STORM_LOG_ASSERT(parameters == valueFunction.gatherVariables(), "Parameters in result and program definition do not coincide.");
         
         if(storm::settings::parametricSettings().exportResultToFile()) {
-            storm::utility::exportParametricMcResult(value, constraintCollector);
+            storm::utility::exportParametricMcResult(valueFunction, constraintCollector);
+        }
+        
+        if (storm::settings::parametricSettings().exportToSmt2File() && std::get<1>(result) && std::get<2>(result) && std::get<3>(result) && std::get<4>(result) && std::get<5>(result)) {
+            storm::modelchecker::reachability::DirectEncoding dec;
+            storm::utility::exportStringStreamToFile(dec.encodeAsSmt2(std::get<1>(result).get(), std::get<2>(result).get(), parameters, std::get<3>(result).get(), carl::rationalize<storm::RationalFunction::CoeffType>(std::get<4>(result).get()), std::get<5>(result).get()), "out.smt");
         }
         
-//        if (storm::settings::parametricSettings().exportToSmt2File()) {
-//            storm::modelchecker::reachability::DirectEncoding dec;
-//            storm::utility::exportStringStreamToFile(dec.encodeAsSmt2(modelChecker.getMatrix(), parameters,));
-//        }        
-
         // All operations have now been performed, so we clean up everything and terminate.
         storm::utility::cli::cleanUp();
         return 0;
-    } catch (storm::exceptions::BaseException const& exception) {
-        STORM_LOG_ERROR("An exception caused StoRM to terminate. The message of the exception is: " << exception.what());
-    } catch (std::exception const& exception) {
-        STORM_LOG_ERROR("An unexpected exception occurred and caused StoRM to terminate. The message of this exception is: " << exception.what());
-    }
+//    } catch (storm::exceptions::BaseException const& exception) {
+//        STORM_LOG_ERROR("An exception caused StoRM to terminate. The message of the exception is: " << exception.what());
+//    } catch (std::exception const& exception) {
+//        STORM_LOG_ERROR("An unexpected exception occurred and caused StoRM to terminate. The message of this exception is: " << exception.what());
+//    }
 }
 
 //#include <memory>
@@ -119,33 +229,33 @@ int main(const int argc, const char** argv) {
 //
 //std::string ParametricStormEntryPoint::reachabilityToSmt2(std::string const& label)
 //{
-//    
+//
 //    storm::storage::BitVector phiStates(mModel->getNumberOfStates(), true);
 //    storm::storage::BitVector initStates = mModel->getInitialStates();
 //    storm::storage::BitVector targetStates = mModel->getLabeledStates(label);
-//    
+//
 //    std::shared_ptr<models::Dtmc<RationalFunction>> dtmc = mModel->as<models::Dtmc<RationalFunction>>();
 //    // 1. make target states absorbing.
 //    dtmc->makeAbsorbing(targetStates);
 //    // 2. throw away anything which does not add to the reachability probability.
 //    // 2a. remove non productive states
 //    storm::storage::BitVector productiveStates = utility::graph::performProbGreater0(*dtmc, dtmc->getBackwardTransitions(), phiStates, targetStates);
-//    // 2b. calculate set of states wich 
+//    // 2b. calculate set of states wich
 //    storm::storage::BitVector almostSurelyReachingTargetStates = ~utility::graph::performProbGreater0(*dtmc, dtmc->getBackwardTransitions(), phiStates, ~productiveStates);
 //    // 2c. Make such states also target states.
 //    dtmc->makeAbsorbing(almostSurelyReachingTargetStates);
-//    // 2d. throw away non reachable states 
+//    // 2d. throw away non reachable states
 //    storm::storage::BitVector reachableStates = utility::graph::performProbGreater0(*dtmc, dtmc->getTransitionMatrix(), phiStates, initStates);
 //    storm::storage::BitVector bv = productiveStates & reachableStates;
 //    dtmc->getStateLabeling().addAtomicProposition("__targets__", targetStates | almostSurelyReachingTargetStates);
 //    models::Dtmc<RationalFunction> subdtmc = dtmc->getSubDtmc(bv);
-//    
+//
 //    phiStates = storm::storage::BitVector(subdtmc.getNumberOfStates(), true);
 //    initStates = subdtmc.getInitialStates();
 //    targetStates = subdtmc.getLabeledStates("__targets__");
 //    storm::storage::BitVector deadlockStates(phiStates);
 //    deadlockStates.set(subdtmc.getNumberOfStates()-1,false);
-//    
+//
 //    // Search for states with only one non-deadlock successor.
 //    std::map<StateId, storage::DeterministicTransition<RationalFunction>> chainedStates;
 //    StateId nrStates = subdtmc.getNumberOfStates();
@@ -197,8 +307,8 @@ int main(const int argc, const char** argv) {
 //            eliminatedStates.set(it->first, true);
 //        }
 //    }
-//    
-//    
+//
+//
 //    for(auto chainedState : chainedStates)
 //    {
 //        if(!eliminatedStates[chainedState.first])
@@ -206,15 +316,15 @@ int main(const int argc, const char** argv) {
 //            std::cout << chainedState.first << " -- " << chainedState.second.probability() << " --> " << chainedState.second.targetState() << std::endl;
 //        }
 //    }
-//    
+//
 //    storage::StronglyConnectedComponentDecomposition<RationalFunction> sccs(subdtmc);
 //    std::cout << sccs << std::endl;
 //
 //    modelchecker::reachability::DirectEncoding dec;
 //    std::vector<carl::Variable> parameters;
-    
+
 //    return dec.encodeAsSmt2(subdtmc, parameters, subdtmc.getLabeledStates("init"), subdtmc.getLabeledStates("__targets__"), mpq_class(1,2));
-//    
+//
 //}
 //
 //
@@ -229,8 +339,8 @@ int main(const int argc, const char** argv) {
 //        fstream << entry.reachabilityToSmt2(s->getOptionByLongName("reachability").getArgument(0).getValueAsString());
 //        fstream.close();
 //    }
-//    
-//    
+//
+//
 //}
 //
 //}
diff --git a/src/utility/ConstantsComparator.cpp b/src/utility/ConstantsComparator.cpp
new file mode 100644
index 000000000..24b86231e
--- /dev/null
+++ b/src/utility/ConstantsComparator.cpp
@@ -0,0 +1,136 @@
+#include "src/utility/ConstantsComparator.h"
+
+#include "src/storage/sparse/StateType.h"
+
+namespace storm {
+    namespace utility {
+        
+        template<typename ValueType>
+        ValueType one() {
+            return ValueType(1);
+        }
+        
+        template<typename ValueType>
+        ValueType zero() {
+            return ValueType(0);
+        }
+        
+        template<typename ValueType>
+        ValueType infinity() {
+            return std::numeric_limits<ValueType>::infinity();
+        }
+        
+        template<typename ValueType>
+        ValueType& simplify(ValueType& value) {
+            // In the general case, we don't to anything here, but merely return the value. If something else is
+            // supposed to happen here, the templated function can be specialized for this particular type.
+            return value;
+        }
+        
+        template<typename ValueType>
+        bool ConstantsComparator<ValueType>::isOne(ValueType const& value) const {
+            return value == one<ValueType>();
+        }
+        
+        template<typename ValueType>
+        bool ConstantsComparator<ValueType>::isZero(ValueType const& value) const {
+            return value == zero<ValueType>();
+        }
+        
+        template<typename ValueType>
+        bool ConstantsComparator<ValueType>::isEqual(ValueType const& value1, ValueType const& value2) const {
+            return value1 == value2;
+        }
+        
+        ConstantsComparator<double>::ConstantsComparator() : precision(storm::settings::generalSettings().getPrecision()) {
+            // Intentionally left empty.
+        }
+        
+        ConstantsComparator<double>::ConstantsComparator(double precision) : precision(precision) {
+            // Intentionally left empty.
+        }
+        
+        bool ConstantsComparator<double>::isOne(double const& value) const {
+            return std::abs(value - one<double>()) <= precision;
+        }
+        
+        bool ConstantsComparator<double>::isZero(double const& value) const {
+            return std::abs(value) <= precision;
+        }
+        
+        bool ConstantsComparator<double>::isEqual(double const& value1, double const& value2) const {
+            return std::abs(value1 - value2) <= precision;
+        }
+        
+        bool ConstantsComparator<double>::isConstant(double const& value) const {
+            return true;
+        }
+        
+#ifdef PARAMETRIC_SYSTEMS
+        template<>
+        RationalFunction& simplify(RationalFunction& value) {
+            value.simplify();
+            return value;
+        }
+        
+        bool ConstantsComparator<storm::RationalFunction>::isOne(storm::RationalFunction const& value) const {
+            return value.isOne();
+        }
+        
+        bool ConstantsComparator<storm::RationalFunction>::isZero(storm::RationalFunction const& value) const {
+            return value.isZero();
+        }
+        
+        bool ConstantsComparator<storm::RationalFunction>::isEqual(storm::RationalFunction const& value1, storm::RationalFunction const& value2) const {
+            return value1 == value2;
+        }
+        
+        bool ConstantsComparator<storm::RationalFunction>::isConstant(storm::RationalFunction const& value) const {
+            return value.isConstant();
+        }
+        
+        bool ConstantsComparator<storm::Polynomial>::isOne(storm::Polynomial const& value) const {
+            return value.isOne();
+        }
+        
+        bool ConstantsComparator<storm::Polynomial>::isZero(storm::Polynomial const& value) const {
+            return value.isZero();
+        }
+        
+        bool ConstantsComparator<storm::Polynomial>::isEqual(storm::Polynomial const& value1, storm::Polynomial const& value2) const {
+            return value1 == value2;
+        }
+        
+        bool ConstantsComparator<storm::Polynomial>::isConstant(storm::Polynomial const& value) const {
+            return value.isConstant();
+        }
+        
+#endif
+        
+        template<typename IndexType, typename ValueType>
+        storm::storage::MatrixEntry<IndexType, ValueType>& simplify(storm::storage::MatrixEntry<IndexType, ValueType>& matrixEntry) {
+            simplify(matrixEntry.getValue());
+            return matrixEntry;
+        }
+        
+        template class ConstantsComparator<double>;
+        template class ConstantsComparator<RationalFunction>;
+        template class ConstantsComparator<Polynomial>;
+        
+        template double one();
+        template double zero();
+        template double infinity();
+
+        template RationalFunction one();
+        template RationalFunction zero();
+
+        template Polynomial one();
+        template Polynomial zero();
+
+        template double& simplify(double& value);
+        template RationalFunction& simplify(RationalFunction& value);
+        
+        template storm::storage::MatrixEntry<storm::storage::sparse::state_type, double>& simplify(storm::storage::MatrixEntry<storm::storage::sparse::state_type, double>& matrixEntry);
+        template storm::storage::MatrixEntry<storm::storage::sparse::state_type, RationalFunction>& simplify(storm::storage::MatrixEntry<storm::storage::sparse::state_type, RationalFunction>& matrixEntry);
+    }
+}
\ No newline at end of file
diff --git a/src/utility/ConstantsComparator.h b/src/utility/ConstantsComparator.h
index d97c4ef93..431bb974f 100644
--- a/src/utility/ConstantsComparator.h
+++ b/src/utility/ConstantsComparator.h
@@ -13,69 +13,49 @@
 #include <cstdint>
 
 #include "src/settings/SettingsManager.h"
+#include "src/storage/SparseMatrix.h"
 
 namespace storm {
     namespace utility {
         
         template<typename ValueType>
-        ValueType one() {
-            return ValueType(1);
-        }
+        ValueType one();
         
         template<typename ValueType>
-        ValueType zero() {
-            return ValueType(0);
-        }
+        ValueType zero();
         
         template<typename ValueType>
-        ValueType infinity() {
-            return std::numeric_limits<ValueType>::infinity();
-        }
+        ValueType infinity();
+        
+        template<typename ValueType>
+        ValueType& simplify(ValueType& value);
         
         // A class that can be used for comparing constants.
         template<typename ValueType>
         class ConstantsComparator {
         public:
-            bool isOne(ValueType const& value) const {
-                return value == one<ValueType>();
-            }
+            bool isOne(ValueType const& value) const;
             
-            bool isZero(ValueType const& value) const {
-                return value == zero<ValueType>();
-            }
+            bool isZero(ValueType const& value) const;
             
-            bool isEqual(ValueType const& value1, ValueType const& value2) const {
-                return value1 == value2;
-            }
+            bool isEqual(ValueType const& value1, ValueType const& value2) const;
         };
         
         // For doubles we specialize this class and consider the comparison modulo some predefined precision.
         template<>
         class ConstantsComparator<double> {
         public:
-            ConstantsComparator() : precision(storm::settings::generalSettings().getPrecision()) {
-                // Intentionally left empty.
-            }
+            ConstantsComparator();
 
-            ConstantsComparator(double precision) : precision(precision) {
-                // Intentionally left empty.
-            }
+            ConstantsComparator(double precision);
             
-            bool isOne(double const& value) const {
-                return std::abs(value - one<double>()) <= precision;
-            }
+            bool isOne(double const& value) const;
             
-            bool isZero(double const& value) const {
-                return std::abs(value) <= precision;
-            }
+            bool isZero(double const& value) const;
             
-            bool isEqual(double const& value1, double const& value2) const {
-                return std::abs(value1 - value2) <= precision;
-            }
+            bool isEqual(double const& value1, double const& value2) const;
             
-            bool isConstant(double const& value) const {
-                return true;
-            }
+            bool isConstant(double const& value) const;
             
         private:
             // The precision used for comparisons.
@@ -83,47 +63,38 @@ namespace storm {
         };
         
 #ifdef PARAMETRIC_SYSTEMS
+        template<>
+        RationalFunction& simplify(RationalFunction& value);
+        
         template<>
         class ConstantsComparator<storm::RationalFunction> {
         public:
-            bool isOne(storm::RationalFunction const& value) const {
-                return value.isOne();
-            }
+            bool isOne(storm::RationalFunction const& value) const;
             
-            bool isZero(storm::RationalFunction const& value) const {
-                return value.isZero();
-            }
+            bool isZero(storm::RationalFunction const& value) const;
             
-            bool isEqual(storm::RationalFunction const& value1, storm::RationalFunction const& value2) const {
-                return value1 == value2;
-            }
+            bool isEqual(storm::RationalFunction const& value1, storm::RationalFunction const& value2) const;
             
-            bool isConstant(storm::RationalFunction const& value) const {
-                return value.isConstant();
-            }
+            bool isConstant(storm::RationalFunction const& value) const;
         };
         
         template<>
         class ConstantsComparator<storm::Polynomial> {
         public:
-            bool isOne(storm::Polynomial const& value) const {
-                return value.isOne();
-            }
+            bool isOne(storm::Polynomial const& value) const;
             
-            bool isZero(storm::Polynomial const& value) const {
-                return value.isZero();
-            }
+            bool isZero(storm::Polynomial const& value) const;
             
-            bool isEqual(storm::Polynomial const& value1, storm::Polynomial const& value2) const {
-                return value1 == value2;
-            }
+            bool isEqual(storm::Polynomial const& value1, storm::Polynomial const& value2) const;
             
-            bool isConstant(storm::Polynomial const& value) const {
-                return value.isConstant();
-            }
+            bool isConstant(storm::Polynomial const& value) const;
         };
+        
+        
 #endif
         
+        template<typename IndexType, typename ValueType>
+        storm::storage::MatrixEntry<IndexType, ValueType>& simplify(storm::storage::MatrixEntry<IndexType, ValueType>& matrixEntry);
     }
 }
 
diff --git a/src/utility/export.h b/src/utility/export.h
index 7d66bd0b3..943da4608 100644
--- a/src/utility/export.h
+++ b/src/utility/export.h
@@ -35,11 +35,11 @@ void exportParametricMcResult(const ValueType& mcresult, storm::modelchecker::re
 	filestream.close();
 }
 
-void exportStringStreamToFile(std::stringstream& stream, std::string filepath) {
+void exportStringStreamToFile(std::string const& str, std::string filepath) {
 	// todo add checks.
 	std::ofstream filestream;
 	filestream.open(filepath);
-	filestream << stream.str();
+	filestream << str;
 	filestream.close();
 }