diff --git a/examples/dtmc/crowds/crowds15_5.pm b/examples/dtmc/crowds/crowds15_5.pm
index 511d962f0..95ab8a910 100644
--- a/examples/dtmc/crowds/crowds15_5.pm
+++ b/examples/dtmc/crowds/crowds15_5.pm
@@ -2,9 +2,9 @@ dtmc
 
 // probability of forwarding
 const double    PF = 0.8;
-const double notPF = .2;  // must be 1-PF
+const double notPF = 0.2;  // must be 1-PF
 // probability that a crowd member is bad
-const double  badC = .167;
+const double  badC = 0.167;
  // probability that a crowd member is good
 const double goodC = 0.833;
 // Total number of protocol runs to analyze
@@ -92,4 +92,4 @@ endmodule
 
 label "observe0Greater1" = observe0 > 1;
 label "observeIGreater1" = observe1 > 1 | observe2 > 1 | observe3 > 1 | observe4 > 1 | observe5 > 1 | observe6 > 1 | observe7 > 1 | observe8 > 1 | observe9 > 1 | observe10 > 1 | observe11 > 1 | observe12 > 1 | observe13 > 1 | observe14 > 1;
-label "observeOnlyTrueSender" = observe0 > 1 & observe1 <= 1 & observe2 <= 1 & observe3 <= 1 & observe4 <= 1 & observe5 <= 1 & observe6 <= 1 & observe7 <= 1 & observe8 <= 1 & observe9 <= 1 & observe10 <= 1 & observe11 <= 1 & observe12 <= 1 & observe13 <= 1 & observe14 <= 1;
\ No newline at end of file
+label "observeOnlyTrueSender" = observe0 > 1 & observe1 <= 1 & observe2 <= 1 & observe3 <= 1 & observe4 <= 1 & observe5 <= 1 & observe6 <= 1 & observe7 <= 1 & observe8 <= 1 & observe9 <= 1 & observe10 <= 1 & observe11 <= 1 & observe12 <= 1 & observe13 <= 1 & observe14 <= 1;
diff --git a/examples/dtmc/crowds/crowds20_5.pm b/examples/dtmc/crowds/crowds20_5.pm
index 31c63770a..1809c22d7 100644
--- a/examples/dtmc/crowds/crowds20_5.pm
+++ b/examples/dtmc/crowds/crowds20_5.pm
@@ -2,9 +2,9 @@ dtmc
 
 // probability of forwarding
 const double    PF = 0.8;
-const double notPF = .2;  // must be 1-PF
+const double notPF = 0.2;  // must be 1-PF
 // probability that a crowd member is bad
-const double  badC = .167;
+const double  badC = 0.167;
  // probability that a crowd member is good
 const double goodC = 0.833;
 // Total number of protocol runs to analyze
@@ -107,4 +107,5 @@ endmodule
 
 label "observe0Greater1" = observe0 > 1;
 label "observeIGreater1" = observe1 > 1 | observe2 > 1 | observe3 > 1 | observe4 > 1 | observe5 > 1 | observe6 > 1 | observe7 > 1 | observe8 > 1 | observe9 > 1 | observe10 > 1 | observe11 > 1 | observe12 > 1 | observe13 > 1 | observe14 > 1 | observe15 > 1 | observe16 > 1 | observe17 > 1 | observe18 > 1 | observe19 > 1;
-label "observeOnlyTrueSender" = observe0 > 1 & observe1 <= 1 & observe2 <= 1 & observe3 <= 1 & observe4 <= 1 & observe5 <= 1 & observe6 <= 1 & observe7 <= 1 & observe8 <= 1 & observe9 <= 1 & observe10 <= 1 & observe11 <= 1 & observe12 <= 1 & observe13 <= 1 & observe14 <= 1 & observe15 <= 1 & observe16 <= 1 & observe17 <= 1 & observe18 <= 1 & observe19 <= 1;
\ No newline at end of file
+label "observeOnlyTrueSender" = observe0 > 1 & observe1 <= 1 & observe2 <= 1 & observe3 <= 1 & observe4 <= 1 & observe5 <= 1 & observe6 <= 1 & observe7 <= 1 & observe8 <= 1 & observe9 <= 1 & observe10 <= 1 & observe11 <= 1 & observe12 <= 1 & observe13 <= 1 & observe14 <= 1 & observe15 <= 1 & observe16 <= 1 & observe17 <= 1 & observe18 <= 1 & observe19 <= 1;
+
diff --git a/resources/3rdparty/cudd-2.5.0/src/cudd/cuddAddApply.c b/resources/3rdparty/cudd-2.5.0/src/cudd/cuddAddApply.c
index fc1ab8c15..e789bfc73 100644
--- a/resources/3rdparty/cudd-2.5.0/src/cudd/cuddAddApply.c
+++ b/resources/3rdparty/cudd-2.5.0/src/cudd/cuddAddApply.c
@@ -1078,7 +1078,7 @@ Cudd_addMod(
     F = *f; G = *g;
     if (cuddIsConstant(F) && cuddIsConstant(G)) {
         // If g is <=0, then result is NaN
-        if (cuddV(G) <= 0) value = (0.0/0.0);
+        if (cuddV(G) <= 0) value = (NAN);
         // Take care of negative case (% is remainder, not modulo)
         else {
             rem = ((int)cuddV(F) % (int)cuddV(G));
@@ -1119,9 +1119,9 @@ Cudd_addLogXY(
     F = *f; G = *g;
     if (cuddIsConstant(F) && cuddIsConstant(G)) {
         // If base is <=0 or ==1 (or +Inf/NaN), then result is NaN
-        if (cuddV(G) <= 0 || cuddV(G) == 1.0 || G==DD_PLUS_INFINITY(dd) || cuddV(G) != cuddV(G)) value = (0.0/0.0);
+        if (cuddV(G) <= 0 || cuddV(G) == 1.0 || G==DD_PLUS_INFINITY(dd) || cuddV(G) != cuddV(G)) value = (NAN);
         // If arg is <0 or NaN, then result is NaN
-        else if (cuddV(F) < 0 || cuddV(F) != cuddV(F)) value = (0.0/0.0);
+        else if (cuddV(F) < 0 || cuddV(F) != cuddV(F)) value = (NAN);
         // If arg is +Inf, then result is +Inf
         else if (F==DD_PLUS_INFINITY(dd)) return DD_PLUS_INFINITY(dd);
         // If arg is (positive/negative) 0, then result is -Inf
diff --git a/src/modelchecker/AbstractModelChecker.cpp b/src/modelchecker/AbstractModelChecker.cpp
index 7ba1de099..2b6d2ef3a 100644
--- a/src/modelchecker/AbstractModelChecker.cpp
+++ b/src/modelchecker/AbstractModelChecker.cpp
@@ -216,8 +216,8 @@ namespace storm {
             }
         }
         
-        std::unique_ptr<CheckResult> AbstractModelChecker::checkExpectedTimeOperatorFormula(storm::logic::ExpectedTimeOperatorFormula const& stateFormula) {
-            STORM_LOG_THROW(stateFormula.getSubformula().isStateFormula(), storm::exceptions::InvalidArgumentException, "The given formula is invalid.");
+		std::unique_ptr<CheckResult> AbstractModelChecker::checkExpectedTimeOperatorFormula(storm::logic::ExpectedTimeOperatorFormula const& stateFormula) {
+			STORM_LOG_THROW(stateFormula.getSubformula().isEventuallyFormula(), storm::exceptions::InvalidArgumentException, "The given formula is invalid.");
             
             // If the reward bound is 0, is suffices to do qualitative model checking.
             bool qualitative = false;
@@ -248,8 +248,8 @@ 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> AbstractModelChecker::checkLongRunAverageOperatorFormula(storm::logic::LongRunAverageOperatorFormula const& stateFormula) {
+			STORM_LOG_THROW(stateFormula.getSubformula().isStateFormula(), storm::exceptions::InvalidArgumentException, "The given formula is invalid.");
             
             std::unique_ptr<CheckResult> result;
             if (stateFormula.hasOptimalityType()) {
diff --git a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
index abec6a929..e3149901f 100644
--- a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
+++ b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
@@ -354,14 +354,7 @@ namespace storm {
                             // Now insert all (cumulative) probability values that target an MEC.
                             for (uint_fast64_t targetMecIndex = 0; targetMecIndex < auxiliaryStateToProbabilityMap.size(); ++targetMecIndex) {
                                 if (auxiliaryStateToProbabilityMap[targetMecIndex] != 0) {
-                                    // If the target MEC is the same as the current one, instead of adding a transition, we need to add the weighted reward
-                                    // to the right-hand side vector of the SSP.
-                                    if (mecIndex == targetMecIndex) {
-                                        b.back() += auxiliaryStateToProbabilityMap[mecIndex] * lraValuesForEndComponents[mecIndex];
-                                    } else {
-                                        // Otherwise, we add a transition to the auxiliary state that is associated with the target MEC.
-                                        sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]);
-                                    }
+                                    sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]);
                                 }
                             }
                             
@@ -390,7 +383,7 @@ namespace storm {
             
             // Set the values for all states in MECs.
             for (auto state : statesInMecs) {
-                result[state] = lraValuesForEndComponents[stateToMecIndexMap[state]];
+				result[state] = x[firstAuxiliaryStateIndex + stateToMecIndexMap[state]];
             }
             
             return result;
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
index a7775139a..8f63110c8 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
@@ -1,6 +1,7 @@
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
 
 #include <vector>
+#include <memory>
 
 #include "src/utility/macros.h"
 #include "src/utility/vector.h"
@@ -11,6 +12,8 @@
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 
 #include "src/exceptions/InvalidStateException.h"
+#include "src/storage/StronglyConnectedComponentDecomposition.h"
+
 #include "src/exceptions/InvalidPropertyException.h"
 
 namespace storm {
@@ -300,6 +303,233 @@ namespace storm {
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative)));
         }
         
+		template<typename ValueType>
+		std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(storm::storage::BitVector const& psiStates, bool qualitative) const {
+			// If there are no goal states, we avoid the computation and directly return zero.
+			auto numOfStates = this->getModel().getNumberOfStates();
+			if (psiStates.empty()) {
+				return std::vector<ValueType>(numOfStates, storm::utility::zero<ValueType>());
+			}
+
+			// Likewise, if all bits are set, we can avoid the computation and set.
+			if ((~psiStates).empty()) {
+				return std::vector<ValueType>(numOfStates, storm::utility::one<ValueType>());
+			}
+
+			// Start by decomposing the DTMC into its BSCCs.
+			storm::storage::StronglyConnectedComponentDecomposition<double> bsccDecomposition(this->getModel(), false, true);
+
+			// Get some data members for convenience.
+			typename storm::storage::SparseMatrix<ValueType> const& transitionMatrix = this->getModel().getTransitionMatrix();
+			ValueType one = storm::utility::one<ValueType>();
+			ValueType zero = storm::utility::zero<ValueType>();
+
+			// First we check which states are in BSCCs. We use this later to speed up reachability analysis
+			storm::storage::BitVector statesInBsccs(numOfStates);
+			
+			std::vector<uint_fast64_t> stateToBsccIndexMap(transitionMatrix.getColumnCount());
+
+			for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
+				storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
+
+				// Gather information for later use.
+				for (auto const& state : bscc) {
+					statesInBsccs.set(state);
+					stateToBsccIndexMap[state] = currentBsccIndex;
+				}
+			}
+
+			storm::storage::BitVector statesNotInBsccs = ~statesInBsccs;
+
+			// calculate steady state distribution for all BSCCs by calculating an eigenvector for the eigenvalue 1 of the transposed transition matrix for the bsccs
+			storm::storage::SparseMatrix<ValueType> bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true);
+
+			//subtract identity matrix
+			for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) {
+				for (auto& entry : bsccEquationSystem.getRow(row)) {
+					if (entry.getColumn() == row) {
+						entry.setValue(entry.getValue() - one);
+					}
+				}
+			}
+			//now transpose, this internally removes all explicit zeros from the matrix that where introduced when subtracting the identity matrix
+			bsccEquationSystem = bsccEquationSystem.transpose();
+
+			std::vector<ValueType> bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero);
+			std::vector<ValueType> bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), one);
+			{
+				auto solver = this->linearEquationSolverFactory->create(bsccEquationSystem);
+				solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide);
+			}
+
+			//calculate LRA Value for each BSCC from steady state distribution in BSCCs
+			// we have to do some scaling, as the probabilities for each BSCC have to sum up to one, which they don't necessarily do in the solution of the equation system
+			std::vector<ValueType> bsccLra(bsccDecomposition.size(), zero);
+			std::vector<ValueType> bsccTotalValue(bsccDecomposition.size(), zero);
+			size_t i = 0;
+			for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) {
+				if (psiStates.get(*stateIter)) {
+					bsccLra[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i];
+				}
+				bsccTotalValue[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i];
+			}
+			for (i = 0; i < bsccLra.size(); ++i) {
+				bsccLra[i] /= bsccTotalValue[i];
+			}
+
+			//calculate LRA for states not in bsccs as expected reachability rewards
+			//target states are states in bsccs, transition reward is the lra of the bscc for each transition into a bscc and 0 otherwise
+			//this corresponds to sum of LRAs in BSCC weighted by the reachability probability of the BSCC
+
+			std::vector<ValueType> rewardRightSide;
+			rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits());
+
+			for (auto state : statesNotInBsccs) {
+				ValueType reward = zero;
+				for (auto entry : transitionMatrix.getRow(state)) {
+					if (statesInBsccs.get(entry.getColumn())) {
+						reward += entry.getValue() * bsccLra[stateToBsccIndexMap[entry.getColumn()]];
+					}
+				}
+				rewardRightSide.push_back(reward);
+			}
+
+			storm::storage::SparseMatrix<ValueType> rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true);
+			rewardEquationSystemMatrix.convertToEquationSystem();
+
+			std::vector<ValueType> rewardSolution(rewardEquationSystemMatrix.getColumnCount(), one);
+
+			{
+				auto solver = this->linearEquationSolverFactory->create(rewardEquationSystemMatrix);
+				solver->solveEquationSystem(rewardSolution, rewardRightSide);
+			}
+
+			// now fill the result vector
+			std::vector<ValueType> result(numOfStates);
+
+			auto rewardSolutionIter = rewardSolution.begin();
+			for (size_t state = 0; state < numOfStates; ++state) {
+				if (statesInBsccs.get(state)) {
+					//assign the value of the bscc the state is in
+					result[state] = bsccLra[stateToBsccIndexMap[state]];
+				} else {
+					assert(rewardSolutionIter != rewardSolution.end());
+					//take the value from the reward computation
+					//since the n-th state not in any bscc is the n-th entry in rewardSolution we can just take the next value from the iterator
+					result[state] = *rewardSolutionIter;
+					++rewardSolutionIter;
+				}
+			}
+
+			return result;
+
+			//old implementeation
+
+			//now we calculate the long run average for each BSCC in isolation and compute the weighted contribution of the BSCC to the LRA value of all states
+			//for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
+			//	storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
+
+			//	storm::storage::BitVector statesInThisBscc(numOfStates);
+			//	for (auto const& state : bscc) {
+			//		statesInThisBscc.set(state);
+			//	}
+
+			//	//ValueType lraForThisBscc = this->computeLraForBSCC(transitionMatrix, psiStates, bscc);
+			//	ValueType lraForThisBscc = bsccLra[currentBsccIndex];
+
+			//	//the LRA value of a BSCC contributes to the LRA value of a state with the probability of reaching that BSCC from that state
+			//	//thus we add Prob(Eventually statesInThisBscc) * lraForThisBscc to the result vector
+
+			//	//the reachability probabilities will be zero in other BSCCs, thus we can set the left operand of the until formula to statesNotInBsccs as an optimization
+			//	std::vector<ValueType> reachProbThisBscc = this->computeUntilProbabilitiesHelper(statesNotInBsccs, statesInThisBscc, false);
+			//	
+			//	storm::utility::vector::scaleVectorInPlace<ValueType>(reachProbThisBscc, lraForThisBscc);
+			//	storm::utility::vector::addVectorsInPlace<ValueType>(result, reachProbThisBscc);
+			//}
+
+			//return result;
+		}
+
+		template<typename ValueType>
+		std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+			std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
+			ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
+
+			return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeLongRunAverageHelper(subResult.getTruthValuesVector(), qualitative)));
+		}
+
+
+		template<typename ValueType>
+		ValueType SparseDtmcPrctlModelChecker<ValueType>::computeLraForBSCC(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+			//if size is one we can immediately derive the result
+			if (bscc.size() == 1){
+				if (psiStates.get(*(bscc.begin()))) {
+					return storm::utility::one<ValueType>();
+				} else{
+					return storm::utility::zero<ValueType>();
+				}
+			}
+
+			storm::storage::BitVector subsystem = storm::storage::BitVector(transitionMatrix.getRowCount());
+			subsystem.set(bscc.begin(), bscc.end());
+
+			//we now have to solve ((P')^t - I) * x = 0, where P' is the submatrix of transitionMatrix,
+			// ^t means transose, and I is the identity matrix.
+			
+			storm::storage::SparseMatrix<ValueType> subsystemMatrix = transitionMatrix.getSubmatrix(false, subsystem, subsystem, true);
+			subsystemMatrix = subsystemMatrix.transpose();
+
+			// subtract 1 on the diagonal and at the same time add a row with all ones to enforce that the result of the equation system is unique
+			storm::storage::SparseMatrixBuilder<ValueType> equationSystemBuilder(subsystemMatrix.getRowCount() + 1, subsystemMatrix.getColumnCount(), subsystemMatrix.getEntryCount() + subsystemMatrix.getColumnCount());
+			ValueType one = storm::utility::one<ValueType>();
+			ValueType zero = storm::utility::zero<ValueType>();
+			bool foundDiagonalElement = false;
+			for (uint_fast64_t row = 0; row < subsystemMatrix.getRowCount(); ++row) {
+				for (auto& entry : subsystemMatrix.getRowGroup(row)) {
+					if (entry.getColumn() == row) {
+						equationSystemBuilder.addNextValue(row, entry.getColumn(), entry.getValue() - one);
+						foundDiagonalElement = true;
+					} else {
+						equationSystemBuilder.addNextValue(row, entry.getColumn(), entry.getValue());
+					}
+				}
+
+				// Throw an exception if a row did not have an element on the diagonal.
+				STORM_LOG_THROW(foundDiagonalElement, storm::exceptions::InvalidOperationException, "Internal Error, no diagonal entry found.");
+			}
+			for (uint_fast64_t column = 0; column + 1 < subsystemMatrix.getColumnCount(); ++column) {
+				equationSystemBuilder.addNextValue(subsystemMatrix.getRowCount(), column, one);
+			}
+			equationSystemBuilder.addNextValue(subsystemMatrix.getRowCount(), subsystemMatrix.getColumnCount() - 1, zero);
+			subsystemMatrix = equationSystemBuilder.build();
+
+			// create x and b for the equation system. setting the last entry of b to one enforces the sum over the unique solution vector is one
+			std::vector<ValueType> steadyStateDist(subsystemMatrix.getRowCount(), zero);
+			std::vector<ValueType> b(subsystemMatrix.getRowCount(), zero);
+			b[subsystemMatrix.getRowCount() - 1] = one;
+
+
+			{
+				auto solver = linearEquationSolverFactory.create(subsystemMatrix);
+				solver->solveEquationSystem(steadyStateDist, b);
+			}
+
+			//remove the last entry of the vector as it was just there to enforce the unique solution
+			steadyStateDist.pop_back();
+			
+			//calculate the fraction we spend in psi states on the long run
+			std::vector<ValueType> steadyStateForPsiStates(transitionMatrix.getRowCount() - 1, zero);
+			storm::utility::vector::setVectorValues(steadyStateForPsiStates, psiStates & subsystem, steadyStateDist);
+
+			ValueType result = zero;
+
+			for (auto value : steadyStateForPsiStates) {
+				result += value;
+			}
+
+			return result;
+		}
+
         template<typename ValueType>
         storm::models::sparse::Dtmc<ValueType> const& SparseDtmcPrctlModelChecker<ValueType>::getModel() const {
             return this->template getModelAs<storm::models::sparse::Dtmc<ValueType>>();
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
index e5407f9b5..b55f00faa 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
@@ -5,6 +5,7 @@
 #include "src/models/sparse/Dtmc.h"
 #include "src/utility/solver.h"
 #include "src/solver/LinearEquationSolver.h"
+#include "src/storage/StronglyConnectedComponent.h"
 
 namespace storm {
     namespace modelchecker {
@@ -27,6 +28,7 @@ namespace storm {
             virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+			virtual std::unique_ptr<CheckResult> computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             
         protected:
             storm::models::sparse::Dtmc<ValueType> const& getModel() const override;
@@ -39,7 +41,10 @@ namespace storm {
             std::vector<ValueType> computeInstantaneousRewardsHelper(uint_fast64_t stepCount) const;
             std::vector<ValueType> computeCumulativeRewardsHelper(uint_fast64_t stepBound) const;
             static std::vector<ValueType> computeReachabilityRewardsHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative);
-            
+			std::vector<ValueType> computeLongRunAverageHelper(storm::storage::BitVector const& psiStates, bool qualitative) const;
+
+			static ValueType computeLraForBSCC(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+
             // An object that is used for retrieving linear equation solvers.
             std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
         };
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
index bc4bb8a8e..0a9007c22 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
@@ -1,6 +1,7 @@
 #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
 
 #include <vector>
+#include <memory>
 
 #include "src/utility/constants.h"
 #include "src/utility/macros.h"
@@ -10,8 +11,13 @@
 #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 
+#include "src/solver/LpSolver.h"
+
 #include "src/exceptions/InvalidStateException.h"
 #include "src/exceptions/InvalidPropertyException.h"
+#include "src/storage/expressions/Expressions.h"
+
+#include "src/storage/MaximalEndComponentDecomposition.h"
 
 namespace storm {
     namespace modelchecker {
@@ -316,6 +322,213 @@ namespace storm {
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->MinMaxLinearEquationSolverFactory, qualitative)));
         }
         
+
+		template<typename ValueType>
+		std::vector<ValueType> SparseMdpPrctlModelChecker<ValueType>::computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const {
+			// If there are no goal states, we avoid the computation and directly return zero.
+			auto numOfStates = this->getModel().getNumberOfStates();
+			if (psiStates.empty()) {
+				return std::vector<ValueType>(numOfStates, storm::utility::zero<ValueType>());
+			}
+
+			// Likewise, if all bits are set, we can avoid the computation and set.
+			if ((~psiStates).empty()) {
+				return std::vector<ValueType>(numOfStates, storm::utility::one<ValueType>());
+			}
+
+			// Start by decomposing the MDP into its MECs.
+			storm::storage::MaximalEndComponentDecomposition<double> mecDecomposition(this->getModel());
+
+			// Get some data members for convenience.
+			typename storm::storage::SparseMatrix<ValueType> const& transitionMatrix = this->getModel().getTransitionMatrix();
+			std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
+			ValueType one = storm::utility::one<ValueType>();
+			ValueType zero = storm::utility::zero<ValueType>();
+
+			//first calculate LRA for the Maximal End Components.
+			storm::storage::BitVector statesInMecs(numOfStates);
+			std::vector<uint_fast64_t> stateToMecIndexMap(transitionMatrix.getColumnCount());
+			std::vector<ValueType> lraValuesForEndComponents(mecDecomposition.size(), zero);
+
+			for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) {
+				storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex];
+
+				lraValuesForEndComponents[currentMecIndex] = computeLraForMaximalEndComponent(minimize, transitionMatrix, psiStates, mec);
+
+				// Gather information for later use.
+				for (auto const& stateChoicesPair : mec) {
+					statesInMecs.set(stateChoicesPair.first);
+					stateToMecIndexMap[stateChoicesPair.first] = currentMecIndex;
+				}
+			}
+
+			// For fast transition rewriting, we build some auxiliary data structures.
+			storm::storage::BitVector statesNotContainedInAnyMec = ~statesInMecs;
+			uint_fast64_t firstAuxiliaryStateIndex = statesNotContainedInAnyMec.getNumberOfSetBits();
+			uint_fast64_t lastStateNotInMecs = 0;
+			uint_fast64_t numberOfStatesNotInMecs = 0;
+			std::vector<uint_fast64_t> statesNotInMecsBeforeIndex;
+			statesNotInMecsBeforeIndex.reserve(this->getModel().getNumberOfStates());
+			for (auto state : statesNotContainedInAnyMec) {
+				while (lastStateNotInMecs <= state) {
+					statesNotInMecsBeforeIndex.push_back(numberOfStatesNotInMecs);
+					++lastStateNotInMecs;
+				}
+				++numberOfStatesNotInMecs;
+			}
+
+			// Finally, we are ready to create the SSP matrix and right-hand side of the SSP.
+			std::vector<ValueType> b;
+			typename storm::storage::SparseMatrixBuilder<ValueType> sspMatrixBuilder(0, 0, 0, false, true, numberOfStatesNotInMecs + mecDecomposition.size());
+
+			// If the source state is not contained in any MEC, we copy its choices (and perform the necessary modifications).
+			uint_fast64_t currentChoice = 0;
+			for (auto state : statesNotContainedInAnyMec) {
+				sspMatrixBuilder.newRowGroup(currentChoice);
+
+				for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice, ++currentChoice) {
+					std::vector<ValueType> auxiliaryStateToProbabilityMap(mecDecomposition.size());
+					b.push_back(storm::utility::zero<ValueType>());
+
+					for (auto element : transitionMatrix.getRow(choice)) {
+						if (statesNotContainedInAnyMec.get(element.getColumn())) {
+							// If the target state is not contained in an MEC, we can copy over the entry.
+							sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue());
+						} else {
+							// If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector
+							// so that we are able to write the cumulative probability to the MEC into the matrix.
+							auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue();
+						}
+					}
+
+					// Now insert all (cumulative) probability values that target an MEC.
+					for (uint_fast64_t mecIndex = 0; mecIndex < auxiliaryStateToProbabilityMap.size(); ++mecIndex) {
+						if (auxiliaryStateToProbabilityMap[mecIndex] != 0) {
+							sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + mecIndex, auxiliaryStateToProbabilityMap[mecIndex]);
+						}
+					}
+				}
+			}
+
+			// Now we are ready to construct the choices for the auxiliary states.
+			for (uint_fast64_t mecIndex = 0; mecIndex < mecDecomposition.size(); ++mecIndex) {
+				storm::storage::MaximalEndComponent const& mec = mecDecomposition[mecIndex];
+				sspMatrixBuilder.newRowGroup(currentChoice);
+
+				for (auto const& stateChoicesPair : mec) {
+					uint_fast64_t state = stateChoicesPair.first;
+					boost::container::flat_set<uint_fast64_t> const& choicesInMec = stateChoicesPair.second;
+
+					for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) {
+						std::vector<ValueType> auxiliaryStateToProbabilityMap(mecDecomposition.size());
+
+						// If the choice is not contained in the MEC itself, we have to add a similar distribution to the auxiliary state.
+						if (choicesInMec.find(choice) == choicesInMec.end()) {
+							b.push_back(storm::utility::zero<ValueType>());
+
+							for (auto element : transitionMatrix.getRow(choice)) {
+								if (statesNotContainedInAnyMec.get(element.getColumn())) {
+									// If the target state is not contained in an MEC, we can copy over the entry.
+									sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue());
+								} else {
+									// If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector
+									// so that we are able to write the cumulative probability to the MEC into the matrix.
+									auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue();
+								}
+							}
+
+							// Now insert all (cumulative) probability values that target an MEC.
+							for (uint_fast64_t targetMecIndex = 0; targetMecIndex < auxiliaryStateToProbabilityMap.size(); ++targetMecIndex) {
+								if (auxiliaryStateToProbabilityMap[targetMecIndex] != 0) {
+									sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]);
+								}
+							}
+
+							++currentChoice;
+						}
+					}
+				}
+
+				// For each auxiliary state, there is the option to achieve the reward value of the LRA associated with the MEC.
+				++currentChoice;
+				b.push_back(lraValuesForEndComponents[mecIndex]);
+			}
+
+			// Finalize the matrix and solve the corresponding system of equations.
+			storm::storage::SparseMatrix<ValueType> sspMatrix = sspMatrixBuilder.build(currentChoice);
+
+			std::vector<ValueType> sspResult(numberOfStatesNotInMecs + mecDecomposition.size());
+			std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory->create(sspMatrix);
+			solver->solveEquationSystem(minimize, sspResult, b);
+
+			// Prepare result vector.
+			std::vector<ValueType> result(this->getModel().getNumberOfStates());
+
+			// Set the values for states not contained in MECs.
+			storm::utility::vector::setVectorValues(result, statesNotContainedInAnyMec, sspResult);
+
+			// Set the values for all states in MECs.
+			for (auto state : statesInMecs) {
+				result[state] = sspResult[firstAuxiliaryStateIndex + stateToMecIndexMap[state]];
+			}
+
+			return result;
+		}
+
+		template<typename ValueType>
+		std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+			STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
+			
+			std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
+			ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
+			
+			return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeLongRunAverageHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector(), qualitative)));
+		}
+
+		template<typename ValueType>
+		ValueType SparseMdpPrctlModelChecker<ValueType>::computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::MaximalEndComponent const& mec) {
+			std::shared_ptr<storm::solver::LpSolver> solver = storm::utility::solver::getLpSolver("LRA for MEC");
+			solver->setModelSense(minimize ? storm::solver::LpSolver::ModelSense::Maximize : storm::solver::LpSolver::ModelSense::Minimize);
+
+			//// First, we need to create the variables for the problem.
+			std::map<uint_fast64_t, storm::expressions::Variable> stateToVariableMap;
+			for (auto const& stateChoicesPair : mec) {
+				std::string variableName = "h" + std::to_string(stateChoicesPair.first);
+				stateToVariableMap[stateChoicesPair.first] = solver->addUnboundedContinuousVariable(variableName);
+			}
+			storm::expressions::Variable lambda = solver->addUnboundedContinuousVariable("L", 1);
+			solver->update();
+
+			// Now we encode the problem as constraints.
+			for (auto const& stateChoicesPair : mec) {
+				uint_fast64_t state = stateChoicesPair.first;
+
+				// Now, based on the type of the state, create a suitable constraint.
+				for (auto choice : stateChoicesPair.second) {
+					storm::expressions::Expression constraint = -lambda;
+					ValueType r = 0;
+
+					for (auto element : transitionMatrix.getRow(choice)) {
+						constraint = constraint + stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue());
+						if (psiStates.get(element.getColumn())) {
+							r += element.getValue();
+						}
+					}
+					constraint = solver->getConstant(r) + constraint;
+
+					if (minimize) {
+						constraint = stateToVariableMap.at(state) <= constraint;
+					} else {
+						constraint = stateToVariableMap.at(state) >= constraint;
+					}
+					solver->addConstraint("state" + std::to_string(state) + "," + std::to_string(choice), constraint);
+				}
+			}
+
+			solver->optimize();
+			return solver->getContinuousValue(lambda);
+		}
+
         template<typename ValueType>
         storm::models::sparse::Mdp<ValueType> const& SparseMdpPrctlModelChecker<ValueType>::getModel() const {
             return this->template getModelAs<storm::models::sparse::Mdp<ValueType>>();
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
index 6246b846c..02be4a47e 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
@@ -5,6 +5,8 @@
 #include "src/models/sparse/Mdp.h"
 #include "src/utility/solver.h"
 #include "src/solver/MinMaxLinearEquationSolver.h"
+#include "src/storage/MaximalEndComponent.h"
+
 
 namespace storm {
     namespace counterexamples {
@@ -39,6 +41,7 @@ namespace storm {
             virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+			virtual std::unique_ptr<CheckResult> computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             
         protected:
             storm::models::sparse::Mdp<ValueType> const& getModel() const override;
@@ -52,7 +55,10 @@ namespace storm {
             std::vector<ValueType> computeInstantaneousRewardsHelper(bool minimize, uint_fast64_t stepCount) const;
             std::vector<ValueType> computeCumulativeRewardsHelper(bool minimize, uint_fast64_t stepBound) const;
             std::vector<ValueType> computeReachabilityRewardsHelper(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& MinMaxLinearEquationSolverFactory, bool qualitative) const;
+			std::vector<ValueType> computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const;
             
+			static ValueType computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::MaximalEndComponent const& mec);
+
             // An object that is used for retrieving solvers for systems of linear equations that are the result of nondeterministic choices.
             std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>> MinMaxLinearEquationSolverFactory;
         };
diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp
index 880d15d82..58a842cd3 100644
--- a/src/storage/SparseMatrix.cpp
+++ b/src/storage/SparseMatrix.cpp
@@ -239,20 +239,12 @@ namespace storm {
         
         template<typename ValueType>
         SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications, std::vector<MatrixEntry<index_type, ValueType>> const& columnsAndValues, std::vector<index_type> const& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), rowGroupIndices(rowGroupIndices) {
-            for (auto const& element : *this) {
-                if (element.getValue() != storm::utility::zero<ValueType>()) {
-                    ++this->nonzeroEntryCount;
-                }
-            }
+			this->updateNonzeroEntryCount();
         }
         
         template<typename ValueType>
         SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type>&& rowIndications, std::vector<MatrixEntry<index_type, ValueType>>&& columnsAndValues, std::vector<index_type>&& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(std::move(columnsAndValues)), rowIndications(std::move(rowIndications)), rowGroupIndices(std::move(rowGroupIndices)) {
-            for (auto const& element : *this) {
-                if (element.getValue() != storm::utility::zero<ValueType>()) {
-                    ++this->nonzeroEntryCount;
-                }
-            }
+			this->updateNonzeroEntryCount();
         }
                 
         template<typename ValueType>
@@ -366,7 +358,23 @@ namespace storm {
         template<typename ValueType>
         typename SparseMatrix<ValueType>::index_type SparseMatrix<ValueType>::getNonzeroEntryCount() const {
             return nonzeroEntryCount;
-        }
+		}
+
+		template<typename ValueType>
+		void SparseMatrix<ValueType>::updateNonzeroEntryCount() const {
+			//SparseMatrix<ValueType>* self = const_cast<SparseMatrix<ValueType>*>(this);
+			this->nonzeroEntryCount = 0;
+			for (auto const& element : *this) {
+				if (element.getValue() != storm::utility::zero<ValueType>()) {
+					++this->nonzeroEntryCount;
+				}
+			}
+		}
+
+		template<typename ValueType>
+		void SparseMatrix<ValueType>::updateNonzeroEntryCount(std::make_signed<index_type>::type difference) {
+			this->nonzeroEntryCount += difference;
+		}
         
         template<typename ValueType>
         typename SparseMatrix<ValueType>::index_type SparseMatrix<ValueType>::getRowGroupCount() const {
@@ -416,6 +424,7 @@ namespace storm {
             columnValuePtr->setValue(storm::utility::one<ValueType>());
             ++columnValuePtr;
             for (; columnValuePtr != columnValuePtrEnd; ++columnValuePtr) {
+				++this->nonzeroEntryCount;
                 columnValuePtr->setColumn(0);
                 columnValuePtr->setValue(storm::utility::zero<ValueType>());
             }
@@ -603,10 +612,15 @@ namespace storm {
         }
         
         template <typename ValueType>
-        SparseMatrix<ValueType> SparseMatrix<ValueType>::transpose(bool joinGroups) const {
+        SparseMatrix<ValueType> SparseMatrix<ValueType>::transpose(bool joinGroups, bool keepZeros) const {
             index_type rowCount = this->getColumnCount();
             index_type columnCount = joinGroups ? this->getRowGroupCount() : this->getRowCount();
-            index_type entryCount = this->getEntryCount();
+			if (keepZeros) {
+				index_type entryCount = this->getEntryCount();
+			} else {
+				this->updateNonzeroEntryCount();
+				index_type entryCount = this->getNonzeroEntryCount();
+			}
             
             std::vector<index_type> rowIndications(rowCount + 1);
             std::vector<MatrixEntry<index_type, ValueType>> columnsAndValues(entryCount);
@@ -614,7 +628,7 @@ namespace storm {
             // First, we need to count how many entries each column has.
             for (index_type group = 0; group < columnCount; ++group) {
                 for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) {
-                    if (transition.getValue() != storm::utility::zero<ValueType>()) {
+					if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) {
                         ++rowIndications[transition.getColumn() + 1];
                     }
                 }
@@ -633,7 +647,7 @@ namespace storm {
             // Now we are ready to actually fill in the values of the transposed matrix.
             for (index_type group = 0; group < columnCount; ++group) {
                 for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) {
-                    if (transition.getValue() != storm::utility::zero<ValueType>()) {
+					if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) {
                         columnsAndValues[nextIndices[transition.getColumn()]] = std::make_pair(group, transition.getValue());
                         nextIndices[transition.getColumn()]++;
                     }
@@ -659,13 +673,22 @@ namespace storm {
         template<typename ValueType>
         void SparseMatrix<ValueType>::invertDiagonal() {
             // Now iterate over all row groups and set the diagonal elements to the inverted value.
-            // If there is a row without the diagonal element, an exception is thrown.
-            ValueType one = storm::utility::one<ValueType>();
+			// If there is a row without the diagonal element, an exception is thrown.
+			ValueType one = storm::utility::one<ValueType>();
+			ValueType zero = storm::utility::zero<ValueType>();
             bool foundDiagonalElement = false;
             for (index_type group = 0; group < this->getRowGroupCount(); ++group) {
                 for (auto& entry : this->getRowGroup(group)) {
                     if (entry.getColumn() == group) {
-                        entry.setValue(one - entry.getValue());
+						if (entry.getValue() == one) {
+							--this->nonzeroEntryCount;
+							entry.setValue(zero);
+						} else if (entry.getValue() == zero) {
+							++this->nonzeroEntryCount;
+							entry.setValue(one);
+						} else {
+							entry.setValue(one - entry.getValue());
+						}
                         foundDiagonalElement = true;
                     }
                 }
@@ -695,6 +718,7 @@ namespace storm {
             for (index_type group = 0; group < this->getRowGroupCount(); ++group) {
                 for (auto& entry : this->getRowGroup(group)) {
                     if (entry.getColumn() == group) {
+						--this->nonzeroEntryCount;
                         entry.setValue(storm::utility::zero<ValueType>());
                     }
                 }
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index 8eded966b..45e7b19ae 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -465,12 +465,26 @@ namespace storm {
 			*/
 			uint_fast64_t getRowGroupEntryCount(uint_fast64_t const group) const;
 
-            /*!
-             * Returns the number of nonzero entries in the matrix.
-             *
-             * @return The number of nonzero entries in the matrix.
-             */
-            index_type getNonzeroEntryCount() const;
+			/*!
+			* Returns the cached number of nonzero entries in the matrix.
+			*
+			* @see updateNonzeroEntryCount()
+			*
+			* @return The number of nonzero entries in the matrix.
+			*/
+			index_type getNonzeroEntryCount() const;
+
+			/*!
+			* Recompute the nonzero entry count
+			*/
+			void updateNonzeroEntryCount() const;
+
+			/*!
+			* Change the nonzero entry count by the provided value.
+			*
+			* @param difference Difference between old and new nonzero entry count.
+			*/
+			void updateNonzeroEntryCount(std::make_signed<index_type>::type difference);
             
             /*!
              * Returns the number of row groups in the matrix.
@@ -574,12 +588,13 @@ namespace storm {
             
             /*!
              * Transposes the matrix.
-             *
-             * @param joinGroups A flag indicating whether the row groups are supposed to be treated as single rows.
+			 *
+			 * @param joinGroups A flag indicating whether the row groups are supposed to be treated as single rows.
+			 * @param keepZeros A flag indicating whether entries with value zero should be kept.
              *
              * @return A sparse matrix that represents the transpose of this matrix.
              */
-            storm::storage::SparseMatrix<value_type> transpose(bool joinGroups = false) const;
+			storm::storage::SparseMatrix<value_type> transpose(bool joinGroups = false, bool keepZeros = false) const;
             
             /*!
              * Transforms the matrix into an equation system. That is, it transforms the matrix A into a matrix (1-A).
@@ -826,7 +841,7 @@ namespace storm {
             index_type entryCount;
             
             // The number of nonzero entries in the matrix.
-            index_type nonzeroEntryCount;
+            mutable index_type nonzeroEntryCount;
             
             // The storage for the columns and values of all entries in the matrix.
             std::vector<MatrixEntry<index_type, value_type>> columnsAndValues;
diff --git a/src/utility/cli.h b/src/utility/cli.h
index 9966f8bc8..f56feaf7c 100644
--- a/src/utility/cli.h
+++ b/src/utility/cli.h
@@ -429,7 +429,7 @@ namespace storm {
                         } else {
                             storm::modelchecker::SparseDtmcEliminationModelChecker<ValueType> modelchecker2(*dtmc);
                             if (modelchecker2.canHandle(*formula.get())) {
-                                modelchecker2.check(*formula.get());
+                                result = modelchecker2.check(*formula.get());
                             }
                         }
                     } else if (model->getType() == storm::models::ModelType::Mdp) {
diff --git a/src/utility/graph.h b/src/utility/graph.h
index 6e950bcd2..3c04feb86 100644
--- a/src/utility/graph.h
+++ b/src/utility/graph.h
@@ -46,7 +46,7 @@ namespace storm {
                     currentState = stack.back();
                     stack.pop_back();
                     
-                    for (auto const& successor : transitionMatrix.getRow(currentState)) {
+                    for (auto const& successor : transitionMatrix.getRowGroup(currentState)) {
                         // Only explore the state if the transition was actually there and the successor has not yet
                         // been visited.
                         if (successor.getValue() != storm::utility::zero<T>() && !reachableStates.get(successor.getColumn())) {
diff --git a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp
index 6a814d63d..5c637afcb 100644
--- a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp
+++ b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp
@@ -22,7 +22,11 @@ TEST(GmmxxCtmcCslModelCheckerTest, Cluster) {
     std::shared_ptr<storm::logic::Formula> formula(nullptr);
     
     // Build the model.
-    typename storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#ifdef WINDOWS
+    storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#else
+	typename storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#endif
     options.buildRewards = true;
     options.rewardModelName = "num_repairs";
     std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>::translateProgram(program, options);
@@ -80,7 +84,12 @@ TEST(GmmxxCtmcCslModelCheckerTest, Embedded) {
     std::shared_ptr<storm::logic::Formula> formula(nullptr);
     
     // Build the model.
-    typename storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#ifdef WINDOWS
+	storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#else
+	typename storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#endif
+
     options.buildRewards = true;
     options.rewardModelName = "up";
     std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>::translateProgram(program, options);
@@ -172,8 +181,12 @@ TEST(GmmxxCtmcCslModelCheckerTest, Tandem) {
     std::shared_ptr<storm::logic::Formula> formula(nullptr);
     
     // Build the model.
-    typename storm::builder::ExplicitPrismModelBuilder<double>::Options options;
-    options.buildRewards = true;
+#ifdef WINDOWS
+	storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#else
+	typename storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#endif
+	options.buildRewards = true;
     options.rewardModelName = "customers";
     std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>::translateProgram(program, options);
     ASSERT_EQ(storm::models::ModelType::Ctmc, model->getType());
diff --git a/test/functional/modelchecker/GmmxxHybridDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxHybridDtmcPrctlModelCheckerTest.cpp
index cd1a014a9..df17ce336 100644
--- a/test/functional/modelchecker/GmmxxHybridDtmcPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/GmmxxHybridDtmcPrctlModelCheckerTest.cpp
@@ -16,7 +16,11 @@ TEST(GmmxxHybridDtmcPrctlModelCheckerTest, Die) {
     storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/die.pm");
 
     // Build the die model with its reward model.
-    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#ifdef WINDOWS
+    storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#else
+	typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#endif
     options.buildRewards = true;
     options.rewardModelName = "coin_flips";
     std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
@@ -118,7 +122,11 @@ TEST(GmmxxHybridDtmcPrctlModelCheckerTest, SynchronousLeader) {
     storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/leader-3-5.pm");
     
     // Build the die model with its reward model.
-    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#ifdef WINDOWS
+    storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#else
+	typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+#endif
     options.buildRewards = true;
     options.rewardModelName = "num_rounds";
     std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
diff --git a/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp
index 3dfa1b6be..9db619115 100644
--- a/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp
+++ b/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp
@@ -22,7 +22,11 @@ TEST(SparseCtmcCslModelCheckerTest, Cluster) {
     std::shared_ptr<storm::logic::Formula> formula(nullptr);
     
     // Build the model.
-    typename storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#ifdef WINDOWS
+    storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#else
+	typename storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#endif
     options.buildRewards = true;
     options.rewardModelName = "num_repairs";
     std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>::translateProgram(program, options);
@@ -80,7 +84,11 @@ TEST(SparseCtmcCslModelCheckerTest, Embedded) {
     std::shared_ptr<storm::logic::Formula> formula(nullptr);
     
     // Build the model.
-    typename storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#ifdef WINDOWS
+    storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#else
+	typename storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#endif
     options.buildRewards = true;
     options.rewardModelName = "up";
     std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>::translateProgram(program, options);
@@ -172,7 +180,11 @@ TEST(SparseCtmcCslModelCheckerTest, Tandem) {
     std::shared_ptr<storm::logic::Formula> formula(nullptr);
     
     // Build the model with the customers reward structure.
-    typename storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#ifdef WINDOWS
+    storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#else
+	typename storm::builder::ExplicitPrismModelBuilder<double>::Options options;
+#endif
     options.buildRewards = true;
     options.rewardModelName = "customers";
     std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>::translateProgram(program, options);
diff --git a/test/functional/modelchecker/SparseDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SparseDtmcPrctlModelCheckerTest.cpp
index c204d36dd..7f53f45b4 100644
--- a/test/functional/modelchecker/SparseDtmcPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/SparseDtmcPrctlModelCheckerTest.cpp
@@ -127,3 +127,148 @@ TEST(SparseDtmcPrctlModelCheckerTest, SynchronousLeader) {
     
     EXPECT_NEAR(1.0448979589010925, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision());
 }
+
+TEST(SparseDtmcPrctlModelCheckerTest, LRA_SingleBscc) {
+	storm::storage::SparseMatrixBuilder<double> matrixBuilder;
+	std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc;
+
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 2);
+		matrixBuilder.addNextValue(0, 1, 1.);
+		matrixBuilder.addNextValue(1, 0, 1.);
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(2);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 1);
+
+		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
+
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 4);
+		matrixBuilder.addNextValue(0, 0, .5);
+		matrixBuilder.addNextValue(0, 1, .5);
+		matrixBuilder.addNextValue(1, 0, .5);
+		matrixBuilder.addNextValue(1, 1, .5);
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(2);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 1);
+
+		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
+
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(3, 3, 3);
+		matrixBuilder.addNextValue(0, 1, 1);
+		matrixBuilder.addNextValue(1, 2, 1);
+		matrixBuilder.addNextValue(2, 0, 1);
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(3);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 2);
+
+		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
+
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(1. / 3., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1. / 3., quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1. / 3., quantitativeResult1[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+}
+
+TEST(SparseDtmcPrctlModelCheckerTest, LRA) {
+	storm::storage::SparseMatrixBuilder<double> matrixBuilder;
+	std::shared_ptr<storm::models::sparse::Dtmc<double>> mdp;
+
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(15, 15, 20, true);
+		matrixBuilder.addNextValue(0, 1, 1);
+		matrixBuilder.addNextValue(1, 4, 0.7);
+		matrixBuilder.addNextValue(1, 6, 0.3);
+		matrixBuilder.addNextValue(2, 0, 1);
+
+		matrixBuilder.addNextValue(3, 5, 0.8);
+		matrixBuilder.addNextValue(3, 9, 0.2);
+		matrixBuilder.addNextValue(4, 3, 1);
+		matrixBuilder.addNextValue(5, 3, 1);
+
+		matrixBuilder.addNextValue(6, 7, 1);
+		matrixBuilder.addNextValue(7, 8, 1);
+		matrixBuilder.addNextValue(8, 6, 1);
+
+		matrixBuilder.addNextValue(9, 10, 1);
+		matrixBuilder.addNextValue(10, 9, 1);
+		matrixBuilder.addNextValue(11, 9, 1);
+
+		matrixBuilder.addNextValue(12, 5, 0.4);
+		matrixBuilder.addNextValue(12, 8, 0.3);
+		matrixBuilder.addNextValue(12, 11, 0.3);
+
+		matrixBuilder.addNextValue(13, 7, 0.7);
+		matrixBuilder.addNextValue(13, 12, 0.3);
+
+		matrixBuilder.addNextValue(14, 12, 1);
+
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(15);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 1);
+		ap.addLabelToState("a", 4);
+		ap.addLabelToState("a", 5);
+		ap.addLabelToState("a", 7);
+		ap.addLabelToState("a", 11);
+		ap.addLabelToState("a", 13);
+		ap.addLabelToState("a", 14);
+
+		mdp.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
+
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.3 / 3., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult1[3], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1. / 3., quantitativeResult1[6], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult1[9], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.3/3., quantitativeResult1[12], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.79 / 3., quantitativeResult1[13], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.3 / 3., quantitativeResult1[14], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+}
\ No newline at end of file
diff --git a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp
index c12baf505..59b2ad23f 100644
--- a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp
@@ -145,52 +145,385 @@ TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) {
 
     storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory<double>()));
 
-    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
-    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
-    auto minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, eventuallyFormula);
+	auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+	auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+	auto minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, eventuallyFormula);
+
+	std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*minProbabilityOperatorFormula);
+	storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
 
-    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*minProbabilityOperatorFormula);
-    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
-    
 	EXPECT_NEAR(1, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
 
-    auto maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, eventuallyFormula);
-    
-    result = checker.check(*maxProbabilityOperatorFormula);
-    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
-    
+	auto maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, eventuallyFormula);
+
+	result = checker.check(*maxProbabilityOperatorFormula);
+	storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+
 	EXPECT_NEAR(1, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision());
 
-    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
-    auto trueFormula = std::make_shared<storm::logic::BooleanLiteralFormula>(true);
-    auto boundedUntilFormula = std::make_shared<storm::logic::BoundedUntilFormula>(trueFormula, labelFormula, 25);
-    minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, boundedUntilFormula);
-    
-    result = checker.check(*minProbabilityOperatorFormula);
-    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult3 = result->asExplicitQuantitativeCheckResult<double>();
-    
+	labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+	auto trueFormula = std::make_shared<storm::logic::BooleanLiteralFormula>(true);
+	auto boundedUntilFormula = std::make_shared<storm::logic::BoundedUntilFormula>(trueFormula, labelFormula, 25);
+	minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, boundedUntilFormula);
+
+	result = checker.check(*minProbabilityOperatorFormula);
+	storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult3 = result->asExplicitQuantitativeCheckResult<double>();
+
 	EXPECT_NEAR(0.0625, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision());
 
-    maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, boundedUntilFormula);
-    
-    result = checker.check(*maxProbabilityOperatorFormula);
-    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult4 = result->asExplicitQuantitativeCheckResult<double>();
-    
+	maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, boundedUntilFormula);
+
+	result = checker.check(*maxProbabilityOperatorFormula);
+	storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult4 = result->asExplicitQuantitativeCheckResult<double>();
+
 	EXPECT_NEAR(0.0625, quantitativeResult4[0], storm::settings::nativeEquationSolverSettings().getPrecision());
 
-    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
-    auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
-    auto minRewardOperatorFormula = std::make_shared<storm::logic::RewardOperatorFormula>(storm::logic::OptimalityType::Minimize, reachabilityRewardFormula);
+	labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+	auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
+	auto minRewardOperatorFormula = std::make_shared<storm::logic::RewardOperatorFormula>(storm::logic::OptimalityType::Minimize, reachabilityRewardFormula);
+
+	result = checker.check(*minRewardOperatorFormula);
+	storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult5 = result->asExplicitQuantitativeCheckResult<double>();
 
-    result = checker.check(*minRewardOperatorFormula);
-    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult5 = result->asExplicitQuantitativeCheckResult<double>();
-    
 	EXPECT_NEAR(4.285689611, quantitativeResult5[0], storm::settings::nativeEquationSolverSettings().getPrecision());
 
-    auto maxRewardOperatorFormula = std::make_shared<storm::logic::RewardOperatorFormula>(storm::logic::OptimalityType::Maximize, reachabilityRewardFormula);
+	auto maxRewardOperatorFormula = std::make_shared<storm::logic::RewardOperatorFormula>(storm::logic::OptimalityType::Maximize, reachabilityRewardFormula);
+
+	result = checker.check(*maxRewardOperatorFormula);
+	storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult6 = result->asExplicitQuantitativeCheckResult<double>();
 
-    result = checker.check(*maxRewardOperatorFormula);
-    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult6 = result->asExplicitQuantitativeCheckResult<double>();
-    
 	EXPECT_NEAR(4.285689611, quantitativeResult6[0], storm::settings::nativeEquationSolverSettings().getPrecision());
 }
+
+TEST(SparseMdpPrctlModelCheckerTest, LRA_SingleMec) {
+	storm::storage::SparseMatrixBuilder<double> matrixBuilder;
+	std::shared_ptr<storm::models::sparse::Mdp<double>> mdp;
+
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 2);
+		matrixBuilder.addNextValue(0, 1, 1.);
+		matrixBuilder.addNextValue(1, 0, 1.);
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(2);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 1);
+
+		mdp.reset(new storm::models::sparse::Mdp<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(.5, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.5, quantitativeResult2[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 4);
+		matrixBuilder.addNextValue(0, 0, .5);
+		matrixBuilder.addNextValue(0, 1, .5);
+		matrixBuilder.addNextValue(1, 0, .5);
+		matrixBuilder.addNextValue(1, 1, .5);
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(2);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 1);
+
+		mdp.reset(new storm::models::sparse::Mdp<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(.5, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.5, quantitativeResult2[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(4, 3, 4, true, true, 3);
+		matrixBuilder.newRowGroup(0);
+		matrixBuilder.addNextValue(0, 1, 1);
+		matrixBuilder.newRowGroup(1);
+		matrixBuilder.addNextValue(1, 0, 1);
+		matrixBuilder.addNextValue(2, 2, 1);
+		matrixBuilder.newRowGroup(3);
+		matrixBuilder.addNextValue(3, 0, 1);
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(3);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 2);
+		ap.addLabel("b");
+		ap.addLabelToState("b", 0);
+		ap.addLabel("c");
+		ap.addLabelToState("c", 0);
+		ap.addLabelToState("c", 2);
+
+		mdp.reset(new storm::models::sparse::Mdp<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(1. / 3., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1. / 3., quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1. / 3., quantitativeResult1[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.0, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult2[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult2[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("b");
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult3 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.5, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.5, quantitativeResult3[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.5, quantitativeResult3[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult4 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(1. / 3., quantitativeResult4[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1. / 3., quantitativeResult4[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1. / 3., quantitativeResult4[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("c");
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult5 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(2. / 3., quantitativeResult5[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(2. / 3., quantitativeResult5[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(2. / 3., quantitativeResult5[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult6 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.5, quantitativeResult6[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.5, quantitativeResult6[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.5, quantitativeResult6[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+}
+
+TEST(SparseMdpPrctlModelCheckerTest, LRA) {
+	storm::storage::SparseMatrixBuilder<double> matrixBuilder;
+	std::shared_ptr<storm::models::sparse::Mdp<double>> mdp;
+
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(4, 3, 4, true, true, 3);
+		matrixBuilder.newRowGroup(0);
+		matrixBuilder.addNextValue(0, 1, 1);
+		matrixBuilder.newRowGroup(1);
+		matrixBuilder.addNextValue(1, 1, 1);
+		matrixBuilder.addNextValue(2, 2, 1);
+		matrixBuilder.newRowGroup(3);
+		matrixBuilder.addNextValue(3, 2, 1);
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(3);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 0);
+		ap.addLabel("b");
+		ap.addLabelToState("b", 1);
+		ap.addLabel("c");
+		ap.addLabelToState("c", 2);
+
+		mdp.reset(new storm::models::sparse::Mdp<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.0, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult1[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.0, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult2[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult2[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("b");
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult3 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(1.0, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1.0, quantitativeResult3[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult3[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult4 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.0, quantitativeResult4[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult4[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult4[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("c");
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult5 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(1.0, quantitativeResult5[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1.0, quantitativeResult5[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1.0, quantitativeResult5[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult6 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.0, quantitativeResult6[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult6[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1.0, quantitativeResult6[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(22, 15, 28, true, true, 15);
+		matrixBuilder.newRowGroup(0);
+		matrixBuilder.addNextValue(0, 1, 1);
+		matrixBuilder.newRowGroup(1);
+		matrixBuilder.addNextValue(1, 0, 1);
+		matrixBuilder.addNextValue(2, 2, 1);
+		matrixBuilder.addNextValue(3, 4, 0.7);
+		matrixBuilder.addNextValue(3, 6, 0.3);
+		matrixBuilder.newRowGroup(4);
+		matrixBuilder.addNextValue(4, 0, 1);
+
+		matrixBuilder.newRowGroup(5);
+		matrixBuilder.addNextValue(5, 4, 1);
+		matrixBuilder.addNextValue(6, 5, 0.8);
+		matrixBuilder.addNextValue(6, 9, 0.2);
+		matrixBuilder.newRowGroup(7);
+		matrixBuilder.addNextValue(7, 3, 1);
+		matrixBuilder.addNextValue(8, 5, 1);
+		matrixBuilder.newRowGroup(9);
+		matrixBuilder.addNextValue(9, 3, 1);
+
+		matrixBuilder.newRowGroup(10);
+		matrixBuilder.addNextValue(10, 7, 1);
+		matrixBuilder.newRowGroup(11);
+		matrixBuilder.addNextValue(11, 6, 1);
+		matrixBuilder.addNextValue(12, 8, 1);
+		matrixBuilder.newRowGroup(13);
+		matrixBuilder.addNextValue(13, 6, 1);
+
+		matrixBuilder.newRowGroup(14);
+		matrixBuilder.addNextValue(14, 10, 1);
+		matrixBuilder.newRowGroup(15);
+		matrixBuilder.addNextValue(15, 9, 1);
+		matrixBuilder.addNextValue(16, 11, 1);
+		matrixBuilder.newRowGroup(17);
+		matrixBuilder.addNextValue(17, 9, 1);
+
+		matrixBuilder.newRowGroup(18);
+		matrixBuilder.addNextValue(18, 5, 0.4);
+		matrixBuilder.addNextValue(18, 8, 0.3);
+		matrixBuilder.addNextValue(18, 11, 0.3);
+
+		matrixBuilder.newRowGroup(19);
+		matrixBuilder.addNextValue(19, 7, 0.7);
+		matrixBuilder.addNextValue(19, 12, 0.3);
+
+		matrixBuilder.newRowGroup(20);
+		matrixBuilder.addNextValue(20, 12, 0.1);
+		matrixBuilder.addNextValue(20, 13, 0.9);
+		matrixBuilder.addNextValue(21, 12, 1);
+
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(15);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 1);
+		ap.addLabelToState("a", 4);
+		ap.addLabelToState("a", 5);
+		ap.addLabelToState("a", 7);
+		ap.addLabelToState("a", 11);
+		ap.addLabelToState("a", 13);
+		ap.addLabelToState("a", 14);
+
+		mdp.reset(new storm::models::sparse::Mdp<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(37./60., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(2./3., quantitativeResult1[3], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.5, quantitativeResult1[6], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1./3., quantitativeResult1[9], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(31./60., quantitativeResult1[12], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(101./200., quantitativeResult1[13], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(31./60., quantitativeResult1[14], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.3 / 3., quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult2[3], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1./3., quantitativeResult2[6], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult2[9], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.3 / 3., quantitativeResult2[12], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.79 / 3., quantitativeResult2[13], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.3 / 3., quantitativeResult2[14], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+}
\ No newline at end of file