diff --git a/prism/src/prism/ModelChecker.java b/prism/src/prism/ModelChecker.java
index 090ac177..44c07a98 100644
--- a/prism/src/prism/ModelChecker.java
+++ b/prism/src/prism/ModelChecker.java
@@ -32,12 +32,8 @@ import parser.ast.*;
 
 public interface ModelChecker
 {
-	public void setEngine(int e);
-	public void setOption(String option, boolean b);
-	public void setOption(String option, int i);
-	public void setOption(String option, double d);
-	public void setOption(String option, String s);
-	public Object check(Expression f) throws PrismException;
+	public Object check(Expression expr) throws PrismException;
+	public Object check(Expression expr, Filter filter) throws PrismException;
 }
 
 //------------------------------------------------------------------------------
diff --git a/prism/src/prism/NondetModelChecker.java b/prism/src/prism/NondetModelChecker.java
index 4af0e215..ee84021f 100644
--- a/prism/src/prism/NondetModelChecker.java
+++ b/prism/src/prism/NondetModelChecker.java
@@ -27,263 +27,78 @@
 package prism;
 
 import jdd.*;
-import odd.*;
 import dv.*;
 import mtbdd.*;
 import sparse.*;
 import hybrid.*;
-import parser.*;
 import parser.ast.*;
 
-// Model checker for MDPs
-
-public class NondetModelChecker implements ModelChecker
+/*
+ * Model checker for MDPs
+ */
+public class NondetModelChecker extends StateModelChecker
 {
-	// logs
-	private PrismLog mainLog;
-	private PrismLog techLog;
-
-	// options
-
-	// which engine to use
-	private int engine;
-	// termination criterion (iterative methods)
-	private int termCrit;
-	// parameter for termination criterion
-	private double termCritParam;
-	// max num iterations (iterative methods)
-	private int maxIters;
-	// flags
-	private boolean verbose; // verbose output?
-	private boolean precomp; // use 0,1 precomputation algorithms?
-	private boolean fairness; // use fairness?
-	// sparse bits info
-	private int SBMaxMem;
-	private int numSBLevels;
-
-	// properties file
-	private PropertiesFile propertiesFile;
-
-	// constant values
-	private Values constantValues;
-
-	// Expression2MTBDD object for translating expressions
-	private Expression2MTBDD expr2mtbdd;
-
-	// class-wide storage for probability in the initial state
-	private double numericalRes;
-
-	// model info
-	private NondetModel model;
-	private JDDNode trans;
-	private JDDNode trans01;
-	private JDDNode start;
-	private JDDNode reach;
-	private ODDNode odd;
-	private JDDNode nondetMask;
-	private JDDVars allDDRowVars;
-	private JDDVars allDDColVars;
-	private JDDVars allDDNondetVars;
-	private JDDVars[] varDDRowVars;
-
-	// constructor - set some defaults
-
-	public NondetModelChecker(PrismLog log1, PrismLog log2, Model m, PropertiesFile pf) throws PrismException
-	{
-		// initialise
-		mainLog = log1;
-		techLog = log2;
-		if (!(m instanceof NondetModel)) {
-			throw new PrismException("Error: wrong model type passed to NondetModelChecker.");
-		}
-		model = (NondetModel) m;
-		propertiesFile = pf;
-		trans = model.getTrans();
-		trans01 = model.getTrans01();
-		start = model.getStart();
-		reach = model.getReach();
-		odd = model.getODD();
-		nondetMask = model.getNondetMask();
-		allDDRowVars = model.getAllDDRowVars();
-		allDDColVars = model.getAllDDColVars();
-		allDDNondetVars = model.getAllDDNondetVars();
-		varDDRowVars = model.getVarDDRowVars();
-
-		// create list of all constant values needed
-		constantValues = new Values();
-		constantValues.addValues(model.getConstantValues());
-		if (pf != null)
-			constantValues.addValues(pf.getConstantValues());
-
-		// create Expression2MTBDD object
-		expr2mtbdd = new Expression2MTBDD(mainLog, techLog, model.getVarList(), varDDRowVars, constantValues);
-		expr2mtbdd.setFilter(reach);
-
-		// set up some default options
-		// (although all should be overridden before model checking)
-		engine = Prism.HYBRID;
-		termCrit = Prism.RELATIVE;
-		termCritParam = 1e-6;
-		maxIters = 10000;
-		verbose = false;
-		precomp = true;
-		fairness = true;
-		SBMaxMem = 1024;
-		numSBLevels = -1;
-	}
+	// Model (MDP)
+	protected NondetModel model;
 
-	// set engine
+	// Extra (MDP) model info
+	protected JDDNode nondetMask;
+	protected JDDVars allDDNondetVars;
 
-	public void setEngine(int e)
-	{
-		engine = e;
-	}
+	// Options (in addition to those inherited from StateModelChecker):
 
-	// set options (generic)
+	// Use 0,1 precomputation algorithms?
+	protected boolean precomp;
+	// Use fairness?
+	protected boolean fairness;
 
-	public void setOption(String option, boolean b)
-	{
-		if (option.equals("verbose")) {
-			verbose = b;
-		} else if (option.equals("precomp")) {
-			precomp = b;
-		} else if (option.equals("fairness")) {
-			fairness = b;
-		} else if (option.equals("compact")) {
-			PrismSparse.setCompact(b);
-			PrismHybrid.setCompact(b);
-		} else {
-			mainLog.println("Warning: option \"" + option + "\" not supported by NondetModelChecker.");
-		}
-	}
+	// Constructor
 
-	public void setOption(String option, int i)
+	public NondetModelChecker(Prism prism, Model m, PropertiesFile pf) throws PrismException
 	{
-		if (option.equals("termcrit")) {
-			termCrit = i;
-			PrismMTBDD.setTermCrit(i);
-			PrismSparse.setTermCrit(i);
-			PrismHybrid.setTermCrit(i);
-		} else if (option.equals("maxiters")) {
-			maxIters = i;
-			PrismMTBDD.setMaxIters(i);
-			PrismSparse.setMaxIters(i);
-			PrismHybrid.setMaxIters(i);
-		} else if (option.equals("sbmaxmem")) {
-			SBMaxMem = i;
-			PrismHybrid.setSBMaxMem(i);
-		} else if (option.equals("numsblevels")) {
-			numSBLevels = i;
-			PrismHybrid.setNumSBLevels(i);
-		} else {
-			mainLog.println("Warning: option \"" + option + "\" not supported by NondetModelChecker.");
+		// Initialise
+		super(prism, m, pf);
+		if (!(m instanceof NondetModel)) {
+			throw new PrismException("Wrong model type passed to NondetModelChecker.");
 		}
-	}
+		model = (NondetModel) m;
+		nondetMask = model.getNondetMask();
+		allDDNondetVars = model.getAllDDNondetVars();
 
-	public void setOption(String option, double d)
-	{
-		if (option.equals("termcritparam")) {
-			termCritParam = d;
-			PrismMTBDD.setTermCritParam(d);
-			PrismSparse.setTermCritParam(d);
-			PrismHybrid.setTermCritParam(d);
-		} else {
-			mainLog.println("Warning: option \"" + option + "\" not supported by NondetModelChecker.");
+		// Inherit some options from parent Prism object.
+		// Store locally and/or pass onto engines.
+		precomp = prism.getPrecomp();
+		fairness = prism.getFairness();
+		switch (engine) {
+		case Prism.MTBDD:
+			PrismMTBDD.setTermCrit(prism.getTermCrit());
+			PrismMTBDD.setTermCritParam(prism.getTermCritParam());
+			PrismMTBDD.setMaxIters(prism.getMaxIters());
+			break;
+		case Prism.SPARSE:
+			PrismSparse.setTermCrit(prism.getTermCrit());
+			PrismSparse.setTermCritParam(prism.getTermCritParam());
+			PrismSparse.setMaxIters(prism.getMaxIters());
+			PrismSparse.setCompact(prism.getCompact());
+		case Prism.HYBRID:
+			PrismHybrid.setTermCrit(prism.getTermCrit());
+			PrismHybrid.setTermCritParam(prism.getTermCritParam());
+			PrismHybrid.setMaxIters(prism.getMaxIters());
+			PrismHybrid.setCompact(prism.getCompact());
+			PrismHybrid.setSBMaxMem(prism.getSBMaxMem());
+			PrismHybrid.setNumSBLevels(prism.getNumSBLevels());
 		}
 	}
 
-	public void setOption(String option, String s)
-	{
-		mainLog.println("Warning: option \"" + option + "\" not supported by NondetModelChecker.");
-	}
-
-	// compute number of places to round solution vector to
-	// (given termination criteria, etc.)
-
-	public int getPlacesToRoundBy()
-	{
-		return 1 - (int) (Math.log(termCritParam) / Math.log(10));
-	}
-
 	// -----------------------------------------------------------------------------------
 	// Check a property, i.e. an expression
 	// -----------------------------------------------------------------------------------
 
-	public Object check(Expression expr) throws PrismException
-	{
-		long timer = 0;
-		JDDNode dd;
-		StateList states;
-		boolean b;
-		Object res;
-
-		// start timer
-		timer = System.currentTimeMillis();
-
-		// do model checking and store result
-		dd = checkExpression(expr);
-		states = new StateListMTBDD(dd, model);
-
-		// stop timer
-		timer = System.currentTimeMillis() - timer;
-
-		// print out model checking time
-		mainLog.println("\nTime for model checking: " + timer / 1000.0 + " seconds.");
-
-		// output and result of model checking depend on return type
-		if (expr.getType() == Expression.BOOLEAN) {
-
-			// print out number of satisfying states
-			mainLog.print("\nNumber of satisfying states: ");
-			mainLog.print(states.size());
-			if (states.size() == model.getNumStates()) {
-				mainLog.print(" (all)");
-			} else if (states.includes(model.getStart())) {
-				mainLog.print((model.getNumStartStates() == 1) ? " (including initial state)"
-						: " (including all initial states)");
-			} else {
-				mainLog.print((model.getNumStartStates() == 1) ? " (initial state not satisfied)"
-						: " (initial states not all satisfied)");
-			}
-			mainLog.print("\n");
-
-			// if "verbose", print out satisfying states
-			if (verbose) {
-				mainLog.print("\nSatisfying states:");
-				if (states.size() > 0) {
-					mainLog.print("\n");
-					states.print(mainLog);
-				} else {
-					mainLog.print(" (none)\n");
-				}
-			}
-
-			// result is true if all states satisfy, false otherwise
-			b = (states.size() == model.getNumStates());
-			res = b ? new Boolean(true) : new Boolean(false);
-
-			// print result
-			mainLog.print("\nResult: " + b + " (property " + (b ? "" : "not ") + "satisfied in all states)\n");
-		} else {
-			// result is value stored earlier
-			res = new Double(numericalRes);
-
-			// print result
-			mainLog.print("\nResult: " + res + "\n");
-		}
-
-		// finished with memory
-		states.clear();
-
-		// return result
-		return res;
-	}
-
 	// Check expression (recursive)
 
-	public JDDNode checkExpression(Expression expr) throws PrismException
+	public StateProbs checkExpression(Expression expr) throws PrismException
 	{
-		JDDNode res;
+		StateProbs res;
 
 		// P operator
 		if (expr instanceof ExpressionProb) {
@@ -293,18 +108,13 @@ public class NondetModelChecker implements ModelChecker
 		else if (expr instanceof ExpressionReward) {
 			res = checkExpressionReward((ExpressionReward) expr);
 		}
-		// Label
-		else if (expr instanceof ExpressionLabel) {
-			res = checkExpressionLabel((ExpressionLabel) expr);
-		}
-		// Otherwise, must be an ordinary expression
+		// Otherwise, use the superclass
 		else {
-			res = expr2mtbdd.translateExpression(expr);
+			res = super.checkExpression(expr);
 		}
 
-		// Filter out non-reachable states from solution
-		JDD.Ref(reach);
-		res = JDD.Apply(JDD.TIMES, res, reach);
+		// Filter out non-reachable states from solution (TODO: not for dv?)
+		res.filter(reach);
 
 		return res;
 	}
@@ -315,18 +125,16 @@ public class NondetModelChecker implements ModelChecker
 
 	// P operator
 
-	private JDDNode checkExpressionProb(ExpressionProb expr) throws PrismException
+	protected StateProbs checkExpressionProb(ExpressionProb expr) throws PrismException
 	{
 		Expression pb; // probability bound (expression)
-		double p = 0; // probability value (actual value)
+		double p = 0; // probability bound (actual value)
 		String relOp; // relational operator
 		boolean min; // are we finding min (true) or max (false) probs
 		PathExpression pe; // path expression
 
-		JDDNode filter, sol, tmp;
+		JDDNode sol;
 		StateProbs probs = null;
-		StateList states = null;
-		double minRes = 0, maxRes = 0;
 
 		// get info from prob operator
 		relOp = expr.getRelOp();
@@ -343,11 +151,11 @@ public class NondetModelChecker implements ModelChecker
 				mainLog.print("\nWarning: checking for probability " + relOp + " " + p
 						+ " - formula trivially satisfies all states\n");
 				JDD.Ref(reach);
-				return reach;
+				return new StateProbsMTBDD(reach, model);
 			} else if ((p == 0 && relOp.equals("<")) || (p == 1 && relOp.equals(">"))) {
 				mainLog.print("\nWarning: checking for probability " + relOp + " " + p
 						+ " - formula trivially satisfies no states\n");
-				return JDD.Constant(0);
+				return new StateProbsMTBDD(JDD.Constant(0), model);
 			}
 		}
 
@@ -362,71 +170,40 @@ public class NondetModelChecker implements ModelChecker
 			throw new PrismException("Can't use \"P=?\" for MDPs; use \"Pmin=?\" or \"Pmax=?\"");
 		}
 
-		// translate filter (if present)
-		filter = null;
-		if (expr.getFilter() != null) {
-			filter = checkExpression(expr.getFilter().getExpression());
-		}
-
-		// check if filter satisfies no states
-		if (filter != null)
-			if (filter.equals(JDD.ZERO)) {
-				// for P=? properties, this is an error
-				if (pb == null) {
-					throw new PrismException("Filter {" + expr.getFilter().getExpression() + "} in P=?[] property satisfies no states");
-				}
-				// otherwise just print a warning
-				else {
-					mainLog.println("\nWarning: Filter {" + expr.getFilter().getExpression()
-							+ "} satisfies no states and is being ignored");
-					JDD.Deref(filter);
-					filter = null;
-				}
-			}
-
 		// compute probabilities
 		pe = expr.getPathExpression();
-		try {
-			if (pe instanceof PathExpressionTemporal) {
-				if (((PathExpressionTemporal) pe).hasBounds()) {
-					switch (((PathExpressionTemporal) pe).getOperator()) {
-					case PathExpressionTemporal.P_U:
-						probs = checkProbBoundedUntil((PathExpressionTemporal) pe, min);
-						break;
-					case PathExpressionTemporal.P_F:
-						probs = checkProbBoundedFuture((PathExpressionTemporal) pe, min);
-						break;
-					case PathExpressionTemporal.P_G:
-						probs = checkProbBoundedGlobal((PathExpressionTemporal) pe, min);
-						break;
-					}
-				} else {
-					switch (((PathExpressionTemporal) pe).getOperator()) {
-					case PathExpressionTemporal.P_X:
-						probs = checkProbNext((PathExpressionTemporal) pe, min);
-						break;
-					case PathExpressionTemporal.P_U:
-						probs = checkProbUntil((PathExpressionTemporal) pe, pb, p, min);
-						break;
-					case PathExpressionTemporal.P_F:
-						probs = checkProbFuture((PathExpressionTemporal) pe, pb, p, min);
-						break;
-					case PathExpressionTemporal.P_G:
-						probs = checkProbGlobal((PathExpressionTemporal) pe, pb, p, min);
-						break;
-					}
+		if (pe instanceof PathExpressionTemporal) {
+			if (((PathExpressionTemporal) pe).hasBounds()) {
+				switch (((PathExpressionTemporal) pe).getOperator()) {
+				case PathExpressionTemporal.P_U:
+					probs = checkProbBoundedUntil((PathExpressionTemporal) pe, min);
+					break;
+				case PathExpressionTemporal.P_F:
+					probs = checkProbBoundedFuture((PathExpressionTemporal) pe, min);
+					break;
+				case PathExpressionTemporal.P_G:
+					probs = checkProbBoundedGlobal((PathExpressionTemporal) pe, min);
+					break;
+				}
+			} else {
+				switch (((PathExpressionTemporal) pe).getOperator()) {
+				case PathExpressionTemporal.P_X:
+					probs = checkProbNext((PathExpressionTemporal) pe, min);
+					break;
+				case PathExpressionTemporal.P_U:
+					probs = checkProbUntil((PathExpressionTemporal) pe, pb, p, min);
+					break;
+				case PathExpressionTemporal.P_F:
+					probs = checkProbFuture((PathExpressionTemporal) pe, pb, p, min);
+					break;
+				case PathExpressionTemporal.P_G:
+					probs = checkProbGlobal((PathExpressionTemporal) pe, pb, p, min);
+					break;
 				}
 			}
-			if (probs == null)
-				throw new PrismException("Unrecognised path operator in P operator");
-		} catch (PrismException e) {
-			if (filter != null)
-				JDD.Deref(filter);
-			throw e;
 		}
-
-		// round off probabilities
-		// probs.roundOff(getPlacesToRoundBy());
+		if (probs == null)
+			throw new PrismException("Unrecognised path operator in P operator");
 
 		// print out probabilities
 		if (verbose) {
@@ -434,141 +211,35 @@ public class NondetModelChecker implements ModelChecker
 			probs.print(mainLog);
 		}
 
-		// if a filter was provided, there is some additional output...
-		if (filter != null) {
-
-			// for non-"P=?"-type properties, print out some probs
-			if (pb != null) {
-				if (!expr.noFilterRequests()) {
-					mainLog
-							.println("\nWarning: \"{min}\", \"{max}\" only apply to \"P=?\" properties; they are ignored here.");
-				}
-				mainLog.print("\n" + (min ? "Minimum" : "Maximum")
-						+ " probabilities (non-zero only) for states satisfying " + expr.getFilter().getExpression() + ":\n");
-				probs.printFiltered(mainLog, filter);
-			}
-
-			// for "P=?"-type properties...
-			if (pb == null) {
-				// compute/print min info
-				if (expr.filterMinRequested()) {
-					minRes = probs.minOverBDD(filter);
-					mainLog.print("\nMinimum " + (min ? "minimum" : "maximum") + " probability for states satisfying "
-							+ expr.getFilter().getExpression() + ": " + minRes + "\n");
-					tmp = probs.getBDDFromInterval(minRes - termCritParam, minRes + termCritParam);
-					JDD.Ref(filter);
-					tmp = JDD.And(tmp, filter);
-					states = new StateListMTBDD(tmp, model);
-					mainLog.print("There are " + states.size() + " states with this minimum probability (+/- "
-							+ termCritParam + ")");
-					if (!verbose && states.size() > 10) {
-						mainLog
-								.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
-						states.print(mainLog, 10);
-					} else {
-						mainLog.print(":\n");
-						states.print(mainLog);
-					}
-					JDD.Deref(tmp);
-				}
-				// compute/print min info
-				if (expr.filterMaxRequested()) {
-					maxRes = probs.maxOverBDD(filter);
-					mainLog.print("\nMaximum " + (min ? "minimum" : "maximum") + " probability for states satisfying "
-							+ expr.getFilter().getExpression() + ": " + maxRes + "\n");
-					tmp = probs.getBDDFromInterval(maxRes - termCritParam, maxRes + termCritParam);
-					JDD.Ref(filter);
-					tmp = JDD.And(tmp, filter);
-					states = new StateListMTBDD(tmp, model);
-					mainLog.print("There are " + states.size() + " states with this maximum probability (+/- "
-							+ termCritParam + ")");
-					if (!verbose && states.size() > 10) {
-						mainLog
-								.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
-						states.print(mainLog, 10);
-					} else {
-						mainLog.print(":\n");
-						states.print(mainLog);
-					}
-					JDD.Deref(tmp);
-				}
-			}
-		}
-
-		// for P=? properties...
-		// if there are multiple initial states or if there is a filter
-		// satisfying
-		// multiple states with no {min}/{max}/etc., print a warning...
+		// For =? properties, just return values
 		if (pb == null) {
-			if (filter == null) {
-				if (model.getNumStartStates() > 1) {
-					mainLog
-							.print("\nWarning: There are multiple initial states; the result of model checking is for the first one: ");
-					model.getStartStates().print(mainLog, 1);
-				}
-			} else if (expr.noFilterRequests()) {
-				StateListMTBDD filterStates = new StateListMTBDD(filter, model);
-				if (filterStates.size() > 1) {
-					mainLog.print("\nWarning: The filter {" + expr.getFilter().getExpression() + "} is satisfied by "
-							+ filterStates.size() + " states.\n");
-					mainLog.print("The result of model checking is for the first of these: ");
-					filterStates.print(mainLog, 1);
-				}
-			}
+			return probs;
 		}
-
-		// compute result of property
-		// if there's a bound, get set of satisfying states
-		if (pb != null) {
+		// Otherwise, compare against bound to get set of satisfying states
+		else {
 			sol = probs.getBDDFromInterval(relOp, p);
 			// remove unreachable states from solution
 			JDD.Ref(reach);
 			sol = JDD.And(sol, reach);
+			// free vector
+			probs.clear();
+			return new StateProbsMTBDD(sol, model);
 		}
-		// if there's no bound, result will be a probability
-		else {
-			// just store empty set for sol
-			sol = JDD.Constant(0);
-			// use filter if present
-			if (filter != null) {
-				if (expr.filterMinRequested())
-					numericalRes = minRes;
-				else if (expr.filterMaxRequested())
-					numericalRes = maxRes;
-				else
-					numericalRes = probs.firstFromBDD(filter);
-			}
-			// otherwise use initial state
-			else {
-				numericalRes = probs.firstFromBDD(start);
-			}
-		}
-
-		// free vector
-		probs.clear();
-
-		// derefs
-		if (filter != null)
-			JDD.Deref(filter);
-
-		return sol;
 	}
 
 	// R operator
 
-	private JDDNode checkExpressionReward(ExpressionReward expr) throws PrismException
+	protected StateProbs checkExpressionReward(ExpressionReward expr) throws PrismException
 	{
 		Object rs; // reward struct index
 		Expression rb; // reward bound (expression)
-		double r = 0; // reward value (actual value)
+		double r = 0; // reward bound (actual value)
 		String relOp; // relational operator
 		boolean min; // are we finding min (true) or max (false) rewards
 		PathExpression pe; // path expression
 
-		JDDNode stateRewards = null, transRewards = null, filter, sol, tmp;
+		JDDNode stateRewards = null, transRewards = null, sol;
 		StateProbs rewards = null;
-		StateList states = null;
-		double minRes = 0, maxRes = 0;
 		int i;
 
 		// get info from reward operator
@@ -605,11 +276,11 @@ public class NondetModelChecker implements ModelChecker
 				mainLog.print("\nWarning: checking for reward " + relOp + " " + r
 						+ " - formula trivially satisfies all states\n");
 				JDD.Ref(reach);
-				return reach;
+				return new StateProbsMTBDD(reach, model);
 			} else if (r == 0 && relOp.equals("<")) {
 				mainLog.print("\nWarning: checking for reward " + relOp + " " + r
 						+ " - formula trivially satisfies no states\n");
-				return JDD.Constant(0);
+				return new StateProbsMTBDD(JDD.Constant(0), model);
 			}
 		}
 
@@ -624,51 +295,20 @@ public class NondetModelChecker implements ModelChecker
 			throw new PrismException("Can't use \"R=?\" for MDPs; use \"Rmin=?\" or \"Rmax=?\"");
 		}
 
-		// translate filter (if present)
-		filter = null;
-		if (expr.getFilter() != null) {
-			filter = checkExpression(expr.getFilter().getExpression());
-		}
-
-		// check if filter satisfies no states
-		if (filter != null)
-			if (filter.equals(JDD.ZERO)) {
-				// for R=? properties, this is an error
-				if (rb == null) {
-					throw new PrismException("Filter {" + expr.getFilter().getExpression() + "} in R=?[] property satisfies no states");
-				}
-				// otherwise just print a warning
-				else {
-					mainLog.println("\nWarning: Filter {" + expr.getFilter().getExpression()
-							+ "} satisfies no states and is being ignored");
-					JDD.Deref(filter);
-					filter = null;
-				}
-			}
-
 		// compute rewards
 		pe = expr.getPathExpression();
-		try {
-			if (pe instanceof PathExpressionTemporal) {
-				switch (((PathExpressionTemporal) pe).getOperator()) {
-				case PathExpressionTemporal.R_I:
-					rewards = checkRewardInst((PathExpressionTemporal) pe, stateRewards, transRewards, min);
-					break;
-				case PathExpressionTemporal.R_F:
-					rewards = checkRewardReach((PathExpressionTemporal) pe, stateRewards, transRewards, min);
-					break;
-				}
+		if (pe instanceof PathExpressionTemporal) {
+			switch (((PathExpressionTemporal) pe).getOperator()) {
+			case PathExpressionTemporal.R_I:
+				rewards = checkRewardInst((PathExpressionTemporal) pe, stateRewards, transRewards, min);
+				break;
+			case PathExpressionTemporal.R_F:
+				rewards = checkRewardReach((PathExpressionTemporal) pe, stateRewards, transRewards, min);
+				break;
 			}
-			if (rewards == null)
-				throw new PrismException("Unrecognised path operator in R operator");
-		} catch (PrismException e) {
-			if (filter != null)
-				JDD.Deref(filter);
-			throw e;
 		}
-
-		// round off rewards
-		// rewards.roundOff(getPlacesToRoundBy());
+		if (rewards == null)
+			throw new PrismException("Unrecognised path operator in R operator");
 
 		// print out rewards
 		if (verbose) {
@@ -676,129 +316,25 @@ public class NondetModelChecker implements ModelChecker
 			rewards.print(mainLog);
 		}
 
-		// if a filter was provided, there is some additional output...
-		if (filter != null) {
-
-			// for non-"R=?"-type properties, print out some rewards
-			if (rb != null) {
-				if (!expr.noFilterRequests()) {
-					mainLog
-							.println("\nWarning: \"{min}\", \"{max}\" only apply to \"R=?\" properties; they are ignored here.");
-				}
-				mainLog.print("\n" + (min ? "Minimum" : "Maximum") + " rewards (non-zero only) for states satisfying "
-						+ expr.getFilter().getExpression() + ":\n");
-				rewards.printFiltered(mainLog, filter);
-			}
-
-			// for "R=?"-type properties...
-			if (rb == null) {
-				// compute/print min info
-				if (expr.filterMinRequested()) {
-					minRes = rewards.minOverBDD(filter);
-					mainLog.print("\nMinimum " + (min ? "minimum" : "maximum") + " reward for states satisfying "
-							+ expr.getFilter().getExpression() + ": " + minRes + "\n");
-					tmp = rewards.getBDDFromInterval(minRes - termCritParam, minRes + termCritParam);
-					JDD.Ref(filter);
-					tmp = JDD.And(tmp, filter);
-					states = new StateListMTBDD(tmp, model);
-					mainLog.print("There are " + states.size() + " states with this minimum reward (+/- "
-							+ termCritParam + ")");
-					if (!verbose && states.size() > 10) {
-						mainLog
-								.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
-						states.print(mainLog, 10);
-					} else {
-						mainLog.print(":\n");
-						states.print(mainLog);
-					}
-					JDD.Deref(tmp);
-				}
-				// compute/print min info
-				if (expr.filterMaxRequested()) {
-					maxRes = rewards.maxOverBDD(filter);
-					mainLog.print("\nMaximum " + (min ? "minimum" : "maximum") + " reward for states satisfying "
-							+ expr.getFilter().getExpression() + ": " + maxRes + "\n");
-					tmp = rewards.getBDDFromInterval(maxRes - termCritParam, maxRes + termCritParam);
-					JDD.Ref(filter);
-					tmp = JDD.And(tmp, filter);
-					states = new StateListMTBDD(tmp, model);
-					mainLog.print("There are " + states.size() + " states with this maximum reward (+/- "
-							+ termCritParam + ")");
-					if (!verbose && states.size() > 10) {
-						mainLog
-								.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
-						states.print(mainLog, 10);
-					} else {
-						mainLog.print(":\n");
-						states.print(mainLog);
-					}
-					JDD.Deref(tmp);
-				}
-			}
-		}
-
-		// for R=? properties...
-		// if there are multiple initial states or if there is a filter
-		// satisfying
-		// multiple states with no {min}/{max}/etc., print a warning...
+		// For =? properties, just return values
 		if (rb == null) {
-			if (filter == null) {
-				if (model.getNumStartStates() > 1) {
-					mainLog
-							.print("\nWarning: There are multiple initial states; the result of model checking is for the first one: ");
-					model.getStartStates().print(mainLog, 1);
-				}
-			} else if (expr.noFilterRequests()) {
-				StateListMTBDD filterStates = new StateListMTBDD(filter, model);
-				if (filterStates.size() > 1) {
-					mainLog.print("\nWarning: The filter {" + expr.getFilter().getExpression() + "} is satisfied by "
-							+ filterStates.size() + " states.\n");
-					mainLog.print("The result of model checking is for the first of these: ");
-					filterStates.print(mainLog, 1);
-				}
-			}
+			return rewards;
 		}
-
-		// compute result of property
-		// if there's a bound, get set of satisfying states
-		if (rb != null) {
+		// Otherwise, compare against bound to get set of satisfying states
+		else {
 			sol = rewards.getBDDFromInterval(relOp, r);
 			// remove unreachable states from solution
 			JDD.Ref(reach);
 			sol = JDD.And(sol, reach);
+			// free vector
+			rewards.clear();
+			return new StateProbsMTBDD(sol, model);
 		}
-		// if there's no bound, result will be a reward
-		else {
-			// just store empty set for sol
-			sol = JDD.Constant(0);
-			// use filter if present
-			if (filter != null) {
-				if (expr.filterMinRequested())
-					numericalRes = minRes;
-				else if (expr.filterMaxRequested())
-					numericalRes = maxRes;
-				else
-					numericalRes = rewards.firstFromBDD(filter);
-			}
-			// otherwise use initial state
-			else {
-				numericalRes = rewards.firstFromBDD(start);
-			}
-		}
-
-		// free vector
-		rewards.clear();
-
-		// derefs
-		if (filter != null)
-			JDD.Deref(filter);
-
-		return sol;
 	}
 
 	// next
 
-	private StateProbs checkProbNext(PathExpressionTemporal pe, boolean min) throws PrismException
+	protected StateProbs checkProbNext(PathExpressionTemporal pe, boolean min) throws PrismException
 	{
 		Expression expr;
 		JDDNode b;
@@ -810,7 +346,7 @@ public class NondetModelChecker implements ModelChecker
 		expr = ((PathExpressionExpr) pe.getOperand2()).getExpression();
 
 		// model check operand first
-		b = checkExpression(expr);
+		b = checkExpressionDD(expr);
 
 		// print out some info about num states
 		// mainLog.print("\nb = " + JDD.GetNumMintermsString(b,
@@ -827,7 +363,7 @@ public class NondetModelChecker implements ModelChecker
 
 	// bounded until
 
-	private StateProbs checkProbBoundedUntil(PathExpressionTemporal pe, boolean min) throws PrismException
+	protected StateProbs checkProbBoundedUntil(PathExpressionTemporal pe, boolean min) throws PrismException
 	{
 		Expression expr1, expr2;
 		int time;
@@ -847,9 +383,9 @@ public class NondetModelChecker implements ModelChecker
 		}
 
 		// model check operands first
-		b1 = checkExpression(expr1);
+		b1 = checkExpressionDD(expr1);
 		try {
-			b2 = checkExpression(expr2);
+			b2 = checkExpressionDD(expr2);
 		} catch (PrismException e) {
 			JDD.Deref(b1);
 			throw e;
@@ -887,7 +423,7 @@ public class NondetModelChecker implements ModelChecker
 
 	// until (unbounded)
 
-	private StateProbs checkProbUntil(PathExpressionTemporal pe, Expression pb, double p, boolean min)
+	protected StateProbs checkProbUntil(PathExpressionTemporal pe, Expression pb, double p, boolean min)
 			throws PrismException
 	{
 		Expression expr1, expr2;
@@ -902,9 +438,9 @@ public class NondetModelChecker implements ModelChecker
 		expr2 = ((PathExpressionExpr) pe.getOperand2()).getExpression();
 
 		// model check operands first
-		b1 = checkExpression(expr1);
+		b1 = checkExpressionDD(expr1);
 		try {
-			b2 = checkExpression(expr2);
+			b2 = checkExpressionDD(expr2);
 		} catch (PrismException e) {
 			JDD.Deref(b1);
 			throw e;
@@ -987,7 +523,7 @@ public class NondetModelChecker implements ModelChecker
 	// bounded future (eventually)
 	// F<=k phi == true U<=k phi
 
-	private StateProbs checkProbBoundedFuture(PathExpressionTemporal pe, boolean min) throws PrismException
+	protected StateProbs checkProbBoundedFuture(PathExpressionTemporal pe, boolean min) throws PrismException
 	{
 		PathExpressionTemporal pe2;
 		pe2 = new PathExpressionTemporal(PathExpressionTemporal.P_U, new PathExpressionExpr(Expression.True()), pe
@@ -998,7 +534,7 @@ public class NondetModelChecker implements ModelChecker
 	// future (eventually)
 	// F phi == true U phi
 
-	private StateProbs checkProbFuture(PathExpressionTemporal pe, Expression pb, double p, boolean min)
+	protected StateProbs checkProbFuture(PathExpressionTemporal pe, Expression pb, double p, boolean min)
 			throws PrismException
 	{
 		PathExpressionTemporal pe2;
@@ -1012,12 +548,13 @@ public class NondetModelChecker implements ModelChecker
 	// P(G<=k phi) == 1-P(true U<=k !phi)
 	// Pmin(G<=k phi) == 1-Pmax(true U<=k !phi)
 
-	private StateProbs checkProbBoundedGlobal(PathExpressionTemporal pe, boolean min) throws PrismException
+	protected StateProbs checkProbBoundedGlobal(PathExpressionTemporal pe, boolean min) throws PrismException
 	{
 		PathExpressionTemporal pe2;
 		StateProbs probs;
 		pe2 = new PathExpressionTemporal(PathExpressionTemporal.P_U, new PathExpressionExpr(Expression.True()),
-				new PathExpressionExpr(Expression.Not(((PathExpressionExpr)pe.getOperand2()).getExpression())), pe.getLowerBound(), pe.getUpperBound());
+				new PathExpressionExpr(Expression.Not(((PathExpressionExpr) pe.getOperand2()).getExpression())), pe
+						.getLowerBound(), pe.getUpperBound());
 		probs = checkProbBoundedUntil(pe2, !min);
 		probs.subtractFromOne();
 		return probs;
@@ -1028,13 +565,13 @@ public class NondetModelChecker implements ModelChecker
 	// P(G phi) == 1-P(true U !phi)
 	// Pmin(G phi) == 1-Pmax(true U !phi)
 
-	private StateProbs checkProbGlobal(PathExpressionTemporal pe, Expression pb, double p, boolean min)
+	protected StateProbs checkProbGlobal(PathExpressionTemporal pe, Expression pb, double p, boolean min)
 			throws PrismException
 	{
 		PathExpressionTemporal pe2;
 		StateProbs probs;
 		pe2 = new PathExpressionTemporal(PathExpressionTemporal.P_U, new PathExpressionExpr(Expression.True()),
-				new PathExpressionExpr(Expression.Not(((PathExpressionExpr)pe.getOperand2()).getExpression())));
+				new PathExpressionExpr(Expression.Not(((PathExpressionExpr) pe.getOperand2()).getExpression())));
 		probs = checkProbUntil(pe2, pb, p, !min);
 		probs.subtractFromOne();
 		return probs;
@@ -1042,7 +579,7 @@ public class NondetModelChecker implements ModelChecker
 
 	// inst reward
 
-	private StateProbs checkRewardInst(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards,
+	protected StateProbs checkRewardInst(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards,
 			boolean min) throws PrismException
 	{
 		int time; // time bound
@@ -1062,7 +599,7 @@ public class NondetModelChecker implements ModelChecker
 
 	// reach reward
 
-	private StateProbs checkRewardReach(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards,
+	protected StateProbs checkRewardReach(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards,
 			boolean min) throws PrismException
 	{
 		Expression expr;
@@ -1075,7 +612,7 @@ public class NondetModelChecker implements ModelChecker
 		expr = ((PathExpressionExpr) pe.getOperand2()).getExpression();
 
 		// model check operand first
-		b = checkExpression(expr);
+		b = checkExpressionDD(expr);
 
 		// print out some info about num states
 		// mainLog.print("\nb = " + JDD.GetNumMintermsString(b,
@@ -1095,41 +632,13 @@ public class NondetModelChecker implements ModelChecker
 		return rewards;
 	}
 
-	// Check label
-
-	private JDDNode checkExpressionLabel(ExpressionLabel expr) throws PrismException
-	{
-		LabelList ll;
-		JDDNode dd;
-		int i;
-
-		// treat special cases
-		if (expr.getName().equals("deadlock")) {
-			dd = model.getFixedDeadlocks();
-			JDD.Ref(dd);
-			return dd;
-		}
-		else if (expr.getName().equals("init")) {
-			dd = start;
-			JDD.Ref(dd);
-			return dd;
-		}
-		// get expression associated with label
-		ll = propertiesFile.getCombinedLabelList();
-		i = ll.getLabelIndex(expr.getName());
-		if (i == -1)
-			throw new PrismException("Unknown label \"" + expr.getName() + "\" in property");
-		// check recursively
-		return checkExpression(ll.getLabel(i));
-	}
-
 	// -----------------------------------------------------------------------------------
 	// probability computation methods
 	// -----------------------------------------------------------------------------------
 
 	// compute probabilities for next
 
-	private StateProbs computeNextProbs(JDDNode tr, JDDNode b, boolean min)
+	protected StateProbs computeNextProbs(JDDNode tr, JDDNode b, boolean min)
 	{
 		JDDNode tmp;
 		StateProbs probs = null;
@@ -1156,7 +665,7 @@ public class NondetModelChecker implements ModelChecker
 
 	// compute probabilities for bounded until
 
-	private StateProbs computeBoundedUntilProbs(JDDNode tr, JDDNode tr01, JDDNode b1, JDDNode b2, int time, boolean min)
+	protected StateProbs computeBoundedUntilProbs(JDDNode tr, JDDNode tr01, JDDNode b1, JDDNode b2, int time, boolean min)
 			throws PrismException
 	{
 		JDDNode yes, no, maybe;
@@ -1260,7 +769,7 @@ public class NondetModelChecker implements ModelChecker
 	// note: this function doesn't need to know anything about fairness
 	// it is just told whether to compute min or max probabilities
 
-	private StateProbs computeUntilProbsQual(JDDNode tr01, JDDNode b1, JDDNode b2, boolean min)
+	protected StateProbs computeUntilProbsQual(JDDNode tr01, JDDNode b1, JDDNode b2, boolean min)
 	{
 		JDDNode yes = null, no = null, maybe;
 		StateProbs probs = null;
@@ -1363,7 +872,7 @@ public class NondetModelChecker implements ModelChecker
 	// note: this function doesn't need to know anything about fairness
 	// it is just told whether to compute min or max probabilities
 
-	private StateProbs computeUntilProbs(JDDNode tr, JDDNode tr01, JDDNode b1, JDDNode b2, boolean min)
+	protected StateProbs computeUntilProbs(JDDNode tr, JDDNode tr01, JDDNode b1, JDDNode b2, boolean min)
 			throws PrismException
 	{
 		JDDNode yes, no, maybe;
@@ -1484,9 +993,6 @@ public class NondetModelChecker implements ModelChecker
 				JDD.Deref(maybe);
 				throw e;
 			}
-
-			// round off probabilities
-			// probs.roundOff(getPlacesToRoundBy());
 		}
 
 		// derefs
@@ -1499,7 +1005,7 @@ public class NondetModelChecker implements ModelChecker
 
 	// compute rewards for inst reward
 
-	private StateProbs computeInstRewards(JDDNode tr, JDDNode sr, int time, boolean min) throws PrismException
+	protected StateProbs computeInstRewards(JDDNode tr, JDDNode sr, int time, boolean min) throws PrismException
 	{
 		JDDNode rewardsMTBDD;
 		DoubleVector rewardsDV;
@@ -1538,7 +1044,7 @@ public class NondetModelChecker implements ModelChecker
 
 	// compute rewards for reach reward
 
-	private StateProbs computeReachRewards(JDDNode tr, JDDNode tr01, JDDNode sr, JDDNode trr, JDDNode b, boolean min)
+	protected StateProbs computeReachRewards(JDDNode tr, JDDNode tr01, JDDNode sr, JDDNode trr, JDDNode b, boolean min)
 			throws PrismException
 	{
 		JDDNode inf, maybe, prob1, no;
diff --git a/prism/src/prism/Prism.java b/prism/src/prism/Prism.java
index 62749220..a2dba115 100644
--- a/prism/src/prism/Prism.java
+++ b/prism/src/prism/Prism.java
@@ -1128,80 +1128,41 @@ public class Prism implements PrismSettingsListener
 	// model checking
 	// returns result or throws an exception in case of error
 	
-	public Object modelCheck(Model model, PropertiesFile propertiesFile, Expression f) throws PrismException, PrismLangException
+	public Object modelCheck(Model model, PropertiesFile propertiesFile, Expression expr) throws PrismException, PrismLangException
 	{
 		Object res;
+		Filter filter;
 		
 		// Check that property is valid for this model type
 		// and create new model checker object
 		if (model instanceof ProbModel)
 		{
-			f.checkValid(ModulesFile.PROBABILISTIC);
+			expr.checkValid(ModulesFile.PROBABILISTIC);
 			mc = new ProbModelChecker(this, model, propertiesFile);
 		}
 		else if (model instanceof NondetModel)
 		{
-			f.checkValid(ModulesFile.NONDETERMINISTIC);
-			mc = new NondetModelChecker(mainLog, techLog, model, propertiesFile);
+			expr.checkValid(ModulesFile.NONDETERMINISTIC);
+			mc = new NondetModelChecker(this, model, propertiesFile);
 		}
 		else if (model instanceof StochModel)
 		{
-			f.checkValid(ModulesFile.STOCHASTIC);
-			mc = new StochModelChecker(mainLog, techLog, model, propertiesFile);
+			expr.checkValid(ModulesFile.STOCHASTIC);
+			mc = new StochModelChecker(this, model, propertiesFile);
 		}
 		else
 		{
 			throw new PrismException("Unknown model type");
 		}
 		
-		// configure model checker
-		mc.setEngine(getEngine());
-		if (!(model instanceof NondetModel))
-		{
-			mc.setOption("lineqmethod", getLinEqMethod());
-			mc.setOption("lineqmethodparam", getLinEqMethodParam());
-		}
-		mc.setOption("termcrit", getTermCrit());
-		mc.setOption("termcritparam", getTermCritParam());
-		mc.setOption("maxiters", getMaxIters());
-		mc.setOption("precomp", getPrecomp());
-		if (model instanceof StochModel)
-		{
-			mc.setOption("bscccomp", getBSCCComp());
-			mc.setOption("dossdetect", getDoSSDetect());
-		}
-		if (model instanceof NondetModel)
-		{
-			mc.setOption("fairness", getFairness());
-		}
-		if (getEngine() == HYBRID || getEngine() == SPARSE)
-		{
-			mc.setOption("compact", getCompact());
-		}
-		if (getEngine() == HYBRID)
-		{
-			mc.setOption("sbmaxmem", getSBMaxMem());
-			mc.setOption("numsblevels", getNumSBLevels());
-		}
-		if ((getEngine() == HYBRID) && !(model instanceof NondetModel))
-		{
-			mc.setOption("sormaxmem", getSORMaxMem());
-			mc.setOption("numsorlevels", getNumSORLevels());
-		}
-		
-		//		// print out which engine we are using
-		//		mainLog.print("\nEngine: ");
-		//		switch (engine) {
-		//		case MTBDD: mainLog.println("MTBDD"); break;
-		//		case SPARSE: mainLog.println("Sparse"); break;
-		//		case HYBRID: mainLog.println("Hybrid"); break;
-		//		default: mainLog.print("?????"); break;
-		//		}
-		
-		// do model checking
-		res = mc.check(f);
+		// Do model checking
+		if (expr instanceof ExpressionProb)
+			filter = ((ExpressionProb)expr).getFilter();
+		else
+			filter = null;
+		res = mc.check(expr, filter);
 		
-		// return result
+		// Return result
 		return res;
 	}
 	
@@ -1314,42 +1275,8 @@ public class Prism implements PrismSettingsListener
 		
 		mainLog.println("\nComputing steady-state probabilities...");
 		
-		// create/configure new model checker object
-		mc = new StochModelChecker(mainLog, techLog, model, null);
-		switch (getEngine())
-		{
-			case MTBDD: mc.setEngine(MTBDD); break;
-			case SPARSE: mc.setEngine(SPARSE); break;
-			case HYBRID: mc.setEngine(HYBRID); break;
-			default: break;
-		}
-		mc.setOption("lineqmethod", getLinEqMethod());
-		mc.setOption("lineqmethodparam", getLinEqMethodParam());
-		mc.setOption("termcrit", getTermCrit());
-		mc.setOption("termcritparam", getTermCritParam());
-		mc.setOption("maxiters", getMaxIters());
-		mc.setOption("bscccomp", getBSCCComp());
-		if (getEngine() == HYBRID || getEngine() == SPARSE)
-		{
-			mc.setOption("compact", getCompact());
-		}
-		if (getEngine() == HYBRID)
-		{
-			mc.setOption("sbmaxmem", getSBMaxMem());
-			mc.setOption("numsblevels", getNumSBLevels());
-			mc.setOption("sormaxmem", getSORMaxMem());
-			mc.setOption("numsorlevels", getNumSORLevels());
-		}
-		
-		//		// print out which engine we are using
-		//		mainLog.print("\nEngine: ");
-		//		switch (engine) {
-		//		case MTBDD: mainLog.println("MTBDD"); break;
-		//		case SPARSE: mainLog.println("Sparse"); break;
-		//		case HYBRID: mainLog.println("Hybrid"); break;
-		//		default: mainLog.print("?????"); break;
-		//		}
-		
+		// create new model checker object
+		mc = new StochModelChecker(this, model, null);
 		// do steady state calculation
 		l = System.currentTimeMillis();
 		probs = ((StochModelChecker)mc).doSteadyState();
@@ -1373,42 +1300,8 @@ public class Prism implements PrismSettingsListener
 		
 		mainLog.println("\nComputing transient probabilities (time = " + time + ")...");
 		
-		// create/configure new model checker object
-		mc = new StochModelChecker(mainLog, techLog, model, null);
-		switch (getEngine())
-		{
-			case MTBDD: mc.setEngine(MTBDD); break;
-			case SPARSE: mc.setEngine(SPARSE); break;
-			case HYBRID: mc.setEngine(HYBRID); break;
-			default: break;
-		}
-		mc.setOption("lineqmethod", getLinEqMethod()); // don't really need
-		mc.setOption("lineqmethodparam", getLinEqMethodParam()); // don't really need
-		mc.setOption("termcrit", getTermCrit());
-		mc.setOption("termcritparam", getTermCritParam());
-		mc.setOption("maxiters", getMaxIters()); // don't really need
-		mc.setOption("bscccomp", getBSCCComp()); // don't really need
-		mc.setOption("dossdetect", getDoSSDetect());
-		if (getEngine() == HYBRID || getEngine() == SPARSE)
-		{
-			mc.setOption("compact", getCompact());
-		}
-		if (getEngine() == HYBRID)
-		{
-			mc.setOption("sbmaxmem", getSBMaxMem());
-			mc.setOption("numsblevels", getNumSBLevels());
-			mc.setOption("sormaxmem", getSORMaxMem()); // don't really need
-			mc.setOption("numsorlevels", getNumSORLevels()); // don't really need
-		}
-		
-		//		// print out which engine we are using
-		//		mainLog.print("\nEngine: ");
-		//		switch (engine) {
-		//		case MTBDD: mainLog.println("MTBDD"); break;
-		//		case SPARSE: mainLog.println("Sparse"); break;
-		//		case HYBRID: mainLog.println("Hybrid"); break;
-		//		default: mainLog.print("?????"); break;
-		//		}
+		// create new model checker object
+		mc = new StochModelChecker(this, model, null);
 		
 		// do steady state calculation
 		l = System.currentTimeMillis();
diff --git a/prism/src/prism/ProbModelChecker.java b/prism/src/prism/ProbModelChecker.java
index 8d11b156..120994bd 100644
--- a/prism/src/prism/ProbModelChecker.java
+++ b/prism/src/prism/ProbModelChecker.java
@@ -26,196 +26,76 @@
 
 package prism;
 
+import java.util.Vector;
+
 import jdd.*;
-import odd.*;
 import dv.*;
 import mtbdd.*;
 import sparse.*;
 import hybrid.*;
-import parser.*;
 import parser.ast.*;
 
-// Model checker for DTMCs
-
+/*
+ * Model checker for DTMCs.
+ */
 public class ProbModelChecker extends StateModelChecker
 {
-	// options
-
-	// which engine to use
-	private int engine;
-	// method for solving linear equation systems
-	private int linEqMethod;
-	// parameter for linear equation solver methods
-	private double linEqMethodParam;
-	// termination criterion (iterative methods)
-	private int termCrit;
-	// parameter for termination criterion
-	private double termCritParam;
-	// max num iterations (iterative methods)
-	private int maxIters;
-	// flags
-	private boolean precomp; // use 0,1 precomputation algorithms?
-	// sparse bits info
-	private int SBMaxMem;
-	private int numSBLevels;
-	// hybrid sor info
-	private int SORMaxMem;
-	private int numSORLevels;
-
-	// properties file
-	private PropertiesFile propertiesFile;
-
-	// constant values
-	private Values constantValues;
-
-	// Expression2MTBDD object for translating expressions
-	private Expression2MTBDD expr2mtbdd;
-
-	// class-wide storage for any numerical result returned
-	private double numericalRes;
-
-	// model info
-	private ProbModel model;
-	private JDDNode trans;
-	private JDDNode trans01;
-	private JDDNode start;
-	private JDDNode reach;
-	private ODDNode odd;
-	private JDDVars allDDRowVars;
-	private JDDVars allDDColVars;
-	private JDDVars[] varDDRowVars;
-
-	// constructor - set some defaults
-
-	public ProbModelChecker(Prism prism, Model m, PropertiesFile pf) throws PrismException
-	{
-		super(prism, m, pf);
-		
-		// initialise
-		if (!(m instanceof ProbModel)) {
-			throw new PrismException("Wrong model type passed to ProbModelChecker.");
-		}
-		model = (ProbModel) m;
-		propertiesFile = pf;
-		trans = model.getTrans();
-		trans01 = model.getTrans01();
-		start = model.getStart();
-		reach = model.getReach();
-		odd = model.getODD();
-		allDDRowVars = model.getAllDDRowVars();
-		allDDColVars = model.getAllDDColVars();
-		varDDRowVars = model.getVarDDRowVars();
-
-		// create list of all constant values needed
-		constantValues = new Values();
-		constantValues.addValues(model.getConstantValues());
-		if (pf != null)
-			constantValues.addValues(pf.getConstantValues());
-
-		// create Expression2MTBDD object
-		expr2mtbdd = new Expression2MTBDD(mainLog, techLog, model.getVarList(), varDDRowVars, constantValues);
-		expr2mtbdd.setFilter(reach);
-
-		// set up some default options
-		// (although all should be overridden before model checking)
-		engine = Prism.HYBRID;
-		linEqMethod = Prism.JACOBI;
-		linEqMethodParam = 0.9;
-		termCrit = Prism.RELATIVE;
-		termCritParam = 1e-6;
-		maxIters = 10000;
-		precomp = true;
-		SBMaxMem = 1024;
-		numSBLevels = -1;
-		SORMaxMem = 1024;
-		numSORLevels = -1;
-	}
+	// SCC computer
+	protected SCCComputer sccComputer;
 
-	// set engine
+	// Options (in addition to those inherited from StateModelChecker):
 
-	public void setEngine(int e)
-	{
-		engine = e;
-	}
+	// Use 0,1 precomputation algorithms?
+	protected boolean precomp;
+	// Do BSCC computation?
+	protected boolean bsccComp;
 
-	// set options (generic)
+	// Constructor
 
-	public void setOption(String option, boolean b)
-	{
-		if (option.equals("precomp")) {
-			precomp = b;
-		} else if (option.equals("compact")) {
-			PrismSparse.setCompact(b);
-			PrismHybrid.setCompact(b);
-		} else {
-			mainLog.println("Warning: option \"" + option + "\" not supported by ProbModelChecker.");
-		}
-	}
-
-	public void setOption(String option, int i)
+	public ProbModelChecker(Prism prism, Model m, PropertiesFile pf) throws PrismException
 	{
-		if (option.equals("lineqmethod")) {
-			linEqMethod = i;
-			PrismMTBDD.setLinEqMethod(i);
-			PrismSparse.setLinEqMethod(i);
-			PrismHybrid.setLinEqMethod(i);
-		} else if (option.equals("termcrit")) {
-			termCrit = i;
-			PrismMTBDD.setTermCrit(i);
-			PrismSparse.setTermCrit(i);
-			PrismHybrid.setTermCrit(i);
-		} else if (option.equals("maxiters")) {
-			maxIters = i;
-			PrismMTBDD.setMaxIters(i);
-			PrismSparse.setMaxIters(i);
-			PrismHybrid.setMaxIters(i);
-		} else if (option.equals("sbmaxmem")) {
-			SBMaxMem = i;
-			PrismHybrid.setSBMaxMem(i);
-		} else if (option.equals("numsblevels")) {
-			numSBLevels = i;
-			PrismHybrid.setNumSBLevels(i);
-		} else if (option.equals("sormaxmem")) {
-			SORMaxMem = i;
-			PrismHybrid.setSORMaxMem(i);
-		} else if (option.equals("numsorlevels")) {
-			numSORLevels = i;
-			PrismHybrid.setNumSORLevels(i);
-		} else {
-			mainLog.println("Warning: option \"" + option + "\" not supported by ProbModelChecker.");
-		}
-	}
+		// Initialise
+		super(prism, m, pf);
 
-	public void setOption(String option, double d)
-	{
-		if (option.equals("lineqmethodparam")) {
-			linEqMethodParam = d;
-			PrismMTBDD.setLinEqMethodParam(d);
-			PrismSparse.setLinEqMethodParam(d);
-			PrismHybrid.setLinEqMethodParam(d);
-		} else if (option.equals("termcritparam")) {
-			termCritParam = d;
-			PrismMTBDD.setTermCritParam(d);
-			PrismSparse.setTermCritParam(d);
-			PrismHybrid.setTermCritParam(d);
-		} else {
-			mainLog.println("Warning: option \"" + option + "\" not supported by ProbModelChecker.");
+		// Create SCCComputer object
+		sccComputer = new SCCComputer(mainLog, techLog, model);
+		bsccComp = prism.getBSCCComp();
+
+		// Inherit some options from parent Prism object.
+		// Store locally and/or pass onto engines.
+		precomp = prism.getPrecomp();
+		switch (engine) {
+		case Prism.MTBDD:
+			PrismMTBDD.setLinEqMethod(prism.getLinEqMethod());
+			PrismMTBDD.setLinEqMethodParam(prism.getLinEqMethodParam());
+			PrismMTBDD.setTermCrit(prism.getTermCrit());
+			PrismMTBDD.setTermCritParam(prism.getTermCritParam());
+			PrismMTBDD.setMaxIters(prism.getMaxIters());
+			PrismMTBDD.setDoSSDetect(prism.getDoSSDetect());
+			break;
+		case Prism.SPARSE:
+			PrismSparse.setLinEqMethod(prism.getLinEqMethod());
+			PrismSparse.setLinEqMethodParam(prism.getLinEqMethodParam());
+			PrismSparse.setTermCrit(prism.getTermCrit());
+			PrismSparse.setTermCritParam(prism.getTermCritParam());
+			PrismSparse.setMaxIters(prism.getMaxIters());
+			PrismSparse.setCompact(prism.getCompact());
+			PrismSparse.setDoSSDetect(prism.getDoSSDetect());
+		case Prism.HYBRID:
+			PrismHybrid.setLinEqMethod(prism.getLinEqMethod());
+			PrismHybrid.setLinEqMethodParam(prism.getLinEqMethodParam());
+			PrismHybrid.setTermCrit(prism.getTermCrit());
+			PrismHybrid.setTermCritParam(prism.getTermCritParam());
+			PrismHybrid.setMaxIters(prism.getMaxIters());
+			PrismHybrid.setCompact(prism.getCompact());
+			PrismHybrid.setSBMaxMem(prism.getSBMaxMem());
+			PrismHybrid.setNumSBLevels(prism.getNumSBLevels());
+			PrismHybrid.setSORMaxMem(prism.getSORMaxMem());
+			PrismHybrid.setNumSORLevels(prism.getNumSORLevels());
+			PrismHybrid.setDoSSDetect(prism.getDoSSDetect());
 		}
 	}
 
-	public void setOption(String option, String s)
-	{
-		mainLog.println("Warning: option \"" + option + "\" not supported by ProbModelChecker.");
-	}
-
-	// compute number of places to round solution vector to
-	// (given termination criteria, etc.)
-
-	public int getPlacesToRoundBy()
-	{
-		return 1 - (int) (Math.log(termCritParam) / Math.log(10));
-	}
-
 	// -----------------------------------------------------------------------------------
 	// Check a property, i.e. an expression
 	// -----------------------------------------------------------------------------------
@@ -232,11 +112,11 @@ public class ProbModelChecker extends StateModelChecker
 		}
 		// R operator
 		else if (expr instanceof ExpressionReward) {
-			res = null;//checkExpressionReward((ExpressionReward) expr);
+			res = checkExpressionReward((ExpressionReward) expr);
 		}
-		// Label
-		else if (expr instanceof ExpressionLabel) {
-			res = null;//checkExpressionLabel((ExpressionLabel) expr);
+		// S operator
+		else if (expr instanceof ExpressionSS) {
+			res = checkExpressionSteadyState((ExpressionSS) expr);
 		}
 		// Otherwise, use the superclass
 		else {
@@ -255,17 +135,15 @@ public class ProbModelChecker extends StateModelChecker
 
 	// P operator
 
-	private StateProbs checkExpressionProb(ExpressionProb expr) throws PrismException
+	protected StateProbs checkExpressionProb(ExpressionProb expr) throws PrismException
 	{
 		Expression pb; // probability bound (expression)
-		double p = 0; // probability value (actual value)
+		double p = 0; // probability bound (actual value)
 		String relOp; // relational operator
 		PathExpression pe; // path expression
 
-		JDDNode filter, sol, tmp;
+		JDDNode sol;
 		StateProbs probs = null;
-		StateList states = null;
-		double minRes = 0, maxRes = 0;
 
 		// get info from prob operator
 		relOp = expr.getRelOp();
@@ -273,7 +151,7 @@ public class ProbModelChecker extends StateModelChecker
 		if (pb != null) {
 			p = pb.evaluateDouble(constantValues, null);
 			if (p < 0 || p > 1)
-				throw new PrismException("Invalid probability bound " + p + " in P[] formula");
+				throw new PrismException("Invalid probability bound " + p + " in P operator");
 		}
 
 		// check for trivial (i.e. stupid) cases
@@ -295,71 +173,40 @@ public class ProbModelChecker extends StateModelChecker
 			mainLog.print("\nWarning: \"Pmin=?\" and \"Pmax=?\" operators are identical to \"P=?\" for DTMCs\n");
 		}
 
-		// translate filter (if present)
-		filter = null;
-		if (expr.getFilter() != null) {
-			filter = checkExpressionDD(expr.getFilter().getExpression());
-		}
-
-		// check if filter satisfies no states
-		if (filter != null)
-			if (filter.equals(JDD.ZERO)) {
-				// for P=? properties, this is an error
-				if (pb == null) {
-					throw new PrismException("Filter {" + expr.getFilter().getExpression() + "} in P=?[] property satisfies no states");
-				}
-				// otherwise just print a warning
-				else {
-					mainLog.println("\nWarning: Filter {" + expr.getFilter().getExpression()
-							+ "} satisfies no states and is being ignored");
-					JDD.Deref(filter);
-					filter = null;
-				}
-			}
-
 		// compute probabilities
 		pe = expr.getPathExpression();
-		try {
-			if (pe instanceof PathExpressionTemporal) {
-				if (((PathExpressionTemporal) pe).hasBounds()) {
-					switch (((PathExpressionTemporal) pe).getOperator()) {
-					case PathExpressionTemporal.P_U:
-						probs = checkProbBoundedUntil((PathExpressionTemporal) pe);
-						break;
-					case PathExpressionTemporal.P_F:
-						probs = checkProbBoundedFuture((PathExpressionTemporal) pe);
-						break;
-					case PathExpressionTemporal.P_G:
-						probs = checkProbBoundedGlobal((PathExpressionTemporal) pe);
-						break;
-					}
-				} else {
-					switch (((PathExpressionTemporal) pe).getOperator()) {
-					case PathExpressionTemporal.P_X:
-						probs = checkProbNext((PathExpressionTemporal) pe);
-						break;
-					case PathExpressionTemporal.P_U:
-						probs = checkProbUntil((PathExpressionTemporal) pe, pb, p);
-						break;
-					case PathExpressionTemporal.P_F:
-						probs = checkProbFuture((PathExpressionTemporal) pe, pb, p);
-						break;
-					case PathExpressionTemporal.P_G:
-						probs = checkProbGlobal((PathExpressionTemporal) pe, pb, p);
-						break;
-					}
+		if (pe instanceof PathExpressionTemporal) {
+			if (((PathExpressionTemporal) pe).hasBounds()) {
+				switch (((PathExpressionTemporal) pe).getOperator()) {
+				case PathExpressionTemporal.P_U:
+					probs = checkProbBoundedUntil((PathExpressionTemporal) pe);
+					break;
+				case PathExpressionTemporal.P_F:
+					probs = checkProbBoundedFuture((PathExpressionTemporal) pe);
+					break;
+				case PathExpressionTemporal.P_G:
+					probs = checkProbBoundedGlobal((PathExpressionTemporal) pe);
+					break;
+				}
+			} else {
+				switch (((PathExpressionTemporal) pe).getOperator()) {
+				case PathExpressionTemporal.P_X:
+					probs = checkProbNext((PathExpressionTemporal) pe);
+					break;
+				case PathExpressionTemporal.P_U:
+					probs = checkProbUntil((PathExpressionTemporal) pe, pb, p);
+					break;
+				case PathExpressionTemporal.P_F:
+					probs = checkProbFuture((PathExpressionTemporal) pe, pb, p);
+					break;
+				case PathExpressionTemporal.P_G:
+					probs = checkProbGlobal((PathExpressionTemporal) pe, pb, p);
+					break;
 				}
 			}
-			if (probs == null)
-				throw new PrismException("Unrecognised path operator in P operator");
-		} catch (PrismException e) {
-			if (filter != null)
-				JDD.Deref(filter);
-			throw e;
 		}
-
-		// round off probabilities
-		// probs.roundOff(getPlacesToRoundBy());
+		if (probs == null)
+			throw new PrismException("Unrecognised path operator in P operator");
 
 		// print out probabilities
 		if (prism.getVerbose()) {
@@ -367,143 +214,34 @@ public class ProbModelChecker extends StateModelChecker
 			probs.print(mainLog);
 		}
 
-		// if a filter was provided, there is some additional output...
-		if (filter != null) {
-
-			// for non-"P=?"-type properties, print out some probs
-			if (pb != null) {
-				if (!expr.noFilterRequests()) {
-					mainLog
-							.println("\nWarning: \"{min}\", \"{max}\" only apply to \"P=?\" properties; they are ignored here.");
-				}
-				mainLog.print("\nProbabilities (non-zero only) for states satisfying " + expr.getFilter().getExpression() + ":\n");
-				probs.printFiltered(mainLog, filter);
-			}
-
-			// for "P=?"-type properties...
-			if (pb == null) {
-				// compute/print min info
-				if (expr.filterMinRequested()) {
-					minRes = probs.minOverBDD(filter);
-					mainLog.print("\nMinimum probability for states satisfying " + expr.getFilter().getExpression() + ": " + minRes
-							+ "\n");
-					tmp = probs.getBDDFromInterval(minRes - termCritParam, minRes + termCritParam);
-					JDD.Ref(filter);
-					tmp = JDD.And(tmp, filter);
-					states = new StateListMTBDD(tmp, model);
-					mainLog.print("There are " + states.size() + " states with this minimum probability (+/- "
-							+ termCritParam + ")");
-					if (!prism.getVerbose() && states.size() > 10) {
-						mainLog
-								.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
-						states.print(mainLog, 10);
-					} else {
-						mainLog.print(":\n");
-						states.print(mainLog);
-					}
-					JDD.Deref(tmp);
-				}
-				// compute/print min info
-				if (expr.filterMaxRequested()) {
-					maxRes = probs.maxOverBDD(filter);
-					mainLog.print("\nMaximum probability for states satisfying " + expr.getFilter().getExpression() + ": " + maxRes
-							+ "\n");
-					tmp = probs.getBDDFromInterval(maxRes - termCritParam, maxRes + termCritParam);
-					JDD.Ref(filter);
-					tmp = JDD.And(tmp, filter);
-					states = new StateListMTBDD(tmp, model);
-					mainLog.print("There are " + states.size() + " states with this maximum probability (+/- "
-							+ termCritParam + ")");
-					if (!prism.getVerbose() && states.size() > 10) {
-						mainLog
-								.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
-						states.print(mainLog, 10);
-					} else {
-						mainLog.print(":\n");
-						states.print(mainLog);
-					}
-					JDD.Deref(tmp);
-				}
-			}
-		}
-
-		// for P=? properties...
-		// if there are multiple initial states or if there is a filter
-		// satisfying
-		// multiple states with no {min}/{max}/etc., print a warning...
+		// For =? properties, just return values
 		if (pb == null) {
-			if (filter == null) {
-				if (model.getNumStartStates() > 1) {
-					mainLog
-							.print("\nWarning: There are multiple initial states; the result of model checking is for the first one: ");
-					model.getStartStates().print(mainLog, 1);
-				}
-			} else if (expr.noFilterRequests()) {
-				StateListMTBDD filterStates = new StateListMTBDD(filter, model);
-				if (filterStates.size() > 1) {
-					mainLog.print("\nWarning: The filter {" + expr.getFilter().getExpression() + "} is satisfied by "
-							+ filterStates.size() + " states.\n");
-					mainLog.print("The result of model checking is for the first of these: ");
-					filterStates.print(mainLog, 1);
-				}
-			}
+			return probs;
 		}
-
-		// compute result of property
-		// if there's a bound, get set of satisfying states
-		if (pb != null) {
+		// Otherwise, compare against bound to get set of satisfying states
+		else {
 			sol = probs.getBDDFromInterval(relOp, p);
 			// remove unreachable states from solution
 			JDD.Ref(reach);
 			sol = JDD.And(sol, reach);
+			// free vector
+			probs.clear();
+			return new StateProbsMTBDD(sol, model);
 		}
-		// if there's no bound, result will be a probability
-		else {
-			return probs;
-			
-			/*
-			 // just store empty set for sol
-			sol = JDD.Constant(0);
-			// use filter if present
-			if (filter != null) {
-				if (expr.filterMinRequested())
-					numericalRes = minRes;
-				else if (expr.filterMaxRequested())
-					numericalRes = maxRes;
-				else
-					numericalRes = probs.firstFromBDD(filter);
-			}
-			// otherwise use initial state
-			else {
-				numericalRes = probs.firstFromBDD(start);
-			}
-			*/
-		}
-
-		// free vector
-		probs.clear();
-
-		// derefs
-		if (filter != null)
-			JDD.Deref(filter);
-
-		return new StateProbsMTBDD(sol, model);
 	}
 
 	// R operator
 
-	private JDDNode checkExpressionReward(ExpressionReward expr) throws PrismException
+	protected StateProbs checkExpressionReward(ExpressionReward expr) throws PrismException
 	{
 		Object rs; // reward struct index
 		Expression rb; // reward bound (expression)
-		double r = 0; // reward value (actual value)
+		double r = 0; // reward bound (actual value)
 		String relOp; // relational operator
 		PathExpression pe; // path expression
 
-		JDDNode stateRewards = null, transRewards = null, filter, sol, tmp;
+		JDDNode stateRewards = null, transRewards = null, sol;
 		StateProbs rewards = null;
-		StateList states = null;
-		double minRes = 0, maxRes = 0;
 		int i;
 
 		// get info from reward operator
@@ -540,11 +278,11 @@ public class ProbModelChecker extends StateModelChecker
 				mainLog.print("\nWarning: checking for reward " + relOp + " " + r
 						+ " - formula trivially satisfies all states\n");
 				JDD.Ref(reach);
-				return reach;
+				return new StateProbsMTBDD(reach, model);
 			} else if (r == 0 && relOp.equals("<")) {
 				mainLog.print("\nWarning: checking for reward " + relOp + " " + r
 						+ " - formula trivially satisfies no states\n");
-				return JDD.Constant(0);
+				return new StateProbsMTBDD(JDD.Constant(0), model);
 			}
 		}
 
@@ -553,54 +291,26 @@ public class ProbModelChecker extends StateModelChecker
 			mainLog.print("\nWarning: \"Rmin=?\" and \"Rmax=?\" operators are identical to \"R=?\" for DTMCs\n");
 		}
 
-		// translate filter (if present)
-		filter = null;
-		if (expr.getFilter() != null) {
-			filter = checkExpressionDD(expr.getFilter().getExpression());
-		}
-
-		// check if filter satisfies no states
-		if (filter != null)
-			if (filter.equals(JDD.ZERO)) {
-				// for R=? properties, this is an error
-				if (rb == null) {
-					throw new PrismException("Filter {" + expr.getFilter().getExpression() + "} in R=?[] property satisfies no states");
-				}
-				// otherwise just print a warning
-				else {
-					mainLog.println("\nWarning: Filter {" + expr.getFilter().getExpression()
-							+ "} satisfies no states and is being ignored");
-					JDD.Deref(filter);
-					filter = null;
-				}
-			}
-
 		// compute rewards
 		pe = expr.getPathExpression();
-		try {
-			if (pe instanceof PathExpressionTemporal) {
-				switch (((PathExpressionTemporal) pe).getOperator()) {
-				case PathExpressionTemporal.R_C:
-					rewards = checkRewardCumul((PathExpressionTemporal) pe, stateRewards, transRewards);
-					break;
-				case PathExpressionTemporal.R_I:
-					rewards = checkRewardInst((PathExpressionTemporal) pe, stateRewards, transRewards);
-					break;
-				case PathExpressionTemporal.R_F:
-					rewards = checkRewardReach((PathExpressionTemporal) pe, stateRewards, transRewards);
-					break;
-				}
+		if (pe instanceof PathExpressionTemporal) {
+			switch (((PathExpressionTemporal) pe).getOperator()) {
+			case PathExpressionTemporal.R_C:
+				rewards = checkRewardCumul((PathExpressionTemporal) pe, stateRewards, transRewards);
+				break;
+			case PathExpressionTemporal.R_I:
+				rewards = checkRewardInst((PathExpressionTemporal) pe, stateRewards, transRewards);
+				break;
+			case PathExpressionTemporal.R_F:
+				rewards = checkRewardReach((PathExpressionTemporal) pe, stateRewards, transRewards);
+				break;
+			case PathExpressionTemporal.R_S:
+				rewards = checkRewardSS((PathExpressionTemporal) pe, stateRewards, transRewards);
+				break;
 			}
-			if (rewards == null)
-				throw new PrismException("Unrecognised path operator in R operator");
-		} catch (PrismException e) {
-			if (filter != null)
-				JDD.Deref(filter);
-			throw e;
 		}
-
-		// round off rewards
-		// rewards.roundOff(getPlacesToRoundBy());
+		if (rewards == null)
+			throw new PrismException("Unrecognised path operator in R operator");
 
 		// print out rewards
 		if (prism.getVerbose()) {
@@ -608,126 +318,225 @@ public class ProbModelChecker extends StateModelChecker
 			rewards.print(mainLog);
 		}
 
-		// if a filter was provided, there is some additional output...
-		if (filter != null) {
+		// For =? properties, just return values
+		if (rb == null) {
+			return rewards;
+		}
+		// Otherwise, compare against bound to get set of satisfying states
+		else {
+			sol = rewards.getBDDFromInterval(relOp, r);
+			// remove unreachable states from solution
+			JDD.Ref(reach);
+			sol = JDD.And(sol, reach);
+			// free vector
+			rewards.clear();
+			return new StateProbsMTBDD(sol, model);
+		}
+	}
 
-			// for non-"R=?"-type properties, print out some rewards
-			if (rb != null) {
-				if (!expr.noFilterRequests()) {
-					mainLog
-							.println("\nWarning: \"{min}\", \"{max}\" only apply to \"R=?\" properties; they are ignored here.");
-				}
-				mainLog.print("\nRewards (non-zero only) for states satisfying " + expr.getFilter().getExpression() + ":\n");
-				rewards.printFiltered(mainLog, filter);
-			}
+	// S operator
 
-			// for "R=?"-type properties...
-			if (rb == null) {
-				// compute/print min info
-				if (expr.filterMinRequested()) {
-					minRes = rewards.minOverBDD(filter);
-					mainLog.print("\nMinimum reward for states satisfying " + expr.getFilter().getExpression() + ": " + minRes + "\n");
-					tmp = rewards.getBDDFromInterval(minRes - termCritParam, minRes + termCritParam);
-					JDD.Ref(filter);
-					tmp = JDD.And(tmp, filter);
-					states = new StateListMTBDD(tmp, model);
-					mainLog.print("There are " + states.size() + " states with this minimum reward (+/- "
-							+ termCritParam + ")");
-					if (!prism.getVerbose() && states.size() > 10) {
-						mainLog
-								.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
-						states.print(mainLog, 10);
-					} else {
-						mainLog.print(":\n");
-						states.print(mainLog);
-					}
-					JDD.Deref(tmp);
-				}
-				// compute/print min info
-				if (expr.filterMaxRequested()) {
-					maxRes = rewards.maxOverBDD(filter);
-					mainLog.print("\nMaximum reward for states satisfying " + expr.getFilter().getExpression() + ": " + maxRes + "\n");
-					tmp = rewards.getBDDFromInterval(maxRes - termCritParam, maxRes + termCritParam);
-					JDD.Ref(filter);
-					tmp = JDD.And(tmp, filter);
-					states = new StateListMTBDD(tmp, model);
-					mainLog.print("There are " + states.size() + " states with this maximum reward (+/- "
-							+ termCritParam + ")");
-					if (!prism.getVerbose() && states.size() > 10) {
-						mainLog
-								.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
-						states.print(mainLog, 10);
-					} else {
-						mainLog.print(":\n");
-						states.print(mainLog);
-					}
-					JDD.Deref(tmp);
-				}
+	protected StateProbs checkExpressionSteadyState(ExpressionSS expr) throws PrismException
+	{
+		Expression pb; // probability bound (expression)
+		double p = 0; // probability bound (actual value)
+		String relOp; // relational operator
+
+		// bscc stuff
+		Vector vectBSCCs;
+		JDDNode notInBSCCs;
+		// mtbdd stuff
+		JDDNode b, bscc, sol, tmp;
+		// other stuff
+		StateProbs probs = null, totalProbs = null;
+		int i, n;
+		double d, probBSCCs[];
+
+		// get info from steady-state operator
+		relOp = expr.getRelOp();
+		pb = expr.getProb();
+		if (pb != null) {
+			p = pb.evaluateDouble(constantValues, null);
+			if (p < 0 || p > 1)
+				throw new PrismException("Invalid probability bound " + p + " in until formula");
+		}
+
+		// check for trivial (i.e. stupid) cases
+		if (pb != null) {
+			if ((p == 0 && relOp.equals(">=")) || (p == 1 && relOp.equals("<="))) {
+				mainLog.print("\nWarning: checking for probability " + relOp + " " + p
+						+ " - formula trivially satisfies all states\n");
+				JDD.Ref(reach);
+				return new StateProbsMTBDD(reach, model);
+			} else if ((p == 0 && relOp.equals("<")) || (p == 1 && relOp.equals(">"))) {
+				mainLog.print("\nWarning: checking for probability " + relOp + " " + p
+						+ " - formula trivially satisfies no states\n");
+				return new StateProbsMTBDD(JDD.Constant(0), model);
 			}
 		}
 
-		// for R=? properties...
-		// if there are multiple initial states or if there is a filter
-		// satisfying
-		// multiple states with no {min}/{max}/etc., print a warning...
-		if (rb == null) {
-			if (filter == null) {
-				if (model.getNumStartStates() > 1) {
-					mainLog
-							.print("\nWarning: There are multiple initial states; the result of model checking is for the first one: ");
-					model.getStartStates().print(mainLog, 1);
-				}
-			} else if (expr.noFilterRequests()) {
-				StateListMTBDD filterStates = new StateListMTBDD(filter, model);
-				if (filterStates.size() > 1) {
-					mainLog.print("\nWarning: The filter {" + expr.getFilter().getExpression() + "} is satisfied by "
-							+ filterStates.size() + " states.\n");
-					mainLog.print("The result of model checking is for the first of these: ");
-					filterStates.print(mainLog, 1);
+		// model check argument
+		b = checkExpressionDD(expr.getExpression());
+
+		// compute bottom strongly connected components (bsccs)
+		if (bsccComp) {
+			mainLog.print("\nComputing (B)SCCs...");
+			sccComputer.computeBSCCs();
+			vectBSCCs = sccComputer.getVectBSCCs();
+			notInBSCCs = sccComputer.getNotInBSCCs();
+			n = vectBSCCs.size();
+		}
+		// unless we've been told to skip it
+		else {
+			mainLog.println("\nSkipping BSCC computation...");
+			vectBSCCs = new Vector();
+			JDD.Ref(reach);
+			vectBSCCs.add(reach);
+			notInBSCCs = JDD.Constant(0);
+			n = 1;
+		}
+
+		// compute steady state for each bscc...
+		probBSCCs = new double[n];
+		for (i = 0; i < n; i++) {
+
+			mainLog.println("\nComputing steady state probabilities for BSCC " + (i + 1));
+
+			// get bscc
+			bscc = (JDDNode) vectBSCCs.elementAt(i);
+
+			// compute steady state probabilities
+			try {
+				probs = computeSteadyStateProbs(trans, bscc);
+			} catch (PrismException e) {
+				JDD.Deref(b);
+				for (i = 0; i < n; i++) {
+					JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
 				}
+				JDD.Deref(notInBSCCs);
+				throw e;
+			}
+
+			// print out probabilities
+			if (verbose) {
+				mainLog.print("\nBSCC " + (i + 1) + " steady-state probabilities: \n");
+				probs.print(mainLog);
 			}
+
+			// sum probabilities over bdd b
+			d = probs.sumOverBDD(b);
+			probBSCCs[i] = d;
+			mainLog.print("\nBSCC " + (i + 1) + " probability: " + d + "\n");
+
+			// free vector
+			probs.clear();
 		}
 
-		// compute result of property
-		// if there's a bound, get set of satisfying states
-		if (rb != null) {
-			sol = rewards.getBDDFromInterval(relOp, r);
-			// remove unreachable states from solution
-			JDD.Ref(reach);
-			sol = JDD.And(sol, reach);
+		// if every state is in a bscc, it's much easier...
+		if (notInBSCCs.equals(JDD.ZERO)) {
+
+			mainLog.println("\nAll states are in a BSCC (so no reachability probabilities computed)");
+
+			// there's more efficient ways to do this if we just create the
+			// solution bdd directly
+			// but we actually build the prob vector so it can be printed out if
+			// necessary
+			tmp = JDD.Constant(0);
+			for (i = 0; i < n; i++) {
+				bscc = (JDDNode) vectBSCCs.elementAt(i);
+				JDD.Ref(bscc);
+				tmp = JDD.Apply(JDD.PLUS, tmp, JDD.Apply(JDD.TIMES, JDD.Constant(probBSCCs[i]), bscc));
+			}
+			totalProbs = new StateProbsMTBDD(tmp, model);
 		}
-		// if there's no bound, result will be a reward
+		// otherwise we have to do more work...
 		else {
-			// just store empty set for sol
-			sol = JDD.Constant(0);
-			// use filter if present
-			if (filter != null) {
-				if (expr.filterMinRequested())
-					numericalRes = minRes;
-				else if (expr.filterMaxRequested())
-					numericalRes = maxRes;
-				else
-					numericalRes = rewards.firstFromBDD(filter);
+
+			// initialise total probabilities vector
+			switch (engine) {
+			case Prism.MTBDD:
+				totalProbs = new StateProbsMTBDD(JDD.Constant(0), model);
+				break;
+			case Prism.SPARSE:
+				totalProbs = new StateProbsDV(new DoubleVector((int) model.getNumStates()), model);
+				break;
+			case Prism.HYBRID:
+				totalProbs = new StateProbsDV(new DoubleVector((int) model.getNumStates()), model);
+				break;
 			}
-			// otherwise use initial state
-			else {
-				numericalRes = rewards.firstFromBDD(start);
+
+			// compute probabilities of reaching each bscc...
+			for (i = 0; i < n; i++) {
+
+				// skip bsccs with zero probability
+				if (probBSCCs[i] == 0.0)
+					continue;
+
+				mainLog.println("\nComputing probabilities of reaching BSCC " + (i + 1));
+
+				// get bscc
+				bscc = (JDDNode) vectBSCCs.elementAt(i);
+
+				// compute probabilities
+				try {
+					probs = computeUntilProbs(trans, trans01, notInBSCCs, bscc);
+				} catch (PrismException e) {
+					JDD.Deref(b);
+					for (i = 0; i < n; i++) {
+						JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
+					}
+					JDD.Deref(notInBSCCs);
+					totalProbs.clear();
+					throw e;
+				}
+
+				// print out probabilities
+				if (verbose) {
+					mainLog.print("\nBSCC " + (i + 1) + " reachability probabilities: \n");
+					probs.print(mainLog);
+				}
+
+				// times by bscc prob, add to total
+				probs.timesConstant(probBSCCs[i]);
+				totalProbs.add(probs);
+
+				// free vector
+				probs.clear();
 			}
 		}
 
-		// free vector
-		rewards.clear();
+		// print out probabilities
+		if (verbose) {
+			mainLog.print("\nS operator probabilities: \n");
+			totalProbs.print(mainLog);
+		}
 
 		// derefs
-		if (filter != null)
-			JDD.Deref(filter);
+		JDD.Deref(b);
+		for (i = 0; i < n; i++) {
+			JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
+		}
+		JDD.Deref(notInBSCCs);
 
-		return sol;
+		// For =? properties, just return values
+		if (pb == null) {
+			return totalProbs;
+		}
+		// Otherwise, compare against bound to get set of satisfying states
+		else {
+			sol = totalProbs.getBDDFromInterval(relOp, p);
+			// remove unreachable states from solution
+			JDD.Ref(reach);
+			sol = JDD.And(sol, reach);
+			// free vector
+			totalProbs.clear();
+			return new StateProbsMTBDD(sol, model);
+		}
 	}
 
 	// next
 
-	private StateProbs checkProbNext(PathExpressionTemporal pe) throws PrismException
+	protected StateProbs checkProbNext(PathExpressionTemporal pe) throws PrismException
 	{
 		Expression expr;
 		JDDNode b;
@@ -756,7 +565,7 @@ public class ProbModelChecker extends StateModelChecker
 
 	// bounded until
 
-	private StateProbs checkProbBoundedUntil(PathExpressionTemporal pe) throws PrismException
+	protected StateProbs checkProbBoundedUntil(PathExpressionTemporal pe) throws PrismException
 	{
 		Expression expr1, expr2;
 		int time;
@@ -816,7 +625,7 @@ public class ProbModelChecker extends StateModelChecker
 
 	// until (unbounded)
 
-	private StateProbs checkProbUntil(PathExpressionTemporal pe, Expression pb, double p) throws PrismException
+	protected StateProbs checkProbUntil(PathExpressionTemporal pe, Expression pb, double p) throws PrismException
 	{
 		Expression expr1, expr2;
 		JDDNode b1, b2;
@@ -873,7 +682,7 @@ public class ProbModelChecker extends StateModelChecker
 	// bounded future (eventually)
 	// F<=k phi == true U<=k phi
 
-	private StateProbs checkProbBoundedFuture(PathExpressionTemporal pe) throws PrismException
+	protected StateProbs checkProbBoundedFuture(PathExpressionTemporal pe) throws PrismException
 	{
 		PathExpressionTemporal pe2;
 		pe2 = new PathExpressionTemporal(PathExpressionTemporal.P_U, new PathExpressionExpr(Expression.True()), pe
@@ -884,7 +693,7 @@ public class ProbModelChecker extends StateModelChecker
 	// future (eventually)
 	// F phi == true U phi
 
-	private StateProbs checkProbFuture(PathExpressionTemporal pe, Expression pb, double p) throws PrismException
+	protected StateProbs checkProbFuture(PathExpressionTemporal pe, Expression pb, double p) throws PrismException
 	{
 		PathExpressionTemporal pe2;
 		pe2 = new PathExpressionTemporal(PathExpressionTemporal.P_U, new PathExpressionExpr(Expression.True()), pe
@@ -896,7 +705,7 @@ public class ProbModelChecker extends StateModelChecker
 	// F<=k phi == true U<=k phi
 	// P(G<=k phi) == 1-P(true U<=k !phi)
 
-	private StateProbs checkProbBoundedGlobal(PathExpressionTemporal pe) throws PrismException
+	protected StateProbs checkProbBoundedGlobal(PathExpressionTemporal pe) throws PrismException
 	{
 		PathExpressionTemporal pe2;
 		StateProbs probs;
@@ -912,7 +721,7 @@ public class ProbModelChecker extends StateModelChecker
 	// G phi == !(true U !phi)
 	// P(G phi) == 1-P(true U !phi)
 
-	private StateProbs checkProbGlobal(PathExpressionTemporal pe, Expression pb, double p) throws PrismException
+	protected StateProbs checkProbGlobal(PathExpressionTemporal pe, Expression pb, double p) throws PrismException
 	{
 		PathExpressionTemporal pe2;
 		StateProbs probs;
@@ -925,7 +734,7 @@ public class ProbModelChecker extends StateModelChecker
 
 	// cumulative reward
 
-	private StateProbs checkRewardCumul(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards)
+	protected StateProbs checkRewardCumul(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards)
 			throws PrismException
 	{
 		int time; // time
@@ -956,10 +765,10 @@ public class ProbModelChecker extends StateModelChecker
 
 	// inst reward
 
-	private StateProbs checkRewardInst(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards)
+	protected StateProbs checkRewardInst(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards)
 			throws PrismException
 	{
-		int time; // time bound
+		int time; // time
 		StateProbs rewards = null;
 
 		// get info from inst reward
@@ -976,7 +785,7 @@ public class ProbModelChecker extends StateModelChecker
 
 	// reach reward
 
-	private StateProbs checkRewardReach(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards)
+	protected StateProbs checkRewardReach(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards)
 			throws PrismException
 	{
 		Expression expr;
@@ -1009,13 +818,339 @@ public class ProbModelChecker extends StateModelChecker
 		return rewards;
 	}
 
+	// steady state reward
+
+	protected StateProbs checkRewardSS(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards)
+			throws PrismException
+	{
+		// bscc stuff
+		Vector vectBSCCs;
+		JDDNode notInBSCCs;
+		// mtbdd stuff
+		JDDNode newStateRewards, bscc, tmp;
+		// other stuff
+		StateProbs probs = null, rewards = null;
+		int i, n;
+		double d, rewBSCCs[];
+
+		// compute rewards corresponding to each state
+		JDD.Ref(trans);
+		JDD.Ref(transRewards);
+		newStateRewards = JDD.SumAbstract(JDD.Apply(JDD.TIMES, trans, transRewards), allDDColVars);
+		JDD.Ref(stateRewards);
+		newStateRewards = JDD.Apply(JDD.PLUS, newStateRewards, stateRewards);
+
+		// compute bottom strongly connected components (bsccs)
+		if (bsccComp) {
+			mainLog.print("\nComputing (B)SCCs...");
+			sccComputer.computeBSCCs();
+			vectBSCCs = sccComputer.getVectBSCCs();
+			notInBSCCs = sccComputer.getNotInBSCCs();
+			n = vectBSCCs.size();
+		}
+		// unless we've been told to skip it
+		else {
+			mainLog.println("\nSkipping BSCC computation...");
+			vectBSCCs = new Vector();
+			JDD.Ref(reach);
+			vectBSCCs.add(reach);
+			notInBSCCs = JDD.Constant(0);
+			n = 1;
+		}
+
+		// compute steady state for each bscc...
+		rewBSCCs = new double[n];
+		for (i = 0; i < n; i++) {
+
+			mainLog.println("\nComputing steady state probabilities for BSCC " + (i + 1));
+
+			// get bscc
+			bscc = (JDDNode) vectBSCCs.elementAt(i);
+
+			// compute steady state probabilities
+			try {
+				probs = computeSteadyStateProbs(trans, bscc);
+			} catch (PrismException e) {
+				JDD.Deref(newStateRewards);
+				for (i = 0; i < n; i++) {
+					JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
+				}
+				JDD.Deref(notInBSCCs);
+				throw e;
+			}
+
+			// print out probabilities
+			if (verbose) {
+				mainLog.print("\nBSCC " + (i + 1) + " steady-state probabilities: \n");
+				probs.print(mainLog);
+			}
+
+			// do weighted sum of probabilities and rewards
+			JDD.Ref(bscc);
+			JDD.Ref(newStateRewards);
+			tmp = JDD.Apply(JDD.TIMES, bscc, newStateRewards);
+			d = probs.sumOverMTBDD(tmp);
+			rewBSCCs[i] = d;
+			mainLog.print("\nBSCC " + (i + 1) + " Reward: " + d + "\n");
+			JDD.Deref(tmp);
+
+			// free vector
+			probs.clear();
+		}
+
+		// if every state is in a bscc, it's much easier...
+		if (notInBSCCs.equals(JDD.ZERO)) {
+
+			mainLog.println("\nAll states are in a BSCC (so no reachability probabilities computed)");
+
+			// build the reward vector
+			tmp = JDD.Constant(0);
+			for (i = 0; i < n; i++) {
+				bscc = (JDDNode) vectBSCCs.elementAt(i);
+				JDD.Ref(bscc);
+				tmp = JDD.Apply(JDD.PLUS, tmp, JDD.Apply(JDD.TIMES, JDD.Constant(rewBSCCs[i]), bscc));
+			}
+			rewards = new StateProbsMTBDD(tmp, model);
+		}
+		// otherwise we have to do more work...
+		else {
+
+			// initialise rewards vector
+			switch (engine) {
+			case Prism.MTBDD:
+				rewards = new StateProbsMTBDD(JDD.Constant(0), model);
+				break;
+			case Prism.SPARSE:
+				rewards = new StateProbsDV(new DoubleVector((int) model.getNumStates()), model);
+				break;
+			case Prism.HYBRID:
+				rewards = new StateProbsDV(new DoubleVector((int) model.getNumStates()), model);
+				break;
+			}
+
+			// compute probabilities of reaching each bscc...
+			for (i = 0; i < n; i++) {
+
+				// skip bsccs with zero reward
+				if (rewBSCCs[i] == 0.0)
+					continue;
+
+				mainLog.println("\nComputing probabilities of reaching BSCC " + (i + 1));
+
+				// get bscc
+				bscc = (JDDNode) vectBSCCs.elementAt(i);
+
+				// compute probabilities
+				probs = computeUntilProbs(trans, trans01, notInBSCCs, bscc);
+
+				// print out probabilities
+				if (verbose) {
+					mainLog.print("\nBSCC " + (i + 1) + " reachability probabilities: \n");
+					probs.print(mainLog);
+				}
+
+				// times by bscc reward, add to total
+				probs.timesConstant(rewBSCCs[i]);
+				rewards.add(probs);
+
+				// free vector
+				probs.clear();
+			}
+		}
+
+		// derefs
+		JDD.Deref(newStateRewards);
+		for (i = 0; i < n; i++) {
+			JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
+		}
+		JDD.Deref(notInBSCCs);
+
+		return rewards;
+	}
+
+	// -----------------------------------------------------------------------------------
+	// do steady state computation
+	// -----------------------------------------------------------------------------------
+
+	// steady state computation (from initial states)
+
+	public StateProbs doSteadyState() throws PrismException
+	{
+		// bscc stuff
+		Vector vectBSCCs;
+		JDDNode notInBSCCs;
+		// mtbdd stuff
+		JDDNode start, bscc, tmp;
+		// other stuff
+		StateProbs probs = null, solnProbs = null;
+		double d, probBSCCs[];
+		int i, n, whichBSCC, bsccCount;
+
+		// compute bottom strongly connected components (bsccs)
+		if (bsccComp) {
+			mainLog.print("\nComputing (B)SCCs...");
+			sccComputer.computeBSCCs();
+			vectBSCCs = sccComputer.getVectBSCCs();
+			notInBSCCs = sccComputer.getNotInBSCCs();
+			n = vectBSCCs.size();
+		}
+		// unless we've been told to skip it
+		else {
+			mainLog.println("\nSkipping BSCC computation...");
+			vectBSCCs = new Vector();
+			JDD.Ref(reach);
+			vectBSCCs.add(reach);
+			notInBSCCs = JDD.Constant(0);
+			n = 1;
+		}
+
+		// get initial states of model
+		start = model.getStart();
+
+		// see how many bsccs contain initial states and, if just one, which one
+		whichBSCC = -1;
+		bsccCount = 0;
+		for (i = 0; i < n; i++) {
+			bscc = (JDDNode) vectBSCCs.elementAt(i);
+			JDD.Ref(bscc);
+			JDD.Ref(start);
+			tmp = JDD.And(bscc, start);
+			if (!tmp.equals(JDD.ZERO)) {
+				bsccCount++;
+				if (bsccCount == 1)
+					whichBSCC = i;
+			}
+			JDD.Deref(tmp);
+		}
+
+		// if all initial states are in a single bscc, it's easy...
+		JDD.Ref(notInBSCCs);
+		JDD.Ref(start);
+		tmp = JDD.And(notInBSCCs, start);
+		if (tmp.equals(JDD.ZERO) && bsccCount == 1) {
+
+			JDD.Deref(tmp);
+
+			mainLog.println("\nInitial states all in one BSCC (so no reachability probabilities computed)");
+
+			// get bscc
+			bscc = (JDDNode) vectBSCCs.elementAt(whichBSCC);
+
+			// compute steady-state probabilities for the bscc
+			try {
+				solnProbs = computeSteadyStateProbs(trans, bscc);
+			} catch (PrismException e) {
+				for (i = 0; i < n; i++) {
+					JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
+				}
+				JDD.Deref(notInBSCCs);
+				throw e;
+			}
+		}
+
+		// otherwise have to consider all the bsccs
+		else {
+
+			JDD.Deref(tmp);
+
+			// initialise total probabilities vector
+			switch (engine) {
+			case Prism.MTBDD:
+				solnProbs = new StateProbsMTBDD(JDD.Constant(0), model);
+				break;
+			case Prism.SPARSE:
+				solnProbs = new StateProbsDV(new DoubleVector((int) model.getNumStates()), model);
+				break;
+			case Prism.HYBRID:
+				solnProbs = new StateProbsDV(new DoubleVector((int) model.getNumStates()), model);
+				break;
+			}
+
+			// compute prob of reaching each bscc from initial state
+			probBSCCs = new double[n];
+			for (i = 0; i < n; i++) {
+
+				mainLog.println("\nComputing probability of reaching BSCC " + (i + 1));
+
+				// get bscc
+				bscc = (JDDNode) vectBSCCs.elementAt(i);
+
+				// compute probabilities
+				try {
+					probs = computeUntilProbs(trans, trans01, notInBSCCs, bscc);
+				} catch (PrismException e) {
+					for (i = 0; i < n; i++) {
+						JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
+					}
+					JDD.Deref(notInBSCCs);
+					solnProbs.clear();
+					throw e;
+				}
+
+				// sum probabilities over bdd for initial state
+				// and then divide by number of start states
+				// (we assume an equiprobable initial probability distribution
+				// over all initial states)
+				d = probs.sumOverBDD(start);
+				d /= model.getNumStartStates();
+				probBSCCs[i] = d;
+				mainLog.print("\nBSCC " + (i + 1) + " Probability: " + d + "\n");
+
+				// free vector
+				probs.clear();
+			}
+
+			// compute steady-state for each bscc
+			for (i = 0; i < n; i++) {
+
+				mainLog.println("\nComputing steady-state probabilities for BSCC " + (i + 1));
+
+				// get bscc
+				bscc = (JDDNode) vectBSCCs.elementAt(i);
+
+				// compute steady-state probabilities for the bscc
+				try {
+					probs = computeSteadyStateProbs(trans, bscc);
+				} catch (PrismException e) {
+					for (i = 0; i < n; i++) {
+						JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
+					}
+					JDD.Deref(notInBSCCs);
+					solnProbs.clear();
+					throw e;
+				}
+
+				// print out probabilities
+				if (verbose) {
+					mainLog.print("\nBSCC " + (i + 1) + " Steady-State Probabilities: \n");
+					probs.print(mainLog);
+				}
+
+				// times by bscc reach prob, add to total
+				probs.timesConstant(probBSCCs[i]);
+				solnProbs.add(probs);
+
+				// free vector
+				probs.clear();
+			}
+		}
+
+		// derefs
+		for (i = 0; i < n; i++) {
+			JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
+		}
+		JDD.Deref(notInBSCCs);
+
+		return solnProbs;
+	}
+
 	// -----------------------------------------------------------------------------------
 	// probability computation methods
 	// -----------------------------------------------------------------------------------
 
 	// compute probabilities for next
 
-	private StateProbs computeNextProbs(JDDNode tr, JDDNode b)
+	protected StateProbs computeNextProbs(JDDNode tr, JDDNode b)
 	{
 		JDDNode tmp;
 		StateProbs probs = null;
@@ -1032,7 +1167,7 @@ public class ProbModelChecker extends StateModelChecker
 
 	// compute probabilities for bounded until
 
-	private StateProbs computeBoundedUntilProbs(JDDNode tr, JDDNode tr01, JDDNode b1, JDDNode b2, int time)
+	protected StateProbs computeBoundedUntilProbs(JDDNode tr, JDDNode tr01, JDDNode b1, JDDNode b2, int time)
 			throws PrismException
 	{
 		JDDNode yes, no, maybe;
@@ -1123,7 +1258,7 @@ public class ProbModelChecker extends StateModelChecker
 
 	// compute probabilities for until (for qualitative properties)
 
-	private StateProbs computeUntilProbsQual(JDDNode tr01, JDDNode b1, JDDNode b2, double p)
+	protected StateProbs computeUntilProbsQual(JDDNode tr01, JDDNode b1, JDDNode b2, double p)
 	{
 		JDDNode yes, no, maybe;
 		StateProbs probs = null;
@@ -1195,7 +1330,7 @@ public class ProbModelChecker extends StateModelChecker
 
 	// compute probabilities for until (general case)
 
-	private StateProbs computeUntilProbs(JDDNode tr, JDDNode tr01, JDDNode b1, JDDNode b2) throws PrismException
+	protected StateProbs computeUntilProbs(JDDNode tr, JDDNode tr01, JDDNode b1, JDDNode b2) throws PrismException
 	{
 		JDDNode yes, no, maybe;
 		JDDNode probsMTBDD;
@@ -1248,8 +1383,18 @@ public class ProbModelChecker extends StateModelChecker
 
 		// if maybe is empty, we have the probabilities already
 		if (maybe.equals(JDD.ZERO)) {
-			JDD.Ref(yes);
-			probs = new StateProbsMTBDD(yes, model);
+			// we make sure to return a vector of the appropriate type
+			// (doublevector for hybrid/sparse, mtbdd for mtbdd)
+			switch (engine) {
+			case Prism.MTBDD:
+				JDD.Ref(yes);
+				probs = new StateProbsMTBDD(yes, model);
+				break;
+			case Prism.SPARSE:
+			case Prism.HYBRID:
+				probs = new StateProbsDV(yes, model);
+				break;
+			}
 		}
 		// otherwise we compute the actual probabilities
 		else {
@@ -1291,7 +1436,7 @@ public class ProbModelChecker extends StateModelChecker
 
 	// compute cumulative rewards
 
-	private StateProbs computeCumulRewards(JDDNode tr, JDDNode tr01, JDDNode sr, JDDNode trr, int time)
+	protected StateProbs computeCumulRewards(JDDNode tr, JDDNode tr01, JDDNode sr, JDDNode trr, int time)
 			throws PrismException
 	{
 		JDDNode rewardsMTBDD;
@@ -1325,7 +1470,7 @@ public class ProbModelChecker extends StateModelChecker
 
 	// compute rewards for inst reward
 
-	private StateProbs computeInstRewards(JDDNode tr, JDDNode sr, int time) throws PrismException
+	protected StateProbs computeInstRewards(JDDNode tr, JDDNode sr, int time) throws PrismException
 	{
 		JDDNode rewardsMTBDD;
 		DoubleVector rewardsDV;
@@ -1358,7 +1503,7 @@ public class ProbModelChecker extends StateModelChecker
 
 	// compute rewards for reach reward
 
-	private StateProbs computeReachRewards(JDDNode tr, JDDNode tr01, JDDNode sr, JDDNode trr, JDDNode b)
+	protected StateProbs computeReachRewards(JDDNode tr, JDDNode tr01, JDDNode sr, JDDNode trr, JDDNode b)
 			throws PrismException
 	{
 		JDDNode inf, maybe;
@@ -1432,6 +1577,85 @@ public class ProbModelChecker extends StateModelChecker
 
 		return rewards;
 	}
+
+	// compute steady-state probabilities
+
+	// tr = the rate matrix for the whole ctmc
+	// states = the subset of reachable states (e.g. bscc) for which
+	// steady-state is to be done
+
+	protected StateProbs computeSteadyStateProbs(JDDNode tr, JDDNode subset) throws PrismException
+	{
+		JDDNode trf, init;
+		long n;
+		JDDNode probsMTBDD;
+		DoubleVector probsDV;
+		StateProbs probs = null;
+
+		// work out number of states in 'subset'
+		if (tr.equals(reach)) {
+			// avoid a call to GetNumMinterms in this simple (and common) case
+			n = model.getNumStates();
+		} else {
+			n = Math.round(JDD.GetNumMinterms(subset, allDDRowVars.n()));
+		}
+
+		// special case - there is only one state in 'subset'
+		// (in fact, we need to check for this special case because the general
+		// solution work breaks)
+		if (n == 1) {
+			// answer is trivially one in the single state
+			switch (engine) {
+			case Prism.MTBDD:
+				JDD.Ref(subset);
+				return new StateProbsMTBDD(subset, model);
+			case Prism.SPARSE:
+				return new StateProbsDV(subset, model);
+			case Prism.HYBRID:
+				return new StateProbsDV(subset, model);
+			}
+		}
+
+		// filter out unwanted states from transition matrix
+		JDD.Ref(tr);
+		JDD.Ref(subset);
+		trf = JDD.Apply(JDD.TIMES, tr, subset);
+		JDD.Ref(subset);
+		trf = JDD.Apply(JDD.TIMES, trf, JDD.PermuteVariables(subset, allDDRowVars, allDDColVars));
+
+		// compute initial solution (equiprobable)
+		JDD.Ref(subset);
+		init = JDD.Apply(JDD.DIVIDE, subset, JDD.Constant(n));
+
+		try {
+			switch (engine) {
+			case Prism.MTBDD:
+				probsMTBDD = PrismMTBDD.StochSteadyState(trf, odd, init, allDDRowVars, allDDColVars);
+				probs = new StateProbsMTBDD(probsMTBDD, model);
+				break;
+			case Prism.SPARSE:
+				probsDV = PrismSparse.StochSteadyState(trf, odd, init, allDDRowVars, allDDColVars);
+				probs = new StateProbsDV(probsDV, model);
+				break;
+			case Prism.HYBRID:
+				probsDV = PrismHybrid.StochSteadyState(trf, odd, init, allDDRowVars, allDDColVars);
+				probs = new StateProbsDV(probsDV, model);
+				break;
+			default:
+				throw new PrismException("Engine does not support this numerical method");
+			}
+		} catch (PrismException e) {
+			JDD.Deref(trf);
+			JDD.Deref(init);
+			throw e;
+		}
+
+		// derefs
+		JDD.Deref(trf);
+		JDD.Deref(init);
+
+		return probs;
+	}
 }
 
 // ------------------------------------------------------------------------------
diff --git a/prism/src/prism/SCCComputer.java b/prism/src/prism/SCCComputer.java
index 246aba7e..d9b3d72d 100644
--- a/prism/src/prism/SCCComputer.java
+++ b/prism/src/prism/SCCComputer.java
@@ -39,7 +39,7 @@ public class SCCComputer
 	private PrismLog techLog;
 
 	// model info
-	private StochModel model;
+	private Model model;
 	private JDDNode trans01;
 	private JDDNode reach;
 	private JDDVars allDDRowVars;
@@ -59,10 +59,7 @@ public class SCCComputer
 		// initialise
 		mainLog = log1;
 		techLog = log2;
-		if (!(m instanceof StochModel)) {	
-			throw new PrismException("Error: Wrong model type passed to SCCComputer.");
-		}
-		model = (StochModel)m;
+		model = m;
 		trans01 = model.getTrans01();
 		reach = model.getReach();
 		allDDRowVars = model.getAllDDRowVars();
diff --git a/prism/src/prism/StateModelChecker.java b/prism/src/prism/StateModelChecker.java
index f7494033..bb676199 100644
--- a/prism/src/prism/StateModelChecker.java
+++ b/prism/src/prism/StateModelChecker.java
@@ -28,9 +28,6 @@ package prism;
 
 import jdd.*;
 import odd.*;
-import mtbdd.*;
-import sparse.*;
-import hybrid.*;
 import parser.*;
 import parser.ast.*;
 
@@ -43,26 +40,14 @@ public class StateModelChecker implements ModelChecker
 	protected PrismLog mainLog;
 	protected PrismLog techLog;
 
-	// options
-
-	// which engine to use
-	protected int engine;
-	// parameter for termination criterion
-	protected double termCritParam;
-	// flags
-	protected boolean precomp; // use 0,1 precomputation algorithms?
-
-	// properties file
+	// Properties file
 	protected PropertiesFile propertiesFile;
 
-	// constant values
+	// Constant values
 	protected Values constantValues;
 
-	// class-wide storage for any numerical result returned
-	protected double numericalRes;
-
-	// model info
-	protected ProbModel model;
+	// Model info
+	protected Model model;
 	protected VarList varList;
 	protected JDDNode trans;
 	protected JDDNode trans01;
@@ -73,19 +58,29 @@ public class StateModelChecker implements ModelChecker
 	protected JDDVars allDDColVars;
 	protected JDDVars[] varDDRowVars;
 
-	// constructor - set some defaults
+	// Options:
+	
+	// Which engine to use
+	protected int engine;
+	// Parameter for termination criterion
+	protected double termCritParam;
+	// Verbose mode?
+	protected boolean verbose;
+
+	// Constructor
 
 	public StateModelChecker(Prism prism, Model m, PropertiesFile pf) throws PrismException
 	{
-		// initialise
+		// Initialise
 		this.prism = prism;
 		mainLog = prism.getMainLog();
 		techLog = prism.getTechLog();
-		if (!(m instanceof ProbModel)) {
-			throw new PrismException("Wrong model type passed to ProbModelChecker.");
-		}
-		model = (ProbModel) m;
+		model = m;
 		propertiesFile = pf;
+		constantValues = new Values();
+		constantValues.addValues(model.getConstantValues());
+		if (pf != null)
+			constantValues.addValues(pf.getConstantValues());
 		varList = model.getVarList();
 		trans = model.getTrans();
 		trans01 = model.getTrans01();
@@ -96,86 +91,11 @@ public class StateModelChecker implements ModelChecker
 		allDDColVars = model.getAllDDColVars();
 		varDDRowVars = model.getVarDDRowVars();
 
-		// create list of all constant values needed
-		constantValues = new Values();
-		constantValues.addValues(model.getConstantValues());
-		if (pf != null)
-			constantValues.addValues(pf.getConstantValues());
-
-		// set up some default options
-		// (although all should be overridden before model checking)
-		engine = Prism.HYBRID;
-		termCritParam = 1e-6;
-		precomp = true;
-	}
-
-	// set engine
-
-	public void setEngine(int e)
-	{
-		engine = e;
-	}
-
-	// set options (generic)
-
-	public void setOption(String option, boolean b)
-	{
-		if (option.equals("precomp")) {
-			precomp = b;
-		} else if (option.equals("compact")) {
-			PrismSparse.setCompact(b);
-			PrismHybrid.setCompact(b);
-		} else {
-			mainLog.println("Warning: option \"" + option + "\" not supported by ProbModelChecker.");
-		}
-	}
-
-	public void setOption(String option, int i)
-	{
-		if (option.equals("lineqmethod")) {
-			PrismMTBDD.setLinEqMethod(i);
-			PrismSparse.setLinEqMethod(i);
-			PrismHybrid.setLinEqMethod(i);
-		} else if (option.equals("termcrit")) {
-			PrismMTBDD.setTermCrit(i);
-			PrismSparse.setTermCrit(i);
-			PrismHybrid.setTermCrit(i);
-		} else if (option.equals("maxiters")) {
-			PrismMTBDD.setMaxIters(i);
-			PrismSparse.setMaxIters(i);
-			PrismHybrid.setMaxIters(i);
-		} else if (option.equals("sbmaxmem")) {
-			PrismHybrid.setSBMaxMem(i);
-		} else if (option.equals("numsblevels")) {
-			PrismHybrid.setNumSBLevels(i);
-		} else if (option.equals("sormaxmem")) {
-			PrismHybrid.setSORMaxMem(i);
-		} else if (option.equals("numsorlevels")) {
-			PrismHybrid.setNumSORLevels(i);
-		} else {
-			mainLog.println("Warning: option \"" + option + "\" not supported by ProbModelChecker.");
-		}
-	}
-
-	public void setOption(String option, double d)
-	{
-		if (option.equals("lineqmethodparam")) {
-			PrismMTBDD.setLinEqMethodParam(d);
-			PrismSparse.setLinEqMethodParam(d);
-			PrismHybrid.setLinEqMethodParam(d);
-		} else if (option.equals("termcritparam")) {
-			termCritParam = d;
-			PrismMTBDD.setTermCritParam(d);
-			PrismSparse.setTermCritParam(d);
-			PrismHybrid.setTermCritParam(d);
-		} else {
-			mainLog.println("Warning: option \"" + option + "\" not supported by ProbModelChecker.");
-		}
-	}
-
-	public void setOption(String option, String s)
-	{
-		mainLog.println("Warning: option \"" + option + "\" not supported by ProbModelChecker.");
+		// Inherit some options from parent Prism object
+		// Store locally and/or pass onto engines
+		engine = prism.getEngine();
+		termCritParam = prism.getTermCritParam();
+		verbose = prism.getVerbose();
 	}
 
 	// -----------------------------------------------------------------------------------
@@ -183,32 +103,46 @@ public class StateModelChecker implements ModelChecker
 	// -----------------------------------------------------------------------------------
 
 	public Object check(Expression expr) throws PrismException
+	{
+		return check(expr, null);
+	}
+
+	public Object check(Expression expr, Filter filter) throws PrismException
 	{
 		long timer = 0;
-		JDDNode dd;
+		JDDNode ddFilter, tmp;
 		StateList states;
 		boolean b;
-		StateProbs res;
-		Object singleRes;
+		StateProbs vals;
+		Object res;
+		double minRes = 0.0, maxRes = 0.0;
+
+		// Translate filter (creating default one first if none present)
+		if (filter == null) {
+			filter = new Filter(new ExpressionLabel("init"));
+			if (model.getNumStartStates() > 1) {
+				mainLog
+						.print("\nWarning: There are multiple initial states; the result of model checking is for the first one: ");
+				model.getStartStates().print(mainLog, 1);
+			}
+		}
+		ddFilter = checkExpressionDD(filter.getExpression());
+		if (ddFilter.equals(JDD.ZERO)) {
+			throw new PrismException("Filter " + filter + " satisfies no states");
+		}
 
-		// start timer
+		// Do model checking and store result
 		timer = System.currentTimeMillis();
-
-		// do model checking and store result
-		res = checkExpression(expr);
-		
-		// stop timer
+		vals = checkExpression(expr);
 		timer = System.currentTimeMillis() - timer;
-
-		// print out model checking time
 		mainLog.println("\nTime for model checking: " + timer / 1000.0 + " seconds.");
 
-		// output and result of model checking depend on return type
+		// Process results of model checking using filter
 		if (expr.getType() == Expression.BOOLEAN) {
 
-			states = new StateListMTBDD(((StateProbsMTBDD)res).getJDDNode(), model);
-			JDD.Ref(((StateProbsMTBDD)res).getJDDNode());
-			
+			states = new StateListMTBDD(((StateProbsMTBDD) vals).getJDDNode(), model);
+			JDD.Ref(((StateProbsMTBDD) vals).getJDDNode());
+
 			// print out number of satisfying states
 			mainLog.print("\nNumber of satisfying states: ");
 			mainLog.print(states.size());
@@ -224,7 +158,7 @@ public class StateModelChecker implements ModelChecker
 			mainLog.print("\n");
 
 			// if "verbose", print out satisfying states
-			if (prism.getVerbose()) {
+			if (verbose) {
 				mainLog.print("\nSatisfying states:");
 				if (states.size() > 0) {
 					mainLog.print("\n");
@@ -236,23 +170,70 @@ public class StateModelChecker implements ModelChecker
 
 			// result is true if all states satisfy, false otherwise
 			b = (states.size() == model.getNumStates());
-			singleRes = b ? new Boolean(true) : new Boolean(false);
+			res = b ? new Boolean(true) : new Boolean(false);
 			states.clear();
 
 			// print result
 			mainLog.print("\nResult: " + b + " (property " + (b ? "" : "not ") + "satisfied in all states)\n");
 		} else {
-			singleRes = new Double(res.firstFromBDD(start));
-
+			if (filter.minRequested()) {
+				minRes = vals.minOverBDD(ddFilter);
+				// TODO probability -> value?
+				// (what about min/max vals for MDPs?)
+				// (see also poss use of getresultname below)
+				mainLog.print("\nMinimum probability for states satisfying " + filter.getExpression() + ": " + minRes
+						+ "\n");
+				tmp = vals.getBDDFromInterval(minRes - termCritParam, minRes + termCritParam);
+				JDD.Ref(ddFilter);
+				tmp = JDD.And(tmp, ddFilter);
+				states = new StateListMTBDD(tmp, model);
+				mainLog.print("There are " + states.size() + " states with this minimum probability (+/- "
+						+ termCritParam + ")");
+				if (!verbose && states.size() > 10) {
+					mainLog.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
+					states.print(mainLog, 10);
+				} else {
+					mainLog.print(":\n");
+					states.print(mainLog);
+				}
+				JDD.Deref(tmp);
+			}
+			if (filter.maxRequested()) {
+				maxRes = vals.maxOverBDD(ddFilter);
+				mainLog.print("\nMaximum probability for states satisfying " + filter.getExpression() + ": " + maxRes
+						+ "\n");
+				tmp = vals.getBDDFromInterval(maxRes - termCritParam, maxRes + termCritParam);
+				JDD.Ref(ddFilter);
+				tmp = JDD.And(tmp, ddFilter);
+				states = new StateListMTBDD(tmp, model);
+				mainLog.print("There are " + states.size() + " states with this maximum probability (+/- "
+						+ termCritParam + ")");
+				if (!verbose && states.size() > 10) {
+					mainLog.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
+					states.print(mainLog, 10);
+				} else {
+					mainLog.print(":\n");
+					states.print(mainLog);
+				}
+				JDD.Deref(tmp);
+			}
+			if (filter.minRequested())
+				res = minRes;
+			else if (filter.maxRequested())
+				res = maxRes;
+			else
+				res = vals.firstFromBDD(ddFilter);
 			// print result
-			mainLog.print("\nResult: " + singleRes + "\n");
+			mainLog.print("\nResult: " + res + "\n");
+			// TODO: use expr.getresultname?
 		}
 
-		// finished with memory
-		res.clear();
+		// clean up
+		JDD.Deref(ddFilter);
+		vals.clear();
 
 		// return result
-		return singleRes;
+		return res;
 	}
 
 	// Check expression (recursive)
@@ -263,50 +244,49 @@ public class StateModelChecker implements ModelChecker
 
 		// If-then-else
 		if (expr instanceof ExpressionITE) {
-			res = checkExpressionITE((ExpressionITE)expr);
+			res = checkExpressionITE((ExpressionITE) expr);
 		}
 		// Binary ops
 		else if (expr instanceof ExpressionBinaryOp) {
-			res = checkExpressionBinaryOp((ExpressionBinaryOp)expr);
+			res = checkExpressionBinaryOp((ExpressionBinaryOp) expr);
 		}
 		// Unary ops
 		else if (expr instanceof ExpressionUnaryOp) {
-			res = checkExpressionUnaryOp((ExpressionUnaryOp)expr);
+			res = checkExpressionUnaryOp((ExpressionUnaryOp) expr);
 		}
 		// Functions
 		else if (expr instanceof ExpressionFunc) {
-			res = checkExpressionFunc((ExpressionFunc)expr);
+			res = checkExpressionFunc((ExpressionFunc) expr);
 		}
 		// Identifiers
 		else if (expr instanceof ExpressionIdent) {
 			// Should never happen
-			throw new PrismException("Unknown identifier \"" + ((ExpressionIdent)expr).getName() + "\"");
+			throw new PrismException("Unknown identifier \"" + ((ExpressionIdent) expr).getName() + "\"");
 		}
 		// Literals
 		else if (expr instanceof ExpressionLiteral) {
-			res = checkExpressionLiteral((ExpressionLiteral)expr);
+			res = checkExpressionLiteral((ExpressionLiteral) expr);
 		}
 		// Constants
 		else if (expr instanceof ExpressionConstant) {
-			res = checkExpressionConstant((ExpressionConstant)expr);
+			res = checkExpressionConstant((ExpressionConstant) expr);
 		}
 		// Formulas
 		else if (expr instanceof ExpressionFormula) {
 			// Should never happen
-			throw new PrismException("Unexpanded formula \"" + ((ExpressionFormula)expr).getName() + "\"");
+			throw new PrismException("Unexpanded formula \"" + ((ExpressionFormula) expr).getName() + "\"");
 		}
 		// Variables
 		else if (expr instanceof ExpressionVar) {
-			res = checkExpressionVar((ExpressionVar)expr);
+			res = checkExpressionVar((ExpressionVar) expr);
 		}
 		// Labels
 		else if (expr instanceof ExpressionLabel) {
 			res = checkExpressionLabel((ExpressionLabel) expr);
-		}
-		else {
+		} else {
 			throw new PrismException("Couldn't check " + expr.getClass());
 		}
-		
+
 		// Filter out non-reachable states from solution (TODO: not for dv?)
 		res.filter(reach);
 
@@ -315,80 +295,95 @@ public class StateModelChecker implements ModelChecker
 
 	public JDDNode checkExpressionDD(Expression expr) throws PrismException
 	{
-		return ((StateProbsMTBDD)checkExpression(expr)).getJDDNode();
+		return ((StateProbsMTBDD) checkExpression(expr)).getJDDNode();
 	}
-	
+
 	// -----------------------------------------------------------------------------------
 	// Check method for each operator
 	// -----------------------------------------------------------------------------------
 
 	// Check an 'if-then-else'
-	
+
 	private StateProbs checkExpressionITE(ExpressionITE expr) throws PrismException
 	{
 		JDDNode dd, dd1, dd2, dd3;
-		
+
 		dd1 = checkExpressionDD(expr.getOperand1());
 		dd2 = checkExpressionDD(expr.getOperand2());
 		dd3 = checkExpressionDD(expr.getOperand3());
 		dd = JDD.ITE(dd1, dd2, dd3);
-		
+
 		return new StateProbsMTBDD(dd, model);
 	}
-	
+
 	// Check a binary operator
-	
+
 	private StateProbs checkExpressionBinaryOp(ExpressionBinaryOp expr) throws PrismException
 	{
 		JDDNode dd, tmp1, tmp2;
 		int op = expr.getOperator();
-		
+
 		// Optimisations are possible for relational operators
 		// (note dubious use of knowledge that op IDs are consecutive)
 		if (op >= ExpressionBinaryOp.EQ && op <= ExpressionBinaryOp.LE) {
 			return checkExpressionRelOp(op, expr.getOperand1(), expr.getOperand2());
 		}
-		
+
 		// Check operands
 		tmp1 = checkExpressionDD(expr.getOperand1());
 		try {
 			tmp2 = checkExpressionDD(expr.getOperand2());
-		}
-		catch (PrismException e) {
+		} catch (PrismException e) {
 			JDD.Deref(tmp1);
 			throw e;
 		}
-		
+
 		// Apply operation
 		switch (op) {
-		case ExpressionBinaryOp.IMPLIES: dd = JDD.Or(JDD.Not(tmp1), tmp2); break;
-		case ExpressionBinaryOp.OR: dd = JDD.Or(tmp1, tmp2); break;
-		case ExpressionBinaryOp.AND: dd = JDD.And(tmp1, tmp2); break;
-		case ExpressionBinaryOp.PLUS: dd = JDD.Apply(JDD.PLUS, tmp1, tmp2); break;
-		case ExpressionBinaryOp.MINUS: dd = JDD.Apply(JDD.MINUS, tmp1, tmp2); break;
-		case ExpressionBinaryOp.TIMES: dd = JDD.Apply(JDD.TIMES, tmp1, tmp2); break;
-		case ExpressionBinaryOp.DIVIDE: dd = JDD.Apply(JDD.DIVIDE, tmp1, tmp2); break;
-		default: throw new PrismException("Unknown binary operator");
+		case ExpressionBinaryOp.IMPLIES:
+			dd = JDD.Or(JDD.Not(tmp1), tmp2);
+			break;
+		case ExpressionBinaryOp.OR:
+			dd = JDD.Or(tmp1, tmp2);
+			break;
+		case ExpressionBinaryOp.AND:
+			dd = JDD.And(tmp1, tmp2);
+			break;
+		case ExpressionBinaryOp.PLUS:
+			dd = JDD.Apply(JDD.PLUS, tmp1, tmp2);
+			break;
+		case ExpressionBinaryOp.MINUS:
+			dd = JDD.Apply(JDD.MINUS, tmp1, tmp2);
+			break;
+		case ExpressionBinaryOp.TIMES:
+			dd = JDD.Apply(JDD.TIMES, tmp1, tmp2);
+			break;
+		case ExpressionBinaryOp.DIVIDE:
+			dd = JDD.Apply(JDD.DIVIDE, tmp1, tmp2);
+			break;
+		default:
+			throw new PrismException("Unknown binary operator");
 		}
 
 		return new StateProbsMTBDD(dd, model);
 	}
-	
+
 	// Check a relational operator (=, !=, >, >=, < <=)
-	
+
 	private StateProbs checkExpressionRelOp(int op, Expression expr1, Expression expr2) throws PrismException
 	{
 		JDDNode dd, tmp1, tmp2;
 		String s;
 
-		// check for some easy (and common) special cases before resorting to the general case
-		
+		// check for some easy (and common) special cases before resorting to
+		// the general case
+
 		// var relop int
-		if (expr1 instanceof ExpressionVar && expr2.isConstant() && expr2.getType()==Expression.INT) {
+		if (expr1 instanceof ExpressionVar && expr2.isConstant() && expr2.getType() == Expression.INT) {
 			ExpressionVar e1;
 			Expression e2;
 			int i, j, l, h, v;
-			e1 = (ExpressionVar)expr1;
+			e1 = (ExpressionVar) expr1;
 			e2 = expr2;
 			// get var's index
 			s = e1.getName();
@@ -404,35 +399,42 @@ public class StateModelChecker implements ModelChecker
 			i = e2.evaluateInt(constantValues, null);
 			switch (op) {
 			case ExpressionBinaryOp.EQ:
-				if (i>=l && i <= h) dd = JDD.SetVectorElement(dd, varDDRowVars[v], i-l, 1);
+				if (i >= l && i <= h)
+					dd = JDD.SetVectorElement(dd, varDDRowVars[v], i - l, 1);
 				break;
 			case ExpressionBinaryOp.NE:
-				if (i>=l && i <= h) dd = JDD.SetVectorElement(dd, varDDRowVars[v], i-l, 1);
+				if (i >= l && i <= h)
+					dd = JDD.SetVectorElement(dd, varDDRowVars[v], i - l, 1);
 				dd = JDD.Not(dd);
 				break;
 			case ExpressionBinaryOp.GT:
-				for (j = i+1; j <= h; j++) dd = JDD.SetVectorElement(dd, varDDRowVars[v], j-l, 1);
+				for (j = i + 1; j <= h; j++)
+					dd = JDD.SetVectorElement(dd, varDDRowVars[v], j - l, 1);
 				break;
 			case ExpressionBinaryOp.GE:
-				for (j = i; j <= h; j++) dd = JDD.SetVectorElement(dd, varDDRowVars[v], j-l, 1);
+				for (j = i; j <= h; j++)
+					dd = JDD.SetVectorElement(dd, varDDRowVars[v], j - l, 1);
 				break;
 			case ExpressionBinaryOp.LT:
-				for (j = i-1; j >= l; j--) dd = JDD.SetVectorElement(dd, varDDRowVars[v], j-l, 1);
+				for (j = i - 1; j >= l; j--)
+					dd = JDD.SetVectorElement(dd, varDDRowVars[v], j - l, 1);
 				break;
 			case ExpressionBinaryOp.LE:
-				for (j = i; j >= l; j--) dd = JDD.SetVectorElement(dd, varDDRowVars[v], j-l, 1);
+				for (j = i; j >= l; j--)
+					dd = JDD.SetVectorElement(dd, varDDRowVars[v], j - l, 1);
 				break;
-			default: throw new PrismException("Unknown relational operator");
+			default:
+				throw new PrismException("Unknown relational operator");
 			}
 			return new StateProbsMTBDD(dd, model);
 		}
 		// int relop var
-		else if (expr1.isConstant() && expr1.getType()==Expression.INT && expr2 instanceof ExpressionVar) {
+		else if (expr1.isConstant() && expr1.getType() == Expression.INT && expr2 instanceof ExpressionVar) {
 			Expression e1;
 			ExpressionVar e2;
 			int i, j, l, h, v;
 			e1 = expr1;
-			e2 = (ExpressionVar)expr2;
+			e2 = (ExpressionVar) expr2;
 			// get var's index
 			s = e2.getName();
 			v = varList.getIndex(s);
@@ -447,25 +449,32 @@ public class StateModelChecker implements ModelChecker
 			i = e1.evaluateInt(constantValues, null);
 			switch (op) {
 			case ExpressionBinaryOp.EQ:
-				if (i>=l && i <= h) dd = JDD.SetVectorElement(dd, varDDRowVars[v], i-l, 1);
+				if (i >= l && i <= h)
+					dd = JDD.SetVectorElement(dd, varDDRowVars[v], i - l, 1);
 				break;
 			case ExpressionBinaryOp.NE:
-				if (i>=l && i <= h) dd = JDD.SetVectorElement(dd, varDDRowVars[v], i-l, 1);
+				if (i >= l && i <= h)
+					dd = JDD.SetVectorElement(dd, varDDRowVars[v], i - l, 1);
 				dd = JDD.Not(dd);
 				break;
 			case ExpressionBinaryOp.GT:
-				for (j = i-1; j >= l; j--) dd = JDD.SetVectorElement(dd, varDDRowVars[v], j-l, 1);
+				for (j = i - 1; j >= l; j--)
+					dd = JDD.SetVectorElement(dd, varDDRowVars[v], j - l, 1);
 				break;
 			case ExpressionBinaryOp.GE:
-				for (j = i; j >= l; j--) dd = JDD.SetVectorElement(dd, varDDRowVars[v], j-l, 1);
+				for (j = i; j >= l; j--)
+					dd = JDD.SetVectorElement(dd, varDDRowVars[v], j - l, 1);
 				break;
 			case ExpressionBinaryOp.LT:
-				for (j = i+1; j <= h; j++) dd = JDD.SetVectorElement(dd, varDDRowVars[v], j-l, 1);
+				for (j = i + 1; j <= h; j++)
+					dd = JDD.SetVectorElement(dd, varDDRowVars[v], j - l, 1);
 				break;
 			case ExpressionBinaryOp.LE:
-				for (j = i; j <= h; j++) dd = JDD.SetVectorElement(dd, varDDRowVars[v], j-l, 1);
+				for (j = i; j <= h; j++)
+					dd = JDD.SetVectorElement(dd, varDDRowVars[v], j - l, 1);
 				break;
-			default: throw new PrismException("Unknown relational operator");
+			default:
+				throw new PrismException("Unknown relational operator");
 			}
 			return new StateProbsMTBDD(dd, model);
 		}
@@ -492,185 +501,214 @@ public class StateModelChecker implements ModelChecker
 		case ExpressionBinaryOp.LE:
 			dd = JDD.Apply(JDD.LESSTHANEQUALS, tmp1, tmp2);
 			break;
-		default: throw new PrismException("Unknown relational operator");
+		default:
+			throw new PrismException("Unknown relational operator");
 		}
-		
+
 		return new StateProbsMTBDD(dd, model);
 	}
-	
+
 	// Check a unary operator
-	
+
 	private StateProbs checkExpressionUnaryOp(ExpressionUnaryOp expr) throws PrismException
 	{
 		JDDNode dd, tmp;
 		int op = expr.getOperator();
-		
+
 		// Check operand
 		tmp = checkExpressionDD(expr.getOperand());
-		
+
 		// Apply operation
 		switch (op) {
-		case ExpressionUnaryOp.NOT: dd = JDD.Not(tmp); break;
-		case ExpressionUnaryOp.MINUS: dd = JDD.Apply(JDD.MINUS, JDD.Constant(0), tmp); break;
-		case ExpressionUnaryOp.PARENTH: dd = tmp; break;
-		default: throw new PrismException("Unknown unary operator");
+		case ExpressionUnaryOp.NOT:
+			dd = JDD.Not(tmp);
+			break;
+		case ExpressionUnaryOp.MINUS:
+			dd = JDD.Apply(JDD.MINUS, JDD.Constant(0), tmp);
+			break;
+		case ExpressionUnaryOp.PARENTH:
+			dd = tmp;
+			break;
+		default:
+			throw new PrismException("Unknown unary operator");
 		}
 
 		return new StateProbsMTBDD(dd, model);
 	}
-	
+
 	// Check a 'function'
-	
+
 	private StateProbs checkExpressionFunc(ExpressionFunc expr) throws PrismException
 	{
 		switch (expr.getNameCode()) {
-			case ExpressionFunc.MIN: return checkExpressionFuncMin(expr);
-			case ExpressionFunc.MAX: return checkExpressionFuncMax(expr);
-			case ExpressionFunc.FLOOR: return checkExpressionFuncFloor(expr);
-			case ExpressionFunc.CEIL: return checkExpressionFuncCeil(expr);
-			case ExpressionFunc.POW: return checkExpressionFuncPow(expr);
-			case ExpressionFunc.MOD: return checkExpressionFuncMod(expr);
-			case ExpressionFunc.LOG: return checkExpressionFuncLog(expr);
-			default: throw new PrismException("Unrecognised function \"" + expr.getName() + "\"");
+		case ExpressionFunc.MIN:
+			return checkExpressionFuncMin(expr);
+		case ExpressionFunc.MAX:
+			return checkExpressionFuncMax(expr);
+		case ExpressionFunc.FLOOR:
+			return checkExpressionFuncFloor(expr);
+		case ExpressionFunc.CEIL:
+			return checkExpressionFuncCeil(expr);
+		case ExpressionFunc.POW:
+			return checkExpressionFuncPow(expr);
+		case ExpressionFunc.MOD:
+			return checkExpressionFuncMod(expr);
+		case ExpressionFunc.LOG:
+			return checkExpressionFuncLog(expr);
+		default:
+			throw new PrismException("Unrecognised function \"" + expr.getName() + "\"");
 		}
 	}
-	
+
 	private StateProbs checkExpressionFuncMin(ExpressionFunc expr) throws PrismException
 	{
 		int i, n;
 		JDDNode dd, tmp;
-		
+
 		dd = checkExpressionDD(expr.getOperand(0));
 		n = expr.getNumOperands();
 		for (i = 1; i < n; i++) {
 			try {
 				tmp = checkExpressionDD(expr.getOperand(i));
 				dd = JDD.Apply(JDD.MIN, dd, tmp);
-			}
-			catch (PrismException e) {
+			} catch (PrismException e) {
 				JDD.Deref(dd);
 				throw e;
 			}
 		}
-		
+
 		return new StateProbsMTBDD(dd, model);
 	}
-	
+
 	private StateProbs checkExpressionFuncMax(ExpressionFunc expr) throws PrismException
 	{
 		int i, n;
 		JDDNode dd, tmp;
-		
+
 		dd = checkExpressionDD(expr.getOperand(0));
 		n = expr.getNumOperands();
 		for (i = 1; i < n; i++) {
 			try {
 				tmp = checkExpressionDD(expr.getOperand(i));
 				dd = JDD.Apply(JDD.MAX, dd, tmp);
-			}
-			catch (PrismException e) {
+			} catch (PrismException e) {
 				JDD.Deref(dd);
 				throw e;
 			}
 		}
-		
+
 		return new StateProbsMTBDD(dd, model);
 	}
-	
+
 	private StateProbs checkExpressionFuncFloor(ExpressionFunc expr) throws PrismException
 	{
 		JDDNode dd;
-		
+
 		dd = checkExpressionDD(expr.getOperand(0));
 		dd = JDD.MonadicApply(JDD.FLOOR, dd);
-		
+
 		return new StateProbsMTBDD(dd, model);
 	}
-	
+
 	private StateProbs checkExpressionFuncCeil(ExpressionFunc expr) throws PrismException
 	{
 		JDDNode dd;
-		
+
 		dd = checkExpressionDD(expr.getOperand(0));
 		dd = JDD.MonadicApply(JDD.CEIL, dd);
-		
+
 		return new StateProbsMTBDD(dd, model);
 	}
-	
+
 	private StateProbs checkExpressionFuncPow(ExpressionFunc expr) throws PrismException
 	{
 		JDDNode dd1, dd2, dd;
-		
+
 		dd1 = checkExpressionDD(expr.getOperand(0));
 		dd2 = checkExpressionDD(expr.getOperand(1));
 		dd = JDD.Apply(JDD.POW, dd1, dd2);
-		
+
 		return new StateProbsMTBDD(dd, model);
 	}
-	
+
 	private StateProbs checkExpressionFuncMod(ExpressionFunc expr) throws PrismException
 	{
 		JDDNode dd1, dd2, dd;
-		
+
 		dd1 = checkExpressionDD(expr.getOperand(0));
 		dd2 = checkExpressionDD(expr.getOperand(1));
 		dd = JDD.Apply(JDD.MOD, dd1, dd2);
-		
+
 		return new StateProbsMTBDD(dd, model);
 	}
 
 	private StateProbs checkExpressionFuncLog(ExpressionFunc expr) throws PrismException
 	{
 		JDDNode dd1, dd2, dd;
-		
+
 		dd1 = checkExpressionDD(expr.getOperand(0));
 		dd2 = checkExpressionDD(expr.getOperand(1));
 		dd = JDD.Apply(JDD.LOGXY, dd1, dd2);
-		
+
 		return new StateProbsMTBDD(dd, model);
 	}
 
 	// Check a literal
-	
+
 	private StateProbs checkExpressionLiteral(ExpressionLiteral expr) throws PrismException
 	{
 		JDDNode dd;
 		switch (expr.getType()) {
-			case Expression.BOOLEAN: dd = JDD.Constant(expr.evaluateBoolean(null, null) ? 1.0 : 0.0); break;
-			case Expression.INT: dd = JDD.Constant(expr.evaluateInt(null, null)); break;
-			case Expression.DOUBLE: dd = JDD.Constant(expr.evaluateDouble(null, null)); break;
-			default: throw new PrismException("Unknown literal type");
+		case Expression.BOOLEAN:
+			dd = JDD.Constant(expr.evaluateBoolean(null, null) ? 1.0 : 0.0);
+			break;
+		case Expression.INT:
+			dd = JDD.Constant(expr.evaluateInt(null, null));
+			break;
+		case Expression.DOUBLE:
+			dd = JDD.Constant(expr.evaluateDouble(null, null));
+			break;
+		default:
+			throw new PrismException("Unknown literal type");
 		}
 		return new StateProbsMTBDD(dd, model);
 	}
-	
+
 	// Check a constant
-	
+
 	private StateProbs checkExpressionConstant(ExpressionConstant expr) throws PrismException
 	{
 		int i;
 		JDDNode dd;
-		
+
 		i = constantValues.getIndexOf(expr.getName());
-		if (i == -1) throw new PrismException("Couldn't evaluate constant \"" + expr.getName() + "\"");
+		if (i == -1)
+			throw new PrismException("Couldn't evaluate constant \"" + expr.getName() + "\"");
 		switch (constantValues.getType(i)) {
-			case Expression.INT: dd = JDD.Constant(constantValues.getIntValue(i)); break;
-			case Expression.DOUBLE: dd = JDD.Constant(constantValues.getDoubleValue(i)); break;
-			case Expression.BOOLEAN: dd = JDD.Constant(constantValues.getBooleanValue(i) ? 1.0 : 0.0); break;
-			default: throw new PrismException("Unknown type for constant \"" + expr.getName() + "\"");
+		case Expression.INT:
+			dd = JDD.Constant(constantValues.getIntValue(i));
+			break;
+		case Expression.DOUBLE:
+			dd = JDD.Constant(constantValues.getDoubleValue(i));
+			break;
+		case Expression.BOOLEAN:
+			dd = JDD.Constant(constantValues.getBooleanValue(i) ? 1.0 : 0.0);
+			break;
+		default:
+			throw new PrismException("Unknown type for constant \"" + expr.getName() + "\"");
 		}
-		
+
 		return new StateProbsMTBDD(dd, model);
 	}
-	
+
 	// Check a variable reference
-	
+
 	private StateProbs checkExpressionVar(ExpressionVar expr) throws PrismException
 	{
 		String s;
 		int v, l, h, i;
 		JDDNode dd;
-		
+
 		s = expr.getName();
 		// get the variable's index
 		v = varList.getIndex(s);
@@ -683,12 +721,12 @@ public class StateModelChecker implements ModelChecker
 		// create dd
 		dd = JDD.Constant(0);
 		for (i = l; i <= h; i++) {
-			dd = JDD.SetVectorElement(dd, varDDRowVars[v], i-l, i);
+			dd = JDD.SetVectorElement(dd, varDDRowVars[v], i - l, i);
 		}
-				
+
 		return new StateProbsMTBDD(dd, model);
 	}
-	
+
 	// Check label
 
 	private StateProbs checkExpressionLabel(ExpressionLabel expr) throws PrismException
@@ -701,12 +739,10 @@ public class StateModelChecker implements ModelChecker
 		if (expr.getName().equals("deadlock")) {
 			dd = model.getFixedDeadlocks();
 			JDD.Ref(dd);
-		}
-		else if (expr.getName().equals("init")) {
+		} else if (expr.getName().equals("init")) {
 			dd = start;
 			JDD.Ref(dd);
-		}
-		else {
+		} else {
 			// get expression associated with label
 			ll = propertiesFile.getCombinedLabelList();
 			i = ll.getLabelIndex(expr.getName());
@@ -715,7 +751,7 @@ public class StateModelChecker implements ModelChecker
 			// check recursively
 			dd = checkExpressionDD(ll.getLabel(i));
 		}
-		
+
 		return new StateProbsMTBDD(dd, model);
 	}
 }
diff --git a/prism/src/prism/StochModelChecker.java b/prism/src/prism/StochModelChecker.java
index 792274b4..b393ca4c 100644
--- a/prism/src/prism/StochModelChecker.java
+++ b/prism/src/prism/StochModelChecker.java
@@ -26,1309 +26,161 @@
 
 package prism;
 
-import java.util.Vector;
-
 import jdd.*;
-import odd.*;
 import dv.*;
 import mtbdd.*;
 import sparse.*;
 import hybrid.*;
-import parser.*;
 import parser.ast.*;
 
-// Model checker for CTMCs
-
-public class StochModelChecker implements ModelChecker
+/*
+ * Model checker for CTMCs.
+ * 
+ * Much of StochModelChecker's functionality is inherited from the parent
+ * class ProbModelChecker (for DTMCs). Main differences are: 
+ *  - bounded until: time bounds are doubles and computation different
+ *  - next/unbounded until: prob computation uses embedded Markov chain
+ *  - cumulative/instantaneous reward: times are doubles, computation different
+ *  - reach rewards: ...
+ *  - doTransient
+ *  
+ *  TODO: finish this doc
+ */
+public class StochModelChecker extends ProbModelChecker
 {
-	// logs
-	private PrismLog mainLog;
-	private PrismLog techLog;
-
-	// options
-
-	// which engine to use
-	private int engine;
-	// method for solving linear equation systems
-	private int linEqMethod;
-	// parameter for linear equation solver methods
-	private double linEqMethodParam;
-	// termination criterion (iterative methods)
-	private int termCrit;
-	// parameter for termination criterion
-	private double termCritParam;
-	// max num iterations (iterative methods)
-	private int maxIters;
-	// flags
-	private boolean verbose; // verbose output?
-	private boolean precomp; // use 0,1 precomputation algorithms?
-	private boolean bsccComp; // do BSCC computation?
-	// sparse bits info
-	private int SBMaxMem;
-	private int numSBLevels;
-	// hybrid sor info
-	private int SORMaxMem;
-	private int numSORLevels;
-
-	// properties file
-	private PropertiesFile propertiesFile;
-
-	// constant values
-	private Values constantValues;
-
-	// Expression2MTBDD object for translating expressions
-	private Expression2MTBDD expr2mtbdd;
-
-	// scc computer
-	private SCCComputer sccComputer;
-
-	// class-wide storage for probability in the initial state
-	private double numericalRes;
-
-	// model info
-	private StochModel model;
-	private JDDNode trans;
-	private JDDNode trans01;
-	private JDDNode start;
-	private JDDNode reach;
-	private ODDNode odd;
-	private JDDVars allDDRowVars;
-	private JDDVars allDDColVars;
-	private JDDVars[] varDDRowVars;
-	private long numStates;
-
-	// constructor - set some defaults
-
-	public StochModelChecker(PrismLog log1, PrismLog log2, Model m, PropertiesFile pf) throws PrismException
-	{
-		// initialise
-		mainLog = log1;
-		techLog = log2;
-		if (!(m instanceof StochModel)) {
-			throw new PrismException("Wrong model type passed to StochModelChecker.");
-		}
-		model = (StochModel) m;
-		propertiesFile = pf;
-		trans = model.getTrans();
-		trans01 = model.getTrans01();
-		start = model.getStart();
-		reach = model.getReach();
-		odd = model.getODD();
-		allDDRowVars = model.getAllDDRowVars();
-		allDDColVars = model.getAllDDColVars();
-		varDDRowVars = model.getVarDDRowVars();
-		numStates = model.getNumStates();
-
-		// create list of all constant values needed
-		constantValues = new Values();
-		constantValues.addValues(model.getConstantValues());
-		if (pf != null)
-			constantValues.addValues(pf.getConstantValues());
-
-		// create Expression2MTBDD object
-		expr2mtbdd = new Expression2MTBDD(mainLog, techLog, model.getVarList(), varDDRowVars, constantValues);
-		expr2mtbdd.setFilter(reach);
-
-		// create SCCComputer object
-		sccComputer = new SCCComputer(log1, log2, model);
-
-		// set up some default options
-		// (although all should be overridden before model checking)
-		engine = Prism.HYBRID;
-		linEqMethod = Prism.JACOBI;
-		linEqMethodParam = 0.9;
-		termCrit = Prism.RELATIVE;
-		termCritParam = 1e-6;
-		maxIters = 10000;
-		verbose = false;
-		precomp = true;
-		bsccComp = true;
-		SBMaxMem = 1024;
-		numSBLevels = -1;
-		SORMaxMem = 1024;
-		numSORLevels = -1;
-	}
-
-	// set engine
-
-	public void setEngine(int e)
-	{
-		engine = e;
-	}
-
-	// set options (generic)
-
-	public void setOption(String option, boolean b)
-	{
-		if (option.equals("verbose")) {
-			verbose = b;
-		} else if (option.equals("precomp")) {
-			precomp = b;
-		} else if (option.equals("bscccomp")) {
-			bsccComp = b;
-		} else if (option.equals("dossdetect")) {
-			PrismMTBDD.setDoSSDetect(b);
-			PrismSparse.setDoSSDetect(b);
-			PrismHybrid.setDoSSDetect(b);
-		} else if (option.equals("compact")) {
-			PrismHybrid.setCompact(b);
-			PrismSparse.setCompact(b);
-		} else {
-			mainLog.println("Warning: option \"" + option + "\" not supported by StochModelChecker.");
-		}
-	}
-
-	public void setOption(String option, int i)
-	{
-		if (option.equals("lineqmethod")) {
-			linEqMethod = i;
-			PrismMTBDD.setLinEqMethod(i);
-			PrismSparse.setLinEqMethod(i);
-			PrismHybrid.setLinEqMethod(i);
-		} else if (option.equals("termcrit")) {
-			termCrit = i;
-			PrismMTBDD.setTermCrit(i);
-			PrismSparse.setTermCrit(i);
-			PrismHybrid.setTermCrit(i);
-		} else if (option.equals("maxiters")) {
-			maxIters = i;
-			PrismMTBDD.setMaxIters(i);
-			PrismSparse.setMaxIters(i);
-			PrismHybrid.setMaxIters(i);
-		} else if (option.equals("sbmaxmem")) {
-			SBMaxMem = i;
-			PrismHybrid.setSBMaxMem(i);
-		} else if (option.equals("numsblevels")) {
-			numSBLevels = i;
-			PrismHybrid.setNumSBLevels(i);
-		} else if (option.equals("sormaxmem")) {
-			SORMaxMem = i;
-			PrismHybrid.setSORMaxMem(i);
-		} else if (option.equals("numsorlevels")) {
-			numSORLevels = i;
-			PrismHybrid.setNumSORLevels(i);
-		} else {
-			mainLog.println("Warning: option \"" + option + "\" not supported by StochModelChecker.");
-		}
-	}
-
-	public void setOption(String option, double d)
-	{
-		if (option.equals("lineqmethodparam")) {
-			linEqMethodParam = d;
-			PrismMTBDD.setLinEqMethodParam(d);
-			PrismSparse.setLinEqMethodParam(d);
-			PrismHybrid.setLinEqMethodParam(d);
-		} else if (option.equals("termcritparam")) {
-			termCritParam = d;
-			PrismMTBDD.setTermCritParam(d);
-			PrismSparse.setTermCritParam(d);
-			PrismHybrid.setTermCritParam(d);
-		} else {
-			mainLog.println("Warning: option \"" + option + "\" not supported by StochModelChecker.");
-		}
-	}
-
-	public void setOption(String option, String s)
-	{
-		mainLog.println("Warning: option \"" + option + "\" not supported by StochModelChecker.");
-	}
-
-	// compute number of places to round solution vector to
-	// (given termination criteria, etc.)
-
-	public int getPlacesToRoundBy()
-	{
-		return 1 - (int) (Math.log(termCritParam) / Math.log(10));
-	}
-
-	// -----------------------------------------------------------------------------------
-	// Check a property, i.e. an expression
-	// -----------------------------------------------------------------------------------
-
-	public Object check(Expression expr) throws PrismException
-	{
-		long timer = 0;
-		JDDNode dd;
-		StateList states;
-		boolean b;
-		Object res;
-
-		// start timer
-		timer = System.currentTimeMillis();
-
-		// do model checking and store result
-		dd = checkExpression(expr);
-		states = new StateListMTBDD(dd, model);
-
-		// stop timer
-		timer = System.currentTimeMillis() - timer;
-
-		// print out model checking time
-		mainLog.println("\nTime for model checking: " + timer / 1000.0 + " seconds.");
-
-		// output and result of model checking depend on return type
-		if (expr.getType() == Expression.BOOLEAN) {
-
-			// print out number of satisfying states
-			mainLog.print("\nNumber of satisfying states: ");
-			mainLog.print(states.size());
-			if (states.size() == numStates) {
-				mainLog.print(" (all)");
-			} else if (states.includes(model.getStart())) {
-				mainLog.print((model.getNumStartStates() == 1) ? " (including initial state)"
-						: " (including all initial states)");
-			} else {
-				mainLog.print((model.getNumStartStates() == 1) ? " (initial state not satisfied)"
-						: " (initial states not all satisfied)");
-			}
-			mainLog.print("\n");
-
-			// if "verbose", print out satisfying states
-			if (verbose) {
-				mainLog.print("\nSatisfying states:");
-				if (states.size() > 0) {
-					mainLog.print("\n");
-					states.print(mainLog);
-				} else {
-					mainLog.print(" (none)\n");
-				}
-			}
-
-			// result is true if all states satisfy, false otherwise
-			b = (states.size() == numStates);
-			res = b ? new Boolean(true) : new Boolean(false);
-
-			// print result
-			mainLog.print("\nResult: " + b + " (property " + (b ? "" : "not ") + "satisfied in all states)\n");
-		} else {
-			// result is value stored earlier
-			res = new Double(numericalRes);
-
-			// print result
-			mainLog.print("\nResult: " + res + "\n");
-		}
+	// Constructor
 
-		// finished with memory
-		states.clear();
-
-		// return result
-		return res;
-	}
-
-	// Check expression (recursive)
-
-	public JDDNode checkExpression(Expression expr) throws PrismException
+	public StochModelChecker(Prism prism, Model m, PropertiesFile pf) throws PrismException
 	{
-		JDDNode res;
-
-		// P operator
-		if (expr instanceof ExpressionProb) {
-			res = checkExpressionProb((ExpressionProb) expr);
-		}
-		// S operator
-		else if (expr instanceof ExpressionSS) {
-			res = checkExpressionSteadyState((ExpressionSS) expr);
-		}
-		// R operator
-		else if (expr instanceof ExpressionReward) {
-			res = checkExpressionReward((ExpressionReward) expr);
-		}
-		// Label
-		else if (expr instanceof ExpressionLabel) {
-			res = checkExpressionLabel((ExpressionLabel) expr);
-		}
-		// Otherwise, must be an ordinary expression
-		else {
-			res = expr2mtbdd.translateExpression(expr);
-		}
-
-		// Filter out non-reachable states from solution
-		JDD.Ref(reach);
-		res = JDD.Apply(JDD.TIMES, res, reach);
-
-		return res;
+		// Initialise
+		super(prism, m, pf);
 	}
 
 	// -----------------------------------------------------------------------------------
 	// Check method for each operator
 	// -----------------------------------------------------------------------------------
 
-	// P operator
+	// bounded until
 
-	private JDDNode checkExpressionProb(ExpressionProb expr) throws PrismException
+	protected StateProbs checkProbBoundedUntil(PathExpressionTemporal pe) throws PrismException
 	{
-		Expression pb; // probability bound (expression)
-		double p = 0; // probability value (actual value)
-		String relOp; // relational operator
-		PathExpression pe; // path expression
-
-		JDDNode filter, sol, tmp;
-		StateProbs probs = null;
-		StateList states = null;
-		double minRes = 0, maxRes = 0;
-
-		// get info from prob operator
-		relOp = expr.getRelOp();
-		pb = expr.getProb();
-		if (pb != null) {
-			p = pb.evaluateDouble(constantValues, null);
-			if (p < 0 || p > 1)
-				throw new PrismException("Invalid probability bound " + p + " in P[] formula");
-		}
-
-		// check for trivial (i.e. stupid) cases
-		if (pb != null) {
-			if ((p == 0 && relOp.equals(">=")) || (p == 1 && relOp.equals("<="))) {
-				mainLog.print("\nWarning: checking for probability " + relOp + " " + p
-						+ " - formula trivially satisfies all states\n");
-				JDD.Ref(reach);
-				return reach;
-			} else if ((p == 0 && relOp.equals("<")) || (p == 1 && relOp.equals(">"))) {
-				mainLog.print("\nWarning: checking for probability " + relOp + " " + p
-						+ " - formula trivially satisfies no states\n");
-				return JDD.Constant(0);
-			}
-		}
-
-		// print a warning if Pmin/Pmax used
-		if (relOp.equals("min=") || relOp.equals("max=")) {
-			mainLog.print("\nWarning: \"Pmin=?\" and \"Pmax=?\" operators are identical to \"P=?\" for DTMCs\n");
-		}
-
-		// translate filter (if present)
-		filter = null;
-		if (expr.getFilter() != null) {
-			filter = checkExpression(expr.getFilter().getExpression());
-		}
-
-		// check if filter satisfies no states
-		if (filter != null)
-			if (filter.equals(JDD.ZERO)) {
-				// for P=? properties, this is an error
-				if (pb == null) {
-					throw new PrismException("Filter {" + expr.getFilter().getExpression() + "} in P=?[] property satisfies no states");
-				}
-				// otherwise just print a warning
-				else {
-					mainLog.println("\nWarning: Filter {" + expr.getFilter().getExpression()
-							+ "} satisfies no states and is being ignored");
-					JDD.Deref(filter);
-					filter = null;
-				}
-			}
-
-		// compute probabilities
-		pe = expr.getPathExpression();
-		try {
-			if (pe instanceof PathExpressionTemporal) {
-				if (((PathExpressionTemporal) pe).hasBounds()) {
-					switch (((PathExpressionTemporal) pe).getOperator()) {
-					case PathExpressionTemporal.P_U:
-						probs = checkProbBoundedUntil((PathExpressionTemporal) pe);
-						break;
-					case PathExpressionTemporal.P_F:
-						probs = checkProbBoundedFuture((PathExpressionTemporal) pe);
-						break;
-					case PathExpressionTemporal.P_G:
-						probs = checkProbBoundedGlobal((PathExpressionTemporal) pe);
-						break;
-					}
-				} else {
-					switch (((PathExpressionTemporal) pe).getOperator()) {
-					case PathExpressionTemporal.P_X:
-						probs = checkProbNext((PathExpressionTemporal) pe);
-						break;
-					case PathExpressionTemporal.P_U:
-						probs = checkProbUntil((PathExpressionTemporal) pe, pb, p);
-						break;
-					case PathExpressionTemporal.P_F:
-						probs = checkProbFuture((PathExpressionTemporal) pe, pb, p);
-						break;
-					case PathExpressionTemporal.P_G:
-						probs = checkProbGlobal((PathExpressionTemporal) pe, pb, p);
-						break;
-					}
-				}
-			}
-			if (probs == null)
-				throw new PrismException("Unrecognised path operator in P operator");
-		} catch (PrismException e) {
-			if (filter != null)
-				JDD.Deref(filter);
-			throw e;
-		}
-
-		// round off probabilities
-		// probs.roundOff(getPlacesToRoundBy());
-
-		// print out probabilities
-		if (verbose) {
-			mainLog.print("\nProbabilities (non-zero only) for all states:\n");
-			probs.print(mainLog);
-		}
-
-		// if a filter was provided, there is some additional output...
-		if (filter != null) {
+		Expression expr1, expr2;
+		double lTime, uTime; // time bounds
+		Expression expr;
+		JDDNode b1, b2, tmp;
+		StateProbs tmpProbs = null, probs = null;
 
-			// for non-"P=?"-type properties, print out some probs
-			if (pb != null) {
-				if (!expr.noFilterRequests()) {
-					mainLog
-							.println("\nWarning: \"{min}\", \"{max}\" only apply to \"P=?\" properties; they are ignored here.");
-				}
-				mainLog.print("\nProbabilities (non-zero only) for states satisfying " + expr.getFilter().getExpression() + ":\n");
-				probs.printFiltered(mainLog, filter);
-			}
+		// get operands
+		if (!(pe.getOperand1() instanceof PathExpressionExpr) || !(pe.getOperand2() instanceof PathExpressionExpr))
+			throw new PrismException("Invalid path formula");
+		expr1 = ((PathExpressionExpr) pe.getOperand1()).getExpression();
+		expr2 = ((PathExpressionExpr) pe.getOperand2()).getExpression();
 
-			// for "P=?"-type properties...
-			if (pb == null) {
-				// compute/print min info
-				if (expr.filterMinRequested()) {
-					minRes = probs.minOverBDD(filter);
-					mainLog.print("\nMinimum probability for states satisfying " + expr.getFilter().getExpression() + ": " + minRes
-							+ "\n");
-					tmp = probs.getBDDFromInterval(minRes - termCritParam, minRes + termCritParam);
-					JDD.Ref(filter);
-					tmp = JDD.And(tmp, filter);
-					states = new StateListMTBDD(tmp, model);
-					mainLog.print("There are " + states.size() + " states with this minimum probability (+/- "
-							+ termCritParam + ")");
-					if (!verbose && states.size() > 10) {
-						mainLog
-								.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
-						states.print(mainLog, 10);
-					} else {
-						mainLog.print(":\n");
-						states.print(mainLog);
-					}
-					JDD.Deref(tmp);
-				}
-				// compute/print min info
-				if (expr.filterMaxRequested()) {
-					maxRes = probs.maxOverBDD(filter);
-					mainLog.print("\nMaximum probability for states satisfying " + expr.getFilter().getExpression() + ": " + maxRes
-							+ "\n");
-					tmp = probs.getBDDFromInterval(maxRes - termCritParam, maxRes + termCritParam);
-					JDD.Ref(filter);
-					tmp = JDD.And(tmp, filter);
-					states = new StateListMTBDD(tmp, model);
-					mainLog.print("There are " + states.size() + " states with this maximum probability (+/- "
-							+ termCritParam + ")");
-					if (!verbose && states.size() > 10) {
-						mainLog
-								.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
-						states.print(mainLog, 10);
-					} else {
-						mainLog.print(":\n");
-						states.print(mainLog);
-					}
-					JDD.Deref(tmp);
-				}
-			}
-		}
+		// get info from bounded until
 
-		// for P=? properties...
-		// if there are multiple initial states or if there is a filter
-		// satisfying
-		// multiple states with no {min}/{max}/etc., print a warning...
-		if (pb == null) {
-			if (filter == null) {
-				if (model.getNumStartStates() > 1) {
-					mainLog
-							.print("\nWarning: There are multiple initial states; the result of model checking is for the first one: ");
-					model.getStartStates().print(mainLog, 1);
-				}
-			} else if (expr.noFilterRequests()) {
-				StateListMTBDD filterStates = new StateListMTBDD(filter, model);
-				if (filterStates.size() > 1) {
-					mainLog.print("\nWarning: The filter {" + expr.getFilter().getExpression() + "} is satisfied by "
-							+ filterStates.size() + " states.\n");
-					mainLog.print("The result of model checking is for the first of these: ");
-					filterStates.print(mainLog, 1);
-				}
+		// lower bound is 0 if not specified
+		// (i.e. if until is of form U<=t)
+		expr = pe.getLowerBound();
+		if (expr != null) {
+			lTime = expr.evaluateDouble(constantValues, null);
+			if (lTime < 0) {
+				throw new PrismException("Invalid lower bound " + lTime + " in time-bounded until formula");
 			}
+		} else {
+			lTime = 0;
 		}
-
-		// compute result of property
-		// if there's a bound, get set of satisfying states
-		if (pb != null) {
-			sol = probs.getBDDFromInterval(relOp, p);
-			// remove unreachable states from solution
-			JDD.Ref(reach);
-			sol = JDD.And(sol, reach);
-		}
-		// if there's no bound, result will be a probability
-		else {
-			// just store empty set for sol
-			sol = JDD.Constant(0);
-			// use filter if present
-			if (filter != null) {
-				if (expr.filterMinRequested())
-					numericalRes = minRes;
-				else if (expr.filterMaxRequested())
-					numericalRes = maxRes;
-				else
-					numericalRes = probs.firstFromBDD(filter);
-			}
-			// otherwise use initial state
-			else {
-				numericalRes = probs.firstFromBDD(start);
+		// upper bound is -1 if not specified
+		// (i.e. if until is of form U>=t)
+		expr = pe.getUpperBound();
+		if (expr != null) {
+			uTime = expr.evaluateDouble(constantValues, null);
+			if (uTime < 0) {
+				throw new PrismException("Invalid upper bound " + uTime + " in time-bounded until formula");
 			}
-		}
-
-		// free vector
-		probs.clear();
-
-		// derefs
-		if (filter != null)
-			JDD.Deref(filter);
-
-		return sol;
-	}
-
-	// R operator
-
-	private JDDNode checkExpressionReward(ExpressionReward expr) throws PrismException
-	{
-		Object rs; // reward struct index
-		Expression rb; // reward bound (expression)
-		double r = 0; // reward value (actual value)
-		String relOp; // relational operator
-		PathExpression pe; // path expression
-
-		JDDNode stateRewards = null, transRewards = null, filter, sol, tmp;
-		StateProbs rewards = null;
-		StateList states = null;
-		double minRes = 0, maxRes = 0;
-		int i;
-
-		// get info from reward operator
-		rs = expr.getRewardStructIndex();
-		relOp = expr.getRelOp();
-		rb = expr.getReward();
-		if (rb != null) {
-			r = rb.evaluateDouble(constantValues, null);
-			if (r < 0)
-				throw new PrismException("Invalid reward bound " + r + " in R[] formula");
-		}
-
-		// get reward info
-		if (model.getNumRewardStructs() == 0)
-			throw new PrismException("Model has no rewards specified");
-		if (rs == null) {
-			stateRewards = model.getStateRewards(0);
-			transRewards = model.getTransRewards(0);
-		} else if (rs instanceof Expression) {
-			i = ((Expression) rs).evaluateInt(constantValues, null);
-			rs = new Integer(i); // for better error reporting below
-			stateRewards = model.getStateRewards(i - 1);
-			transRewards = model.getTransRewards(i - 1);
-		} else if (rs instanceof String) {
-			stateRewards = model.getStateRewards((String) rs);
-			transRewards = model.getTransRewards((String) rs);
-		}
-		if (stateRewards == null || transRewards == null)
-			throw new PrismException("Invalid reward structure index \"" + rs + "\"");
-
-		// check for trivial (i.e. stupid) cases
-		if (rb != null) {
-			if (r == 0 && relOp.equals(">=")) {
-				mainLog.print("\nWarning: checking for reward " + relOp + " " + r
-						+ " - formula trivially satisfies all states\n");
-				JDD.Ref(reach);
-				return reach;
-			} else if (r == 0 && relOp.equals("<")) {
-				mainLog.print("\nWarning: checking for reward " + relOp + " " + r
-						+ " - formula trivially satisfies no states\n");
-				return JDD.Constant(0);
+			if (uTime < lTime) {
+				throw new PrismException("Upper bound must exceed lower bound in time-bounded until formula");
 			}
+		} else {
+			uTime = -1;
 		}
 
-		// print a warning if Rmin/Rmax used
-		if (relOp.equals("min=") || relOp.equals("max=")) {
-			mainLog.print("\nWarning: \"Rmin=?\" and \"Rmax=?\" operators are identical to \"R=?\" for DTMCs\n");
-		}
-
-		// translate filter (if present)
-		filter = null;
-		if (expr.getFilter() != null) {
-			filter = checkExpression(expr.getFilter().getExpression());
-		}
-
-		// check if filter satisfies no states
-		if (filter != null)
-			if (filter.equals(JDD.ZERO)) {
-				// for R=? properties, this is an error
-				if (rb == null) {
-					throw new PrismException("Filter {" + expr.getFilter().getExpression() + "} in R=?[] property satisfies no states");
-				}
-				// otherwise just print a warning
-				else {
-					mainLog.println("\nWarning: Filter {" + expr.getFilter().getExpression()
-							+ "} satisfies no states and is being ignored");
-					JDD.Deref(filter);
-					filter = null;
-				}
-			}
-
-		// compute rewards
-		pe = expr.getPathExpression();
+		// model check operands first
+		b1 = checkExpressionDD(expr1);
 		try {
-			if (pe instanceof PathExpressionTemporal) {
-				switch (((PathExpressionTemporal) pe).getOperator()) {
-				case PathExpressionTemporal.R_C:
-					rewards = checkRewardCumul((PathExpressionTemporal) pe, stateRewards, transRewards);
-					break;
-				case PathExpressionTemporal.R_I:
-					rewards = checkRewardInst((PathExpressionTemporal) pe, stateRewards, transRewards);
-					break;
-				case PathExpressionTemporal.R_F:
-					rewards = checkRewardReach((PathExpressionTemporal) pe, stateRewards, transRewards);
-					break;
-				case PathExpressionTemporal.R_S:
-					rewards = checkRewardSS((PathExpressionTemporal) pe, stateRewards, transRewards, filter, expr.getFilter().getExpression());
-					break;
-				}
-			}
-			if (rewards == null)
-				throw new PrismException("Unrecognised path operator in R operator");
+			b2 = checkExpressionDD(expr2);
 		} catch (PrismException e) {
-			if (filter != null)
-				JDD.Deref(filter);
+			JDD.Deref(b1);
 			throw e;
 		}
 
-		// round off rewards
-		// rewards.roundOff(getPlacesToRoundBy());
+		// print out some info about num states
+		// mainLog.print("\nb1 = " + JDD.GetNumMintermsString(b1,
+		// allDDRowVars.n()));
+		// mainLog.print(" states, b2 = " + JDD.GetNumMintermsString(b2,
+		// allDDRowVars.n()) + " states\n");
 
-		// print out rewards
-		if (verbose) {
-			mainLog.print("\nRewards (non-zero only) for all states:\n");
-			rewards.print(mainLog);
-		}
+		// compute probabilities
 
-		// if a filter was provided, there is some additional output...
-		if (filter != null) {
+		// a trivial case: "U<=0"
+		if (lTime == 0 && uTime == 0) {
+			// prob is 1 in b2 states, 0 otherwise
+			JDD.Ref(b2);
+			probs = new StateProbsMTBDD(b2, model);
+		} else {
 
-			// for non-"R=?"-type properties, print out some rewards
-			if (rb != null) {
-				if (!expr.noFilterRequests()) {
-					mainLog
-							.println("\nWarning: \"{min}\", \"{max}\" only apply to \"R=?\" properties; they are ignored here.");
-				}
-				mainLog.print("\nRewards (non-zero only) for states satisfying " + expr.getFilter().getExpression() + ":\n");
-				rewards.printFiltered(mainLog, filter);
-			}
+			// break down into different cases to compute probabilities
 
-			// for "R=?"-type properties...
-			if (rb == null) {
-				// compute/print min info
-				if (expr.filterMinRequested()) {
-					minRes = rewards.minOverBDD(filter);
-					mainLog.print("\nMinimum reward for states satisfying " + expr.getFilter().getExpression() + ": " + minRes + "\n");
-					tmp = rewards.getBDDFromInterval(minRes - termCritParam, minRes + termCritParam);
-					JDD.Ref(filter);
-					tmp = JDD.And(tmp, filter);
-					states = new StateListMTBDD(tmp, model);
-					mainLog.print("There are " + states.size() + " states with this minimum reward (+/- "
-							+ termCritParam + ")");
-					if (!verbose && states.size() > 10) {
-						mainLog
-								.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
-						states.print(mainLog, 10);
-					} else {
-						mainLog.print(":\n");
-						states.print(mainLog);
+			// >= lTime
+			if (uTime == -1) {
+				// check for special case of lTime == 0, this is actually an
+				// unbounded until
+				if (lTime == 0) {
+					// compute probs
+					try {
+						probs = computeUntilProbs(trans, trans01, b1, b2);
+					} catch (PrismException e) {
+						JDD.Deref(b1);
+						JDD.Deref(b2);
+						throw e;
 					}
-					JDD.Deref(tmp);
-				}
-				// compute/print min info
-				if (expr.filterMaxRequested()) {
-					maxRes = rewards.maxOverBDD(filter);
-					mainLog.print("\nMaximum reward for states satisfying " + expr.getFilter().getExpression() + ": " + maxRes + "\n");
-					tmp = rewards.getBDDFromInterval(maxRes - termCritParam, maxRes + termCritParam);
-					JDD.Ref(filter);
-					tmp = JDD.And(tmp, filter);
-					states = new StateListMTBDD(tmp, model);
-					mainLog.print("There are " + states.size() + " states with this maximum reward (+/- "
-							+ termCritParam + ")");
-					if (!verbose && states.size() > 10) {
-						mainLog
-								.print(".\nThe first 10 states are displayed below. To view them all, use verbose mode.\n");
-						states.print(mainLog, 10);
-					} else {
-						mainLog.print(":\n");
-						states.print(mainLog);
+				} else {
+					// compute unbounded until probs
+					try {
+						tmpProbs = computeUntilProbs(trans, trans01, b1, b2);
+					} catch (PrismException e) {
+						JDD.Deref(b1);
+						JDD.Deref(b2);
+						throw e;
 					}
-					JDD.Deref(tmp);
+					// compute bounded until probs
+					try {
+						probs = computeBoundedUntilProbs(trans, trans01, b1, b1, lTime, tmpProbs);
+					} catch (PrismException e) {
+						tmpProbs.clear();
+						JDD.Deref(b1);
+						JDD.Deref(b2);
+						throw e;
+					}
+					tmpProbs.clear();
 				}
 			}
-		}
-
-		// for R=? properties...
-		// if there are multiple initial states or if there is a filter
-		// satisfying
-		// multiple states with no {min}/{max}/etc., print a warning...
-		if (rb == null) {
-			if (filter == null) {
-				if (model.getNumStartStates() > 1) {
-					mainLog
-							.print("\nWarning: There are multiple initial states; the result of model checking is for the first one: ");
-					model.getStartStates().print(mainLog, 1);
-				}
-			} else if (expr.noFilterRequests()) {
-				StateListMTBDD filterStates = new StateListMTBDD(filter, model);
-				if (filterStates.size() > 1) {
-					mainLog.print("\nWarning: The filter {" + expr.getFilter().getExpression() + "} is satisfied by "
-							+ filterStates.size() + " states.\n");
-					mainLog.print("The result of model checking is for the first of these: ");
-					filterStates.print(mainLog, 1);
+			// <= uTime
+			else if (lTime == 0) {
+				// nb: uTime != 0 since would be caught above (trivial case)
+				JDD.Ref(b1);
+				JDD.Ref(b2);
+				tmp = JDD.And(b1, JDD.Not(b2));
+				try {
+					probs = computeBoundedUntilProbs(trans, trans01, b2, tmp, uTime, null);
+				} catch (PrismException e) {
+					JDD.Deref(tmp);
+					JDD.Deref(b1);
+					JDD.Deref(b2);
+					throw e;
 				}
+				JDD.Deref(tmp);
 			}
-		}
-
-		// compute result of property
-		// if there's a bound, get set of satisfying states
-		if (rb != null) {
-			sol = rewards.getBDDFromInterval(relOp, r);
-			// remove unreachable states from solution
-			JDD.Ref(reach);
-			sol = JDD.And(sol, reach);
-		}
-		// if there's no bound, result will be a reward
-		else {
-			// just store empty set for sol
-			sol = JDD.Constant(0);
-			// use filter if present
-			if (filter != null) {
-				if (expr.filterMinRequested())
-					numericalRes = minRes;
-				else if (expr.filterMaxRequested())
-					numericalRes = maxRes;
-				else
-					numericalRes = rewards.firstFromBDD(filter);
-			}
-			// otherwise use initial state
-			else {
-				numericalRes = rewards.firstFromBDD(start);
-			}
-		}
-
-		// free vector
-		rewards.clear();
-
-		// derefs
-		if (filter != null)
-			JDD.Deref(filter);
-
-		return sol;
-	}
-
-	// S operator
-
-	private JDDNode checkExpressionSteadyState(ExpressionSS expr) throws PrismException
-	{
-		// formula stuff
-		Expression pb; // probability bound (expression)
-		double p = 0; // probability value (actual value)
-		String relOp; // relational operator
-		// bscc stuff
-		Vector vectBSCCs;
-		JDDNode notInBSCCs;
-		// mtbdd stuff
-		JDDNode b, filter, bscc, diags, emb, sol, tmp;
-		// other stuff
-		StateProbs probs = null, totalProbs = null;
-		int i, n;
-		double d, probBSCCs[];
-		boolean sat;
-
-		// get info from steady-state operator
-		relOp = expr.getRelOp();
-		pb = expr.getProb();
-		if (pb != null) {
-			p = pb.evaluateDouble(constantValues, null);
-			if (p < 0 || p > 1)
-				throw new PrismException("Invalid probability bound " + p + " in until formula");
-		}
-
-		// check for trivial (i.e. stupid) cases
-		if (pb != null) {
-			if ((p == 0 && relOp.equals(">=")) || (p == 1 && relOp.equals("<="))) {
-				mainLog.print("\nWarning: checking for probability " + relOp + " " + p
-						+ " - formula trivially satisfies all states\n");
-				JDD.Ref(reach);
-				return reach;
-			} else if ((p == 0 && relOp.equals("<")) || (p == 1 && relOp.equals(">"))) {
-				mainLog.print("\nWarning: checking for probability " + relOp + " " + p
-						+ " - formula trivially satisfies no states\n");
-				return JDD.Constant(0);
-			}
-		}
-
-		// translate b into a bdd
-		// (i.e. model check argument first)
-		b = checkExpression(expr.getExpression());
-
-		// translate filter (if present)
-		filter = null;
-		if (expr.getFilter() != null) {
-			filter = checkExpression(expr.getFilter().getExpression());
-			// if filter is empty, warn and ignore
-			if (filter.equals(JDD.ZERO)) {
-				mainLog.println("\nWarning: Filter \"" + expr.getFilter().getExpression()
-						+ "\" satisfies no states and is being ignored");
-				JDD.Deref(filter);
-				filter = null;
-			}
-		}
-
-		// print out some basic info
-		mainLog.print("\nCSL Steady State:\n");
-		mainLog.print("\nb = " + JDD.GetNumMintermsString(b, allDDRowVars.n()) + " states\n");
-
-		// compute bottom strongly connected components (bsccs)
-		if (bsccComp) {
-			mainLog.print("\nComputing (B)SCCs...");
-			sccComputer.computeBSCCs();
-			vectBSCCs = sccComputer.getVectBSCCs();
-			notInBSCCs = sccComputer.getNotInBSCCs();
-			n = vectBSCCs.size();
-		}
-		// unless we've been told to skip it
-		else {
-			mainLog.println("\nSkipping BSCC computation...");
-			vectBSCCs = new Vector();
-			JDD.Ref(reach);
-			vectBSCCs.add(reach);
-			notInBSCCs = JDD.Constant(0);
-			n = 1;
-		}
-
-		// compute steady state for each bscc...
-		probBSCCs = new double[n];
-		for (i = 0; i < n; i++) {
-
-			mainLog.println("\nComputing steady state probabilities for BSCC " + (i + 1));
-
-			// get bscc
-			bscc = (JDDNode) vectBSCCs.elementAt(i);
-
-			// compute steady state probabilities
-			try {
-				probs = computeSteadyStateProbs(trans, bscc);
-			} catch (PrismException e) {
-				JDD.Deref(b);
-				if (filter != null)
-					JDD.Deref(filter);
-				for (i = 0; i < n; i++) {
-					JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
-				}
-				JDD.Deref(notInBSCCs);
-				throw e;
-			}
-
-			// round off probabilities
-			// probs.roundOff(getPlacesToRoundBy());
-
-			// print out probabilities
-			if (verbose) {
-				mainLog.print("\nBSCC " + (i + 1) + " Steady-State Probabilities: \n");
-				probs.print(mainLog);
-			}
-
-			// print out filtered probabilities
-			if (filter != null) {
-				mainLog.print("\nBSCC " + (i + 1) + " Steady-State Probabilities (" + expr.getFilter().getExpression() + "): \n");
-				probs.printFiltered(mainLog, filter);
-			}
-
-			// sum probabilities over bdd b
-			d = probs.sumOverBDD(b);
-			probBSCCs[i] = d;
-			mainLog.print("\nBSCC " + (i + 1) + " Probability: " + d + "\n");
-
-			// free vector
-			probs.clear();
-		}
-
-		// if every state is in a bscc, it's much easier...
-		if (notInBSCCs.equals(JDD.ZERO)) {
-
-			mainLog.println("\nAll states are in a BSCC (so no reachability probabilities computed)");
-
-			// there's more efficient ways to do this if we just create the
-			// solution bdd directly
-			// but we actually build the prob vector so it can be printed out if
-			// necessary
-			tmp = JDD.Constant(0);
-			for (i = 0; i < n; i++) {
-				bscc = (JDDNode) vectBSCCs.elementAt(i);
-				JDD.Ref(bscc);
-				tmp = JDD.Apply(JDD.PLUS, tmp, JDD.Apply(JDD.TIMES, JDD.Constant(probBSCCs[i]), bscc));
-			}
-			probs = new StateProbsMTBDD(tmp, model);
-
-			// print out probabilities
-			if (verbose) {
-				mainLog.print("\nProbabilities for CSL steady-state formula: \n");
-				probs.print(mainLog);
-			}
-
-			// print out filtered probabilities
-			if (filter != null) {
-				mainLog.print("\nProbabilities for CSL steady-state formula (" + expr.getFilter().getExpression() + "): \n");
-				probs.printFiltered(mainLog, filter);
-			}
-
-			// compute bdd of satisfying states
-			if (pb != null) {
-				sol = probs.getBDDFromInterval(relOp, p);
-				// remove unreachable states from solution
-				JDD.Ref(reach);
-				sol = JDD.And(sol, reach);
-			}
-			// if there's no bound, result will be a probability (just store
-			// empty set for sol)
-			else {
-				sol = JDD.Constant(0);
-				// use filter if present
-				if (filter != null) {
-					numericalRes = probs.firstFromBDD(filter);
-				}
-				// otherwise use initial state
-				else {
-					numericalRes = probs.firstFromBDD(start);
-				}
-			}
-
-			// free vector
-			probs.clear();
-		}
-		// otherwise we have to do more work...
-		else {
-
-			// initialise total probabilities vector
-			switch (engine) {
-			case Prism.MTBDD:
-				totalProbs = new StateProbsMTBDD(JDD.Constant(0), model);
-				break;
-			case Prism.SPARSE:
-				totalProbs = new StateProbsDV(new DoubleVector((int) numStates), model);
-				break;
-			case Prism.HYBRID:
-				totalProbs = new StateProbsDV(new DoubleVector((int) numStates), model);
-				break;
-			}
-
-			// compute diagonals
-			JDD.Ref(trans);
-			diags = JDD.SumAbstract(trans, allDDColVars);
-
-			// compute embedded markov chain
-			JDD.Ref(trans);
-			JDD.Ref(diags);
-			emb = JDD.Apply(JDD.DIVIDE, trans, diags);
-
-			mainLog.println("\nDiagonals vector: " + JDD.GetInfoString(diags, allDDRowVars.n()));
-			mainLog.println("Embedded Markov chain: " + JDD.GetInfoString(emb, allDDRowVars.n() * 2));
-
-			// compute probabilities of reaching each bscc...
-			for (i = 0; i < n; i++) {
-
-				// skip bsccs with zero probability
-				if (probBSCCs[i] == 0.0)
-					continue;
-
-				mainLog.println("\nComputing probabilities of reaching BSCC " + (i + 1));
-
-				// get bscc
-				bscc = (JDDNode) vectBSCCs.elementAt(i);
-
-				// compute probabilities
-				try {
-					probs = computeUntilProbs(emb, trans01, notInBSCCs, bscc);
-				} catch (PrismException e) {
-					JDD.Deref(diags);
-					JDD.Deref(emb);
-					JDD.Deref(b);
-					if (filter != null)
-						JDD.Deref(filter);
-					for (i = 0; i < n; i++) {
-						JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
-					}
-					JDD.Deref(notInBSCCs);
-					totalProbs.clear();
-					throw e;
-				}
-
-				// round off probabilities
-				// probs.roundOff(getPlacesToRoundBy());
-
-				// print out probabilities
-				if (verbose) {
-					mainLog.print("\nBSCC " + (i + 1) + " Reachability Probabilities: \n");
-					probs.print(mainLog);
-				}
-
-				// print out filtered probabilities
-				if (filter != null) {
-					mainLog.print("\nBSCC " + (i + 1) + " Reachability Probabilities (" + expr.getFilter().getExpression() + "): \n");
-					probs.printFiltered(mainLog, filter);
-				}
-
-				// times by bscc prob, add to total
-				probs.timesConstant(probBSCCs[i]);
-				totalProbs.add(probs);
-
-				// free vector
-				probs.clear();
-			}
-
-			// print out probabilities
-			if (verbose) {
-				mainLog.print("\nProbabilities for CSL steady-state formula: \n");
-				totalProbs.print(mainLog);
-			}
-
-			// print out filtered probabilities
-			if (filter != null) {
-				mainLog.print("\nProbabilities for CSL steady-state formula (" + expr.getFilter().getExpression() + "): \n");
-				totalProbs.printFiltered(mainLog, filter);
-			}
-
-			// compute bdd of satisfying states
-			if (pb != null) {
-				sol = totalProbs.getBDDFromInterval(relOp, p);
-				// remove unreachable states from solution
-				JDD.Ref(reach);
-				sol = JDD.And(sol, reach);
-			}
-			// if there's no bound, result will be a probability (just store
-			// empty set for sol)
-			else {
-				sol = JDD.Constant(0);
-				// use filter if present
-				if (filter != null) {
-					numericalRes = totalProbs.firstFromBDD(filter);
-				}
-				// otherwise use initial state
-				else {
-					numericalRes = totalProbs.firstFromBDD(start);
-				}
-			}
-
-			// free vector
-			totalProbs.clear();
-
-			// derefs
-			JDD.Deref(diags);
-			JDD.Deref(emb);
-		}
-
-		// derefs
-		JDD.Deref(b);
-		if (filter != null)
-			JDD.Deref(filter);
-		for (i = 0; i < n; i++) {
-			JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
-		}
-		JDD.Deref(notInBSCCs);
-
-		return sol;
-	}
-
-	// next
-
-	private StateProbs checkProbNext(PathExpressionTemporal pe) throws PrismException
-	{
-		Expression expr;
-		JDDNode b, diags, emb;
-		StateProbs probs = null;
-
-		// get operand
-		if (!(pe.getOperand2() instanceof PathExpressionExpr))
-			throw new PrismException("Invalid path formula");
-		expr = ((PathExpressionExpr) pe.getOperand2()).getExpression();
-
-		// model check operand first
-		b = checkExpression(expr);
-
-		// print out some info about num states
-		// mainLog.print("\nb = " + JDD.GetNumMintermsString(b,
-		// allDDRowVars.n()) + " states\n");
-
-		// compute diagonals
-		JDD.Ref(trans);
-		diags = JDD.SumAbstract(trans, allDDColVars);
-
-		// compute embedded markov chain
-		JDD.Ref(trans);
-		JDD.Ref(diags);
-		emb = JDD.Apply(JDD.DIVIDE, trans, diags);
-
-		mainLog.println("\nDiagonals vector: " + JDD.GetInfoString(diags, allDDRowVars.n()));
-		mainLog.println("Embedded Markov chain: " + JDD.GetInfoString(emb, allDDRowVars.n() * 2));
-
-		// compute probabilities
-		probs = computeNextProbs(emb, b);
-
-		// derefs
-		JDD.Deref(b);
-		JDD.Deref(diags);
-		JDD.Deref(emb);
-
-		return probs;
-	}
-
-	// bounded until
-
-	private StateProbs checkProbBoundedUntil(PathExpressionTemporal pe) throws PrismException
-	{
-		Expression expr1, expr2;
-		double lTime, uTime; // time bounds
-		Expression expr;
-		JDDNode b1, b2, tmp, diags, emb;
-		StateProbs tmpProbs = null, probs = null;
-
-		// get operands
-		if (!(pe.getOperand1() instanceof PathExpressionExpr) || !(pe.getOperand2() instanceof PathExpressionExpr))
-			throw new PrismException("Invalid path formula");
-		expr1 = ((PathExpressionExpr) pe.getOperand1()).getExpression();
-		expr2 = ((PathExpressionExpr) pe.getOperand2()).getExpression();
-
-		// get info from bounded until
-
-		// lower bound is 0 if not specified
-		// (i.e. if until is of form U<=t)
-		expr = pe.getLowerBound();
-		if (expr != null) {
-			lTime = expr.evaluateDouble(constantValues, null);
-			if (lTime < 0) {
-				throw new PrismException("Invalid lower bound " + lTime + " in time-bounded until formula");
-			}
-		} else {
-			lTime = 0;
-		}
-		// upper bound is -1 if not specified
-		// (i.e. if until is of form U>=t)
-		expr = pe.getUpperBound();
-		if (expr != null) {
-			uTime = expr.evaluateDouble(constantValues, null);
-			if (uTime < 0) {
-				throw new PrismException("Invalid upper bound " + uTime + " in time-bounded until formula");
-			}
-			if (uTime < lTime) {
-				throw new PrismException("Upper bound must exceed lower bound in time-bounded until formula");
-			}
-		} else {
-			uTime = -1;
-		}
-
-		// model check operands first
-		b1 = checkExpression(expr1);
-		try {
-			b2 = checkExpression(expr2);
-		} catch (PrismException e) {
-			JDD.Deref(b1);
-			throw e;
-		}
-
-		// print out some info about num states
-		// mainLog.print("\nb1 = " + JDD.GetNumMintermsString(b1,
-		// allDDRowVars.n()));
-		// mainLog.print(" states, b2 = " + JDD.GetNumMintermsString(b2,
-		// allDDRowVars.n()) + " states\n");
-
-		// compute probabilities
-
-		// a trivial case: "U<=0"
-		if (lTime == 0 && uTime == 0) {
-			// prob is 1 in b2 states, 0 otherwise
-			JDD.Ref(b2);
-			probs = new StateProbsMTBDD(b2, model);
-		} else {
-
-			// break down into different cases to compute probabilities
-
-			// >= lTime
-			if (uTime == -1) {
-				// check for special case of lTime == 0, this is actually an
-				// unbounded until
-				if (lTime == 0) {
-					// compute diagonals
-					JDD.Ref(trans);
-					diags = JDD.SumAbstract(trans, allDDColVars);
-					// compute embedded markov chain
-					JDD.Ref(trans);
-					JDD.Ref(diags);
-					emb = JDD.Apply(JDD.DIVIDE, trans, diags);
-					mainLog.println("\nDiagonals vector: " + JDD.GetInfoString(diags, allDDRowVars.n()));
-					mainLog.println("Embedded Markov chain: " + JDD.GetInfoString(emb, allDDRowVars.n() * 2));
-					// compute probs
-					try {
-						probs = computeUntilProbs(emb, trans01, b1, b2);
-					} catch (PrismException e) {
-						JDD.Deref(diags);
-						JDD.Deref(emb);
-						JDD.Deref(b1);
-						JDD.Deref(b2);
-						throw e;
-					}
-					JDD.Deref(diags);
-					JDD.Deref(emb);
-				} else {
-					// compute diagonals
-					JDD.Ref(trans);
-					diags = JDD.SumAbstract(trans, allDDColVars);
-					// compute embedded markov chain
-					JDD.Ref(trans);
-					JDD.Ref(diags);
-					emb = JDD.Apply(JDD.DIVIDE, trans, diags);
-					mainLog.println("\nDiagonals vector: " + JDD.GetInfoString(diags, allDDRowVars.n()));
-					mainLog.println("Embedded Markov chain: " + JDD.GetInfoString(emb, allDDRowVars.n() * 2));
-					// compute unbounded until probs
-					try {
-						tmpProbs = computeUntilProbs(emb, trans01, b1, b2);
-					} catch (PrismException e) {
-						JDD.Deref(diags);
-						JDD.Deref(emb);
-						JDD.Deref(b1);
-						JDD.Deref(b2);
-						throw e;
-					}
-					JDD.Deref(diags);
-					JDD.Deref(emb);
-					// compute bounded until probs
-					try {
-						probs = computeBoundedUntilProbs(trans, trans01, b1, b1, lTime, tmpProbs);
-					} catch (PrismException e) {
-						tmpProbs.clear();
-						JDD.Deref(b1);
-						JDD.Deref(b2);
-						throw e;
-					}
-					tmpProbs.clear();
-				}
-			}
-			// <= uTime
-			else if (lTime == 0) {
-				// nb: uTime != 0 since would be caught above (trivial case)
-				JDD.Ref(b1);
-				JDD.Ref(b2);
-				tmp = JDD.And(b1, JDD.Not(b2));
-				try {
-					probs = computeBoundedUntilProbs(trans, trans01, b2, tmp, uTime, null);
-				} catch (PrismException e) {
-					JDD.Deref(tmp);
-					JDD.Deref(b1);
-					JDD.Deref(b2);
-					throw e;
-				}
-				JDD.Deref(tmp);
-			}
-			// [lTime,uTime] (including where lTime == uTime)
+			// [lTime,uTime] (including where lTime == uTime)
 			else {
 				JDD.Ref(b1);
 				JDD.Ref(b2);
@@ -1354,697 +206,92 @@ public class StochModelChecker implements ModelChecker
 			}
 		}
 
-		// derefs
-		JDD.Deref(b1);
-		JDD.Deref(b2);
-
-		return probs;
-	}
-
-	// until (unbounded)
-
-	private StateProbs checkProbUntil(PathExpressionTemporal pe, Expression pb, double p) throws PrismException
-	{
-		Expression expr1, expr2;
-		JDDNode b1, b2, diags, emb;
-		StateProbs probs = null;
-
-		// get operands
-		if (!(pe.getOperand1() instanceof PathExpressionExpr) || !(pe.getOperand2() instanceof PathExpressionExpr))
-			throw new PrismException("Invalid path formula");
-		expr1 = ((PathExpressionExpr) pe.getOperand1()).getExpression();
-		expr2 = ((PathExpressionExpr) pe.getOperand2()).getExpression();
-
-		// model check operands first
-		b1 = checkExpression(expr1);
-		try {
-			b2 = checkExpression(expr2);
-		} catch (PrismException e) {
-			JDD.Deref(b1);
-			throw e;
-		}
-
-		// print out some info about num states
-		mainLog.print("\nb1 = " + JDD.GetNumMintermsString(b1, allDDRowVars.n()));
-		mainLog.print(" states, b2 = " + JDD.GetNumMintermsString(b2, allDDRowVars.n()) + " states\n");
-
-		// if p is 0 or 1 and precomputation algorithms are enabled, compute
-		// probabilities qualitatively
-		if (pb != null && ((p == 0) || (p == 1)) && precomp) {
-			mainLog.print("\nWarning: probability bound in formula is " + p
-					+ " so exact probabilities may not be computed\n");
-			probs = computeUntilProbsQual(trans01, b1, b2, p);
-		}
-		// otherwise actually compute probabilities
-		else {
-			// compute diagonals
-			JDD.Ref(trans);
-			diags = JDD.SumAbstract(trans, allDDColVars);
-
-			// compute embedded markov chain
-			JDD.Ref(trans);
-			JDD.Ref(diags);
-			emb = JDD.Apply(JDD.DIVIDE, trans, diags);
-
-			mainLog.println("\nDiagonals vector: " + JDD.GetInfoString(diags, allDDRowVars.n()));
-			mainLog.println("Embedded Markov chain: " + JDD.GetInfoString(emb, allDDRowVars.n() * 2));
-
-			// compute probabilities
-			try {
-				probs = computeUntilProbs(emb, trans01, b1, b2);
-			} catch (PrismException e) {
-				JDD.Deref(b1);
-				JDD.Deref(b2);
-				JDD.Deref(diags);
-				JDD.Deref(emb);
-				throw e;
-			}
-
-			// derefs
-			JDD.Deref(diags);
-			JDD.Deref(emb);
-		}
-
-		// derefs
-		JDD.Deref(b1);
-		JDD.Deref(b2);
-
-		return probs;
-	}
-
-	// bounded future (eventually)
-	// F[t1,t2] phi == true U[t1,t2] phi
-
-	private StateProbs checkProbBoundedFuture(PathExpressionTemporal pe) throws PrismException
-	{
-		PathExpressionTemporal pe2;
-		pe2 = new PathExpressionTemporal(PathExpressionTemporal.P_U, new PathExpressionExpr(Expression.True()), pe
-				.getOperand2(), pe.getLowerBound(), pe.getUpperBound());
-		return checkProbBoundedUntil(pe2);
-	}
-
-	// future (eventually)
-	// F phi == true U phi
-
-	private StateProbs checkProbFuture(PathExpressionTemporal pe, Expression pb, double p) throws PrismException
-	{
-		PathExpressionTemporal pe2;
-		pe2 = new PathExpressionTemporal(PathExpressionTemporal.P_U, new PathExpressionExpr(Expression.True()), pe
-				.getOperand2());
-		return checkProbUntil(pe2, pb, p);
-	}
-
-	// bounded global (always)
-	// F<=k phi == true U<=k phi
-	// P(G<=k phi) == 1-P(true U<=k !phi)
-
-	private StateProbs checkProbBoundedGlobal(PathExpressionTemporal pe) throws PrismException
-	{
-		PathExpressionTemporal pe2;
-		StateProbs probs;
-		pe2 = new PathExpressionTemporal(PathExpressionTemporal.P_U, new PathExpressionExpr(Expression.True()),
-				new PathExpressionExpr(Expression.Not(((PathExpressionExpr) pe.getOperand2()).getExpression())), pe
-						.getLowerBound(), pe.getUpperBound());
-		probs = checkProbBoundedUntil(pe2);
-		probs.subtractFromOne();
-		return probs;
-	}
-
-	// global (always)
-	// G phi == !(true U !phi)
-	// P(G phi) == 1-P(true U !phi)
-
-	private StateProbs checkProbGlobal(PathExpressionTemporal pe, Expression pb, double p) throws PrismException
-	{
-		PathExpressionTemporal pe2;
-		StateProbs probs;
-		pe2 = new PathExpressionTemporal(PathExpressionTemporal.P_U, new PathExpressionExpr(Expression.True()),
-				new PathExpressionExpr(Expression.Not(((PathExpressionExpr) pe.getOperand2()).getExpression())));
-		probs = checkProbUntil(pe2, pb, p);
-		probs.subtractFromOne();
-		return probs;
-	}
-
-	// cumulative reward
-
-	private StateProbs checkRewardCumul(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards)
-			throws PrismException
-	{
-		double time; // time
-		StateProbs rewards = null;
-
-		// get info from inst reward
-		time = pe.getUpperBound().evaluateDouble(constantValues, null);
-		if (time < 0) {
-			throw new PrismException("Invalid time bound " + time + " in cumulative reward formula");
-		}
-
-		// compute rewards
-
-		// a trivial case: "<=0"
-		if (time == 0) {
-			rewards = new StateProbsMTBDD(JDD.Constant(0), model);
-		} else {
-			// compute rewards
-			try {
-				rewards = computeCumulRewards(trans, trans01, stateRewards, transRewards, time);
-			} catch (PrismException e) {
-				throw e;
-			}
-		}
-
-		return rewards;
-	}
-
-	// inst reward
-
-	private StateProbs checkRewardInst(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards)
-			throws PrismException
-	{
-		double time; // time
-		StateProbs sr = null, rewards = null;
-
-		// get info from inst reward
-		time = pe.getUpperBound().evaluateDouble(constantValues, null);
-		if (time < 0) {
-			throw new PrismException("Invalid bound " + time + " in instantaneous reward property");
-		}
-
-		// compute rewards
-
-		// a trivial case: "=0"
-		if (time == 0) {
-			JDD.Ref(stateRewards);
-			rewards = new StateProbsMTBDD(stateRewards, model);
-		} else {
-			// convert state rewards vector to appropriate type (depending on
-			// engine)
-			switch (engine) {
-			case Prism.MTBDD:
-				JDD.Ref(stateRewards);
-				sr = new StateProbsMTBDD(stateRewards, model);
-				break;
-			case Prism.SPARSE:
-				sr = new StateProbsDV(stateRewards, model);
-				break;
-			case Prism.HYBRID:
-				sr = new StateProbsDV(stateRewards, model);
-				break;
-			}
-			// and for the computation, we can reuse the computation for
-			// time-bounded until formulae
-			// which is nice
-			try {
-				rewards = computeBoundedUntilProbs(trans, trans01, reach, reach, time, sr);
-			} catch (PrismException e) {
-				sr.clear();
-				throw e;
-			}
-			sr.clear();
-		}
-
-		return rewards;
-	}
-
-	// reach reward
-
-	private StateProbs checkRewardReach(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards)
-			throws PrismException
-	{
-		Expression expr;
-		JDDNode b, diags, emb, newStateRewards;
-		StateProbs rewards = null;
-
-		// get operand
-		if (!(pe.getOperand2() instanceof PathExpressionExpr))
-			throw new PrismException("Invalid path formula");
-		expr = ((PathExpressionExpr) pe.getOperand2()).getExpression();
-
-		// model check operand first
-		b = checkExpression(expr);
-
-		// print out some info about num states
-		// mainLog.print("\nb = " + JDD.GetNumMintermsString(b,
-		// allDDRowVars.n()) + " states\n");
-
-		// compute diagonals
-		JDD.Ref(trans);
-		diags = JDD.SumAbstract(trans, allDDColVars);
-
-		// compute embedded markov chain
-		JDD.Ref(trans);
-		JDD.Ref(diags);
-		emb = JDD.Apply(JDD.DIVIDE, trans, diags);
-
-		mainLog.println("\nDiagonals vector: " + JDD.GetInfoString(diags, allDDRowVars.n()));
-		mainLog.println("Embedded Markov chain: " + JDD.GetInfoString(emb, allDDRowVars.n() * 2));
-
-		// convert rewards
-		JDD.Ref(stateRewards);
-		JDD.Ref(diags);
-		newStateRewards = JDD.Apply(JDD.DIVIDE, stateRewards, diags);
-
-		// compute rewards
-		try {
-			rewards = computeReachRewards(emb, trans01, newStateRewards, transRewards, b);
-		} catch (PrismException e) {
-			JDD.Deref(b);
-			JDD.Deref(emb);
-			JDD.Deref(diags);
-			JDD.Deref(newStateRewards);
-			throw e;
-		}
-
-		// derefs
-		JDD.Deref(b);
-		JDD.Deref(emb);
-		JDD.Deref(diags);
-		JDD.Deref(newStateRewards);
-
-		return rewards;
-	}
-
-	// steady state reward
-
-	private StateProbs checkRewardSS(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards, JDDNode filter,
-			Expression filterName) throws PrismException
-	{
-		// bscc stuff
-		Vector vectBSCCs;
-		JDDNode notInBSCCs;
-		// mtbdd stuff
-		JDDNode newStateRewards, bscc, diags, emb, tmp;
-		// other stuff
-		StateProbs probs = null, rewards = null;
-		int i, n;
-		double d, rewBSCCs[];
-
-		// compute rewards corresponding to each state
-		JDD.Ref(trans);
-		JDD.Ref(transRewards);
-		newStateRewards = JDD.SumAbstract(JDD.Apply(JDD.TIMES, trans, transRewards), allDDColVars);
-		JDD.Ref(stateRewards);
-		newStateRewards = JDD.Apply(JDD.PLUS, newStateRewards, stateRewards);
-
-		// compute bottom strongly connected components (bsccs)
-		if (bsccComp) {
-			mainLog.print("\nComputing (B)SCCs...");
-			sccComputer.computeBSCCs();
-			vectBSCCs = sccComputer.getVectBSCCs();
-			notInBSCCs = sccComputer.getNotInBSCCs();
-			n = vectBSCCs.size();
-		}
-		// unless we've been told to skip it
-		else {
-			mainLog.println("\nSkipping BSCC computation...");
-			vectBSCCs = new Vector();
-			JDD.Ref(reach);
-			vectBSCCs.add(reach);
-			notInBSCCs = JDD.Constant(0);
-			n = 1;
-		}
-
-		// compute steady state for each bscc...
-		rewBSCCs = new double[n];
-		for (i = 0; i < n; i++) {
-
-			mainLog.println("\nComputing steady state probabilities for BSCC " + (i + 1));
-
-			// get bscc
-			bscc = (JDDNode) vectBSCCs.elementAt(i);
-
-			// compute steady state probabilities
-			try {
-				probs = computeSteadyStateProbs(trans, bscc);
-			} catch (PrismException e) {
-				JDD.Deref(newStateRewards);
-				for (i = 0; i < n; i++) {
-					JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
-				}
-				JDD.Deref(notInBSCCs);
-				throw e;
-			}
-
-			// round off probabilities
-			// probs.roundOff(getPlacesToRoundBy());
-
-			// print out probabilities
-			if (verbose) {
-				mainLog.print("\nBSCC " + (i + 1) + " Steady-State Probabilities: \n");
-				probs.print(mainLog);
-			}
-
-			// print out filtered probabilities
-			if (filter != null) {
-				mainLog.print("\nBSCC " + (i + 1) + " Steady-State Probabilities (" + filterName + "): \n");
-				probs.printFiltered(mainLog, filter);
-			}
-
-			// do weighted sum of probabilities and rewards
-			JDD.Ref(bscc);
-			JDD.Ref(newStateRewards);
-			tmp = JDD.Apply(JDD.TIMES, bscc, newStateRewards);
-			d = probs.sumOverMTBDD(tmp);
-			rewBSCCs[i] = d;
-			mainLog.print("\nBSCC " + (i + 1) + " Reward: " + d + "\n");
-			JDD.Deref(tmp);
-
-			// free vector
-			probs.clear();
-		}
-
-		// if every state is in a bscc, it's much easier...
-		if (notInBSCCs.equals(JDD.ZERO)) {
-
-			mainLog.println("\nAll states are in a BSCC (so no reachability probabilities computed)");
-
-			// build the reward vector
-			tmp = JDD.Constant(0);
-			for (i = 0; i < n; i++) {
-				bscc = (JDDNode) vectBSCCs.elementAt(i);
-				JDD.Ref(bscc);
-				tmp = JDD.Apply(JDD.PLUS, tmp, JDD.Apply(JDD.TIMES, JDD.Constant(rewBSCCs[i]), bscc));
-			}
-			rewards = new StateProbsMTBDD(tmp, model);
-		}
-		// otherwise we have to do more work...
-		else {
-
-			// initialise rewards vector
-			switch (engine) {
-			case Prism.MTBDD:
-				rewards = new StateProbsMTBDD(JDD.Constant(0), model);
-				break;
-			case Prism.SPARSE:
-				rewards = new StateProbsDV(new DoubleVector((int) numStates), model);
-				break;
-			case Prism.HYBRID:
-				rewards = new StateProbsDV(new DoubleVector((int) numStates), model);
-				break;
-			}
-
-			// compute diagonals
-			JDD.Ref(trans);
-			diags = JDD.SumAbstract(trans, allDDColVars);
-
-			// compute embedded markov chain
-			JDD.Ref(trans);
-			JDD.Ref(diags);
-			emb = JDD.Apply(JDD.DIVIDE, trans, diags);
-
-			mainLog.println("\nDiagonals vector: " + JDD.GetInfoString(diags, allDDRowVars.n()));
-			mainLog.println("Embedded Markov chain: " + JDD.GetInfoString(emb, allDDRowVars.n() * 2));
-
-			// compute probabilities of reaching each bscc...
-			for (i = 0; i < n; i++) {
-
-				// skip bsccs with zero reward
-				if (rewBSCCs[i] == 0.0)
-					continue;
-
-				mainLog.println("\nComputing probabilities of reaching BSCC " + (i + 1));
-
-				// get bscc
-				bscc = (JDDNode) vectBSCCs.elementAt(i);
-
-				// compute probabilities
-				try {
-					probs = computeUntilProbs(emb, trans01, notInBSCCs, bscc);
-				} catch (PrismException e) {
-					JDD.Deref(diags);
-					JDD.Deref(emb);
-					throw e;
-				}
-
-				// round off probabilities
-				// probs.roundOff(getPlacesToRoundBy());
-
-				// print out probabilities
-				if (verbose) {
-					mainLog.print("\nBSCC " + (i + 1) + " Reachability Probabilities: \n");
-					probs.print(mainLog);
-				}
-
-				// print out filtered probabilities
-				if (filter != null) {
-					mainLog.print("\nBSCC " + (i + 1) + " Reachability Probabilities (" + filterName + "): \n");
-					probs.printFiltered(mainLog, filter);
-				}
-
-				// times by bscc reward, add to total
-				probs.timesConstant(rewBSCCs[i]);
-				rewards.add(probs);
-
-				// free vector
-				probs.clear();
-			}
-
-			// derefs
-			JDD.Deref(diags);
-			JDD.Deref(emb);
-		}
-
-		// derefs
-		JDD.Deref(newStateRewards);
-		if (filter != null)
-			JDD.Deref(filter);
-		for (i = 0; i < n; i++) {
-			JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
-		}
-		JDD.Deref(notInBSCCs);
-
-		return rewards;
-	}
-
-	// Check label
-
-	private JDDNode checkExpressionLabel(ExpressionLabel expr) throws PrismException
-	{
-		LabelList ll;
-		JDDNode dd;
-		int i;
-
-		// treat special cases
-		if (expr.getName().equals("deadlock")) {
-			dd = model.getFixedDeadlocks();
-			JDD.Ref(dd);
-			return dd;
-		}
-		else if (expr.getName().equals("init")) {
-			dd = start;
-			JDD.Ref(dd);
-			return dd;
-		}
-		// get expression associated with label
-		ll = propertiesFile.getCombinedLabelList();
-		i = ll.getLabelIndex(expr.getName());
-		if (i == -1)
-			throw new PrismException("Unknown label \"" + expr.getName() + "\" in property");
-		// check recursively
-		return checkExpression(ll.getLabel(i));
-	}
-
-	// -----------------------------------------------------------------------------------
-	// do steady state computation
-	// -----------------------------------------------------------------------------------
-
-	// steady state computation (from initial states)
-
-	public StateProbs doSteadyState() throws PrismException
-	{
-		// bscc stuff
-		Vector vectBSCCs;
-		JDDNode notInBSCCs;
-		// mtbdd stuff
-		JDDNode start, bscc, diags, emb, tmp;
-		// other stuff
-		StateProbs probs = null, solnProbs = null;
-		double d, probBSCCs[];
-		int i, n, whichBSCC, bsccCount;
-
-		// compute bottom strongly connected components (bsccs)
-		if (bsccComp) {
-			mainLog.print("\nComputing (B)SCCs...");
-			sccComputer.computeBSCCs();
-			vectBSCCs = sccComputer.getVectBSCCs();
-			notInBSCCs = sccComputer.getNotInBSCCs();
-			n = vectBSCCs.size();
-		}
-		// unless we've been told to skip it
-		else {
-			mainLog.println("\nSkipping BSCC computation...");
-			vectBSCCs = new Vector();
-			JDD.Ref(reach);
-			vectBSCCs.add(reach);
-			notInBSCCs = JDD.Constant(0);
-			n = 1;
-		}
-
-		// get initial states of model
-		start = model.getStart();
-
-		// see how many bsccs contain initial states and, if just one, which one
-		whichBSCC = -1;
-		bsccCount = 0;
-		for (i = 0; i < n; i++) {
-			bscc = (JDDNode) vectBSCCs.elementAt(i);
-			JDD.Ref(bscc);
-			JDD.Ref(start);
-			tmp = JDD.And(bscc, start);
-			if (!tmp.equals(JDD.ZERO)) {
-				bsccCount++;
-				if (bsccCount == 1)
-					whichBSCC = i;
-			}
-			JDD.Deref(tmp);
-		}
-
-		// if all initial states are in a single bscc, it's easy...
-		JDD.Ref(notInBSCCs);
-		JDD.Ref(start);
-		tmp = JDD.And(notInBSCCs, start);
-		if (tmp.equals(JDD.ZERO) && bsccCount == 1) {
+		// derefs
+		JDD.Deref(b1);
+		JDD.Deref(b2);
+
+		return probs;
+	}
+
+	// cumulative reward
 
-			JDD.Deref(tmp);
+	protected StateProbs checkRewardCumul(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards)
+			throws PrismException
+	{
+		double time; // time
+		StateProbs rewards = null;
 
-			mainLog.println("\nInitial states all in one BSCC (so no reachability probabilities computed)");
+		// get info from inst reward
+		time = pe.getUpperBound().evaluateDouble(constantValues, null);
+		if (time < 0) {
+			throw new PrismException("Invalid time bound " + time + " in cumulative reward formula");
+		}
 
-			// get bscc
-			bscc = (JDDNode) vectBSCCs.elementAt(whichBSCC);
+		// compute rewards
 
-			// compute steady-state probabilities for the bscc
+		// a trivial case: "<=0"
+		if (time == 0) {
+			rewards = new StateProbsMTBDD(JDD.Constant(0), model);
+		} else {
+			// compute rewards
 			try {
-				solnProbs = computeSteadyStateProbs(trans, bscc);
+				rewards = computeCumulRewards(trans, trans01, stateRewards, transRewards, time);
 			} catch (PrismException e) {
-				for (i = 0; i < n; i++) {
-					JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
-				}
-				JDD.Deref(notInBSCCs);
 				throw e;
 			}
-
-			// round off probabilities
-			// solnProbs.roundOff(getPlacesToRoundBy());
 		}
 
-		// otherwise have to consider all the bsccs
-		else {
+		return rewards;
+	}
+
+	// inst reward
 
-			JDD.Deref(tmp);
+	protected StateProbs checkRewardInst(PathExpressionTemporal pe, JDDNode stateRewards, JDDNode transRewards)
+			throws PrismException
+	{
+		double time; // time
+		StateProbs sr = null, rewards = null;
+
+		// get info from inst reward
+		time = pe.getUpperBound().evaluateDouble(constantValues, null);
+		if (time < 0) {
+			throw new PrismException("Invalid bound " + time + " in instantaneous reward property");
+		}
+
+		// compute rewards
 
-			// initialise total probabilities vector
+		// a trivial case: "=0"
+		if (time == 0) {
+			JDD.Ref(stateRewards);
+			rewards = new StateProbsMTBDD(stateRewards, model);
+		} else {
+			// convert state rewards vector to appropriate type (depending on
+			// engine)
 			switch (engine) {
 			case Prism.MTBDD:
-				solnProbs = new StateProbsMTBDD(JDD.Constant(0), model);
+				JDD.Ref(stateRewards);
+				sr = new StateProbsMTBDD(stateRewards, model);
 				break;
 			case Prism.SPARSE:
-				solnProbs = new StateProbsDV(new DoubleVector((int) numStates), model);
+				sr = new StateProbsDV(stateRewards, model);
 				break;
 			case Prism.HYBRID:
-				solnProbs = new StateProbsDV(new DoubleVector((int) numStates), model);
+				sr = new StateProbsDV(stateRewards, model);
 				break;
 			}
-
-			// compute diagonals
-			JDD.Ref(trans);
-			diags = JDD.SumAbstract(trans, allDDColVars);
-
-			// compute embedded markov chain
-			JDD.Ref(trans);
-			JDD.Ref(diags);
-			emb = JDD.Apply(JDD.DIVIDE, trans, diags);
-
-			mainLog.println("\nDiagonals vector: " + JDD.GetInfoString(diags, allDDRowVars.n()));
-			mainLog.println("Embedded Markov chain: " + JDD.GetInfoString(emb, allDDRowVars.n() * 2));
-
-			// compute prob of reaching each bscc from initial state
-			probBSCCs = new double[n];
-			for (i = 0; i < n; i++) {
-
-				mainLog.println("\nComputing probability of reaching BSCC " + (i + 1));
-
-				// get bscc
-				bscc = (JDDNode) vectBSCCs.elementAt(i);
-
-				// compute probabilities
-				try {
-					probs = computeUntilProbs(emb, trans01, notInBSCCs, bscc);
-				} catch (PrismException e) {
-					JDD.Deref(diags);
-					JDD.Deref(emb);
-					for (i = 0; i < n; i++) {
-						JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
-					}
-					JDD.Deref(notInBSCCs);
-					solnProbs.clear();
-					throw e;
-				}
-
-				// round off probabilities
-				// probs.roundOff(getPlacesToRoundBy());
-
-				// sum probabilities over bdd for initial state
-				// and then divide by number of start states
-				// (we assume an equiprobable initial probability distribution
-				// over all initial states)
-				d = probs.sumOverBDD(start);
-				d /= model.getNumStartStates();
-				probBSCCs[i] = d;
-				mainLog.print("\nBSCC " + (i + 1) + " Probability: " + d + "\n");
-
-				// free vector
-				probs.clear();
-			}
-
-			// derefs
-			JDD.Deref(diags);
-			JDD.Deref(emb);
-
-			// compute steady-state for each bscc
-			for (i = 0; i < n; i++) {
-
-				mainLog.println("\nComputing steady-state probabilities for BSCC " + (i + 1));
-
-				// get bscc
-				bscc = (JDDNode) vectBSCCs.elementAt(i);
-
-				// compute steady-state probabilities for the bscc
-				try {
-					probs = computeSteadyStateProbs(trans, bscc);
-				} catch (PrismException e) {
-					for (i = 0; i < n; i++) {
-						JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
-					}
-					JDD.Deref(notInBSCCs);
-					solnProbs.clear();
-					throw e;
-				}
-
-				// round off probabilities
-				// probs.roundOff(getPlacesToRoundBy());
-
-				// print out probabilities
-				if (verbose) {
-					mainLog.print("\nBSCC " + (i + 1) + " Steady-State Probabilities: \n");
-					probs.print(mainLog);
-				}
-
-				// times by bscc reach prob, add to total
-				probs.timesConstant(probBSCCs[i]);
-				solnProbs.add(probs);
-
-				// free vector
-				probs.clear();
+			// and for the computation, we can reuse the computation for
+			// time-bounded until formulae
+			// which is nice
+			try {
+				rewards = computeBoundedUntilProbs(trans, trans01, reach, reach, time, sr);
+			} catch (PrismException e) {
+				sr.clear();
+				throw e;
 			}
+			sr.clear();
 		}
 
-		// derefs
-		for (i = 0; i < n; i++) {
-			JDD.Deref((JDDNode) vectBSCCs.elementAt(i));
-		}
-		JDD.Deref(notInBSCCs);
-
-		return solnProbs;
+		return rewards;
 	}
 
 	// -----------------------------------------------------------------------------------
@@ -2082,9 +329,6 @@ public class StochModelChecker implements ModelChecker
 			throw e;
 		}
 
-		// round off probabilities
-		// probs.roundOff(getPlacesToRoundBy());
-
 		// derefs
 		JDD.Deref(init);
 
@@ -2097,17 +341,27 @@ public class StochModelChecker implements ModelChecker
 
 	// compute probabilities for next
 
-	private StateProbs computeNextProbs(JDDNode tr, JDDNode b)
+	protected StateProbs computeNextProbs(JDDNode tr, JDDNode b)
 	{
-		JDDNode tmp;
+		JDDNode diags, emb;
 		StateProbs probs = null;
 
-		// matrix multiply: trans * b
-		JDD.Ref(b);
-		tmp = JDD.PermuteVariables(b, allDDRowVars, allDDColVars);
+		// Compute embedded Markov chain
 		JDD.Ref(tr);
-		tmp = JDD.MatrixMultiply(tr, tmp, allDDColVars, JDD.BOULDER);
-		probs = new StateProbsMTBDD(tmp, model);
+		diags = JDD.SumAbstract(tr, allDDColVars);
+		JDD.Ref(tr);
+		JDD.Ref(diags);
+		emb = JDD.Apply(JDD.DIVIDE, trans, diags);
+		mainLog.println("\nDiagonals vector: " + JDD.GetInfoString(diags, allDDRowVars.n()));
+		mainLog.println("Embedded Markov chain: " + JDD.GetInfoString(emb, allDDRowVars.n() * 2));
+
+		// And then use superclass (ProbModelChecker)
+		// to compute probabilities
+		probs = super.computeNextProbs(emb, b);
+
+		// derefs
+		JDD.Deref(diags);
+		JDD.Deref(emb);
 
 		return probs;
 	}
@@ -2127,7 +381,7 @@ public class StochModelChecker implements ModelChecker
 	// the current engine
 	// i.e. DV for sparse/hybrid, MTBDD for mtbdd
 
-	private StateProbs computeBoundedUntilProbs(JDDNode tr, JDDNode tr01, JDDNode b2, JDDNode nonabs, double time,
+	protected StateProbs computeBoundedUntilProbs(JDDNode tr, JDDNode tr01, JDDNode b2, JDDNode nonabs, double time,
 			StateProbs multProbs) throws PrismException
 	{
 		JDDNode multProbsMTBDD, probsMTBDD;
@@ -2183,188 +437,36 @@ public class StochModelChecker implements ModelChecker
 		return probs;
 	}
 
-	// compute probabilities for until (for qualitative properties)
-
-	private StateProbs computeUntilProbsQual(JDDNode tr01, JDDNode b1, JDDNode b2, double p)
-	{
-		JDDNode yes, no, maybe;
-		StateProbs probs = null;
-
-		// note: we know precomputation is enabled else this function wouldn't
-		// have been called
-
-		// compute yes/no/maybe states
-		if (b2.equals(JDD.ZERO)) {
-			yes = JDD.Constant(0);
-			JDD.Ref(reach);
-			no = reach;
-			maybe = JDD.Constant(0);
-		} else if (b1.equals(JDD.ZERO)) {
-			JDD.Ref(b2);
-			yes = b2;
-			JDD.Ref(reach);
-			JDD.Ref(b2);
-			no = JDD.And(reach, JDD.Not(b2));
-			maybe = JDD.Constant(0);
-		} else {
-			// yes
-			yes = PrismMTBDD.Prob1(tr01, allDDRowVars, allDDColVars, b1, b2);
-			if (yes.equals(reach)) {
-				no = JDD.Constant(0);
-			} else {
-				no = PrismMTBDD.Prob0(tr01, reach, allDDRowVars, allDDColVars, b1, yes);
-			}
-			// maybe
-			JDD.Ref(reach);
-			JDD.Ref(yes);
-			JDD.Ref(no);
-			maybe = JDD.And(reach, JDD.Not(JDD.Or(yes, no)));
-		}
-
-		// print out yes/no/maybe
-		mainLog.print("\nyes = " + JDD.GetNumMintermsString(yes, allDDRowVars.n()));
-		mainLog.print(", no = " + JDD.GetNumMintermsString(no, allDDRowVars.n()));
-		mainLog.print(", maybe = " + JDD.GetNumMintermsString(maybe, allDDRowVars.n()) + "\n");
-
-		// if maybe is empty, we have the probabilities already
-		if (maybe.equals(JDD.ZERO)) {
-			JDD.Ref(yes);
-			probs = new StateProbsMTBDD(yes, model);
-		}
-		// p = 0
-		else if (p == 0) {
-			// anything that's unknown but definitely > 0
-			// may as well be 1
-			JDD.Ref(yes);
-			JDD.Ref(maybe);
-			probs = new StateProbsMTBDD(JDD.Or(yes, maybe), model);
-		}
-		// p = 1
-		else {
-			// anything that's unknown but definitely < 1
-			// may as well be 0
-			JDD.Ref(yes);
-			probs = new StateProbsMTBDD(yes, model);
-		}
-
-		// derefs
-		JDD.Deref(yes);
-		JDD.Deref(no);
-		JDD.Deref(maybe);
-
-		return probs;
-	}
-
 	// compute probabilities for until (general case)
 
-	private StateProbs computeUntilProbs(JDDNode tr, JDDNode tr01, JDDNode b1, JDDNode b2) throws PrismException
+	protected StateProbs computeUntilProbs(JDDNode tr, JDDNode tr01, JDDNode b1, JDDNode b2) throws PrismException
 	{
-		JDDNode yes, no, maybe;
-		JDDNode probsMTBDD;
-		DoubleVector probsDV;
+		JDDNode diags, emb;
 		StateProbs probs = null;
 
-		// compute yes/no/maybe states
-		if (b2.equals(JDD.ZERO)) {
-			yes = JDD.Constant(0);
-			JDD.Ref(reach);
-			no = reach;
-			maybe = JDD.Constant(0);
-		} else if (b1.equals(JDD.ZERO)) {
-			JDD.Ref(b2);
-			yes = b2;
-			JDD.Ref(reach);
-			JDD.Ref(b2);
-			no = JDD.And(reach, JDD.Not(b2));
-			maybe = JDD.Constant(0);
-		} else {
-			// yes
-			if (precomp) {
-				yes = PrismMTBDD.Prob1(tr01, allDDRowVars, allDDColVars, b1, b2);
-			} else {
-				JDD.Ref(b2);
-				yes = b2;
-			}
-			// no
-			if (yes.equals(reach)) {
-				no = JDD.Constant(0);
-			} else if (precomp) {
-				no = PrismMTBDD.Prob0(tr01, reach, allDDRowVars, allDDColVars, b1, yes);
-			} else {
-				JDD.Ref(reach);
-				JDD.Ref(b1);
-				JDD.Ref(b2);
-				no = JDD.And(reach, JDD.Not(JDD.Or(b1, b2)));
-			}
-			// maybe
-			JDD.Ref(reach);
-			JDD.Ref(yes);
-			JDD.Ref(no);
-			maybe = JDD.And(reach, JDD.Not(JDD.Or(yes, no)));
-		}
-
-		// print out yes/no/maybe
-		mainLog.print("\nyes = " + JDD.GetNumMintermsString(yes, allDDRowVars.n()));
-		mainLog.print(", no = " + JDD.GetNumMintermsString(no, allDDRowVars.n()));
-		mainLog.print(", maybe = " + JDD.GetNumMintermsString(maybe, allDDRowVars.n()) + "\n");
-
-		// if maybe is empty, we have the probabilities already
-		if (maybe.equals(JDD.ZERO)) {
-			// nb: we make sure to return a vector of the appropriate type
-			// (doublevector for hybrid/sparse, mtbdd for mtbdd)
-			// this is because bounded until uses this method too
-			switch (engine) {
-			case Prism.MTBDD:
-				JDD.Ref(yes);
-				probs = new StateProbsMTBDD(yes, model);
-				break;
-			case Prism.SPARSE:
-			case Prism.HYBRID:
-				probs = new StateProbsDV(yes, model);
-				break;
-			}
-		}
-		// otherwise we compute the actual probabilities
-		else {
-			// compute probabilities
-			mainLog.println("\nComputing remaining probabilities...");
+		// Compute embedded Markov chain
+		JDD.Ref(tr);
+		diags = JDD.SumAbstract(tr, allDDColVars);
+		JDD.Ref(tr);
+		JDD.Ref(diags);
+		emb = JDD.Apply(JDD.DIVIDE, trans, diags);
+		mainLog.println("\nDiagonals vector: " + JDD.GetInfoString(diags, allDDRowVars.n()));
+		mainLog.println("Embedded Markov chain: " + JDD.GetInfoString(emb, allDDRowVars.n() * 2));
 
-			try {
-				switch (engine) {
-				case Prism.MTBDD:
-					probsMTBDD = PrismMTBDD.ProbUntil(tr, odd, allDDRowVars, allDDColVars, yes, maybe);
-					probs = new StateProbsMTBDD(probsMTBDD, model);
-					break;
-				case Prism.SPARSE:
-					probsDV = PrismSparse.ProbUntil(tr, odd, allDDRowVars, allDDColVars, yes, maybe);
-					probs = new StateProbsDV(probsDV, model);
-					break;
-				case Prism.HYBRID:
-					probsDV = PrismHybrid.ProbUntil(tr, odd, allDDRowVars, allDDColVars, yes, maybe);
-					probs = new StateProbsDV(probsDV, model);
-					break;
-				default:
-					throw new PrismException("Engine does not support this numerical method");
-				}
-			} catch (PrismException e) {
-				JDD.Deref(yes);
-				JDD.Deref(no);
-				JDD.Deref(maybe);
-				throw e;
-			}
-		}
+		// And then use superclass (ProbModelChecker)
+		// to compute probabilities
+		probs = super.computeUntilProbs(emb, tr01, b1, b2);
 
 		// derefs
-		JDD.Deref(yes);
-		JDD.Deref(no);
-		JDD.Deref(maybe);
+		JDD.Deref(diags);
+		JDD.Deref(emb);
 
 		return probs;
 	}
 
 	// compute cumulative rewards
 
-	private StateProbs computeCumulRewards(JDDNode tr, JDDNode tr01, JDDNode sr, JDDNode trr, double time)
+	protected StateProbs computeCumulRewards(JDDNode tr, JDDNode tr01, JDDNode sr, JDDNode trr, double time)
 			throws PrismException
 	{
 		JDDNode rewardsMTBDD;
@@ -2398,163 +500,35 @@ public class StochModelChecker implements ModelChecker
 
 	// compute rewards for reach reward
 
-	private StateProbs computeReachRewards(JDDNode tr, JDDNode tr01, JDDNode sr, JDDNode trr, JDDNode b)
+	protected StateProbs computeReachRewards(JDDNode tr, JDDNode tr01, JDDNode sr, JDDNode trr, JDDNode b)
 			throws PrismException
 	{
-		JDDNode inf, maybe;
-		JDDNode rewardsMTBDD;
-		DoubleVector rewardsDV;
+		JDDNode diags, emb;
 		StateProbs rewards = null;
 
-		// compute states which can't reach goal with probability 1
-		if (b.equals(JDD.ZERO)) {
-			JDD.Ref(reach);
-			inf = reach;
-			maybe = JDD.Constant(0);
-		} else if (b.equals(reach)) {
-			inf = JDD.Constant(0);
-			maybe = JDD.Constant(0);
-		} else {
-			JDDNode prob1 = PrismMTBDD.Prob1(tr01, allDDRowVars, allDDColVars, reach, b);
-			JDD.Ref(reach);
-			inf = JDD.And(reach, JDD.Not(prob1));
-			JDD.Ref(reach);
-			JDD.Ref(inf);
-			JDD.Ref(b);
-			maybe = JDD.And(reach, JDD.Not(JDD.Or(inf, b)));
-		}
-
-		// print out yes/no/maybe
-		mainLog.print("\ngoal = " + JDD.GetNumMintermsString(b, allDDRowVars.n()));
-		mainLog.print(", inf = " + JDD.GetNumMintermsString(inf, allDDRowVars.n()));
-		mainLog.print(", maybe = " + JDD.GetNumMintermsString(maybe, allDDRowVars.n()) + "\n");
-
-		// if maybe is empty, we have the rewards already
-		if (maybe.equals(JDD.ZERO)) {
-			JDD.Ref(inf);
-			rewards = new StateProbsMTBDD(JDD.ITE(inf, JDD.PlusInfinity(), JDD.Constant(0)), model);
-		}
-		// otherwise we compute the actual rewards
-		else {
-			// compute the rewards
-			mainLog.println("\nComputing remaining rewards...");
-
-			try {
-				switch (engine) {
-				case Prism.MTBDD:
-					rewardsMTBDD = PrismMTBDD.ProbReachReward(tr, sr, trr, odd, allDDRowVars, allDDColVars, b, inf,
-							maybe);
-					rewards = new StateProbsMTBDD(rewardsMTBDD, model);
-					break;
-				case Prism.SPARSE:
-					rewardsDV = PrismSparse
-							.ProbReachReward(tr, sr, trr, odd, allDDRowVars, allDDColVars, b, inf, maybe);
-					rewards = new StateProbsDV(rewardsDV, model);
-					break;
-				case Prism.HYBRID:
-					rewardsDV = PrismHybrid
-							.ProbReachReward(tr, sr, trr, odd, allDDRowVars, allDDColVars, b, inf, maybe);
-					rewards = new StateProbsDV(rewardsDV, model);
-					break;
-				default:
-					throw new PrismException("Engine does not support this numerical method");
-				}
-			} catch (PrismException e) {
-				JDD.Deref(inf);
-				JDD.Deref(maybe);
-				throw e;
-			}
-		}
-
-		// derefs
-		JDD.Deref(inf);
-		JDD.Deref(maybe);
-
-		return rewards;
-	}
-
-	// compute steady-state probabilities
-
-	// tr = the rate matrix for the whole ctmc
-	// states = the subset of reachable states (e.g. bscc) for which
-	// steady-state is to be done
-
-	private StateProbs computeSteadyStateProbs(JDDNode tr, JDDNode subset) throws PrismException
-	{
-		JDDNode trf, init;
-		long n;
-		JDDNode probsMTBDD;
-		DoubleVector probsDV;
-		StateProbs probs = null;
-
-		// work out number of states in 'subset'
-		if (tr.equals(reach)) {
-			// avoid a call to GetNumMinterms in this simple (and common) case
-			n = numStates;
-		} else {
-			n = Math.round(JDD.GetNumMinterms(subset, allDDRowVars.n()));
-		}
-
-		// special case - there is only one state in 'subset'
-		// (in fact, we need to check for this special case because the general
-		// solution work breaks)
-		if (n == 1) {
-			// answer is trivially one in the single state
-			switch (engine) {
-			case Prism.MTBDD:
-				JDD.Ref(subset);
-				return new StateProbsMTBDD(subset, model);
-			case Prism.SPARSE:
-				return new StateProbsDV(subset, model);
-			case Prism.HYBRID:
-				return new StateProbsDV(subset, model);
-			}
-		}
-
-		// filter out unwanted states from transition matrix
+		// Compute embedded Markov chain
 		JDD.Ref(tr);
-		JDD.Ref(subset);
-		trf = JDD.Apply(JDD.TIMES, tr, subset);
-		JDD.Ref(subset);
-		trf = JDD.Apply(JDD.TIMES, trf, JDD.PermuteVariables(subset, allDDRowVars, allDDColVars));
-
-		// compute initial solution (equiprobable)
-		JDD.Ref(subset);
-		init = JDD.Apply(JDD.DIVIDE, subset, JDD.Constant(n));
+		diags = JDD.SumAbstract(tr, allDDColVars);
+		JDD.Ref(tr);
+		JDD.Ref(diags);
+		emb = JDD.Apply(JDD.DIVIDE, trans, diags);
+		mainLog.println("\nDiagonals vector: " + JDD.GetInfoString(diags, allDDRowVars.n()));
+		mainLog.println("Embedded Markov chain: " + JDD.GetInfoString(emb, allDDRowVars.n() * 2));
 
-		try {
-			switch (engine) {
-			case Prism.MTBDD:
-				probsMTBDD = PrismMTBDD.StochSteadyState(trf, odd, init, allDDRowVars, allDDColVars);
-				probs = new StateProbsMTBDD(probsMTBDD, model);
-				break;
-			case Prism.SPARSE:
-				probsDV = PrismSparse.StochSteadyState(trf, odd, init, allDDRowVars, allDDColVars);
-				probs = new StateProbsDV(probsDV, model);
-				break;
-			case Prism.HYBRID:
-				probsDV = PrismHybrid.StochSteadyState(trf, odd, init, allDDRowVars, allDDColVars);
-				probs = new StateProbsDV(probsDV, model);
-				break;
-			default:
-				throw new PrismException("Engine does not support this numerical method");
-			}
-		} catch (PrismException e) {
-			JDD.Deref(trf);
-			JDD.Deref(init);
-			throw e;
-		}
+		// And then use superclass (ProbModelChecker)
+		// to compute probabilities
+		rewards = super.computeReachRewards(emb, tr01, sr, trr, b);
 
 		// derefs
-		JDD.Deref(trf);
-		JDD.Deref(init);
+		JDD.Deref(diags);
+		JDD.Deref(emb);
 
-		return probs;
+		return rewards;
 	}
 
 	// compute transient probabilities
 
-	private StateProbs computeTransientProbs(JDDNode tr, JDDNode init, double time) throws PrismException
+	protected StateProbs computeTransientProbs(JDDNode tr, JDDNode init, double time) throws PrismException
 	{
 		JDDNode probsMTBDD;
 		DoubleVector probsDV;