diff --git a/prism/src/explicit/CTMCModelChecker.java b/prism/src/explicit/CTMCModelChecker.java
index dff63be3..6436d1e2 100644
--- a/prism/src/explicit/CTMCModelChecker.java
+++ b/prism/src/explicit/CTMCModelChecker.java
@@ -384,6 +384,232 @@ public class CTMCModelChecker extends DTMCModelChecker
 		return res;
 	}
 
+	/**
+	 * Perform cumulative reward computation.
+	 * Compute, for each state of {@ctmc}, the expected rewards accumulated until {@code t}
+	 * when starting in this state and using reward structure {@code mcRewards}.
+	 * @param ctmc The CTMC
+	 * @param mcRewards The rewards
+	 * @param t Time bound
+	 */
+	@Override
+	public ModelCheckerResult computeCumulativeRewards(DTMC dtmc, MCRewards mcRewards, double t) throws PrismException
+	{
+		ModelCheckerResult res = null;
+		int i, n, iters;
+		double soln[], soln2[], tmpsoln[], sum[];
+		long timer;
+		// Fox-Glynn stuff
+		FoxGlynn fg;
+		int left, right;
+		double q, qt, acc, weights[], totalWeight;
+
+		// Optimisation: If t = 0, this is easy.
+		if (t == 0) {
+			res = new ModelCheckerResult();
+			res.soln = new double[dtmc.getNumStates()];
+			return res;
+		}
+
+		// Start backwards transient computation
+		timer = System.currentTimeMillis();
+		mainLog.println("Starting backwards cumulative rewards computation...");
+
+		CTMC ctmc = (CTMC) dtmc;
+
+		// Store num states
+		n = ctmc.getNumStates();
+
+		// Get uniformisation rate; do Fox-Glynn
+		q = ctmc.getDefaultUniformisationRate();
+		qt = q * t;
+		mainLog.println("\nUniformisation: q.t = " + q + " x " + t + " = " + qt);
+		acc = termCritParam / 8.0;
+		fg = new FoxGlynn(qt, 1e-300, 1e+300, acc);
+		left = fg.getLeftTruncationPoint();
+		right = fg.getRightTruncationPoint();
+		if (right < 0) {
+			throw new PrismException("Overflow in Fox-Glynn computation (time bound too big?)");
+		}
+		weights = fg.getWeights();
+		totalWeight = fg.getTotalWeight();
+		for (i = left; i <= right; i++) {
+			weights[i - left] /= totalWeight;
+		}
+
+		// modify the poisson probabilities to what we need for this computation
+		// first make the kth value equal to the sum of the values for 0...k
+		for (i = left+1; i <= right; i++) {
+			weights[i - left] += weights[i - 1 - left];
+		}
+		// then subtract from 1 and divide by uniformisation constant (q) to give mixed poisson probabilities
+		for (i = left; i <= right; i++) {
+			weights[i - left] = (1 - weights[i - left]) / q;
+		}
+		mainLog.println("Fox-Glynn (" + acc + "): left = " + left + ", right = " + right);
+
+		// Build (implicit) uniformised DTMC
+		dtmc = ctmc.buildImplicitUniformisedDTMC(q);
+
+		// Create solution vector(s)
+		soln = new double[n];
+		soln2 = new double[n];
+
+		// Initialise solution vectors.
+		for (i = 0; i < n; i++)
+			soln[i] = mcRewards.getStateReward(i);
+
+		// do 0th element of summation (doesn't require any matrix powers)
+		sum = new double[n];
+		if (left == 0) {
+			for (i = 0; i < n; i++)
+				sum[i] += weights[0] * soln[i];
+		} else {
+			for (i = 0; i < n; i++)
+				sum[i] += soln[i] / q;
+		}
+
+		// Start iterations
+		iters = 1;
+		while (iters <= right) {
+			// Matrix-vector multiply
+			dtmc.mvMult(soln, soln2, null, false);
+			// Swap vectors for next iter
+			tmpsoln = soln;
+			soln = soln2;
+			soln2 = tmpsoln;
+			// Add to sum
+			if (iters >= left) {
+				for (i = 0; i < n; i++)
+					sum[i] += weights[iters - left] * soln[i];
+			} else {
+				for (i = 0; i < n; i++)
+					sum[i] += soln[i] / q;
+			}
+			iters++;
+		}
+
+		// Finished backwards transient computation
+		timer = System.currentTimeMillis() - timer;
+		mainLog.print("Backwards transient cumulative rewards computation");
+		mainLog.println(" took " + iters + " iters and " + timer / 1000.0 + " seconds.");
+
+		// Return results
+		res = new ModelCheckerResult();
+		res.soln = sum;
+		res.lastSoln = soln2;
+		res.numIters = iters;
+		res.timeTaken = timer / 1000.0;
+		res.timePre = 0.0;
+		return res;
+	}
+
+	/**
+	 * Perform instantaneous reward computation.
+	 * Compute, for each state of {@ctmc}, the expected rewards at time {@code t}
+	 * when starting in this state and using reward structure {@code mcRewards}.
+	 * @param ctmc The CTMC
+	 * @param mcRewards The rewards
+	 * @param t Time bound
+	 */
+	@Override
+	public ModelCheckerResult computeInstantaneousRewards(DTMC dtmc, MCRewards mcRewards, double t) throws PrismException
+	{
+		ModelCheckerResult res = null;
+		int i, n, iters;
+		double soln[], soln2[], tmpsoln[], sum[];
+		long timer;
+		// Fox-Glynn stuff
+		FoxGlynn fg;
+		int left, right;
+		double q, qt, acc, weights[], totalWeight;
+
+		// Store num states
+		n = dtmc.getNumStates();
+
+		// Optimisation: If t = 0, this is easy.
+		if (t == 0) {
+			res = new ModelCheckerResult();
+			res.soln = new double[dtmc.getNumStates()];
+			for (i = 0; i < n; i++)
+				res.soln[i] = mcRewards.getStateReward(i);
+			return res;
+		}
+
+		// Start backwards transient computation
+		timer = System.currentTimeMillis();
+		mainLog.println("Starting backwards instantaneous rewards computation...");
+
+		CTMC ctmc = (CTMC) dtmc;
+
+		// Get uniformisation rate; do Fox-Glynn
+		q = ctmc.getDefaultUniformisationRate();
+		qt = q * t;
+		mainLog.println("\nUniformisation: q.t = " + q + " x " + t + " = " + qt);
+		acc = termCritParam / 8.0;
+		fg = new FoxGlynn(qt, 1e-300, 1e+300, acc);
+		left = fg.getLeftTruncationPoint();
+		right = fg.getRightTruncationPoint();
+		if (right < 0) {
+			throw new PrismException("Overflow in Fox-Glynn computation (time bound too big?)");
+		}
+		weights = fg.getWeights();
+		totalWeight = fg.getTotalWeight();
+		for (i = left; i <= right; i++) {
+			weights[i - left] /= totalWeight;
+		}
+
+		mainLog.println("Fox-Glynn (" + acc + "): left = " + left + ", right = " + right);
+
+		// Build (implicit) uniformised DTMC
+		dtmc = ctmc.buildImplicitUniformisedDTMC(q);
+
+		// Create solution vector(s)
+		soln = new double[n];
+		soln2 = new double[n];
+
+		// Initialise solution vectors.
+		for (i = 0; i < n; i++)
+			soln[i] = mcRewards.getStateReward(i);
+
+		// do 0th element of summation (doesn't require any matrix powers)
+		sum = new double[n];
+		if (left == 0)
+			for (i = 0; i < n; i++)
+				sum[i] += weights[0] * soln[i];
+
+		// Start iterations
+		iters = 1;
+		while (iters <= right) {
+			// Matrix-vector multiply
+			dtmc.mvMult(soln, soln2, null, false);
+			// Swap vectors for next iter
+			tmpsoln = soln;
+			soln = soln2;
+			soln2 = tmpsoln;
+			// Add to sum
+			if (iters >= left) {
+				for (i = 0; i < n; i++)
+					sum[i] += weights[iters - left] * soln[i];
+			}
+			iters++;
+		}
+
+		// Finished backwards transient computation
+		timer = System.currentTimeMillis() - timer;
+		mainLog.print("Backwards transient instantaneous rewards computation");
+		mainLog.println(" took " + iters + " iters and " + timer / 1000.0 + " seconds.");
+
+		// Return results
+		res = new ModelCheckerResult();
+		res.soln = sum;
+		res.lastSoln = soln2;
+		res.numIters = iters;
+		res.timeTaken = timer / 1000.0;
+		res.timePre = 0.0;
+		return res;
+	}
+
 	/**
 	 * Compute expected reachability rewards.
 	 * @param dtmc The CTMC
diff --git a/prism/src/explicit/DTMCModelChecker.java b/prism/src/explicit/DTMCModelChecker.java
index d803b242..f9e0f416 100644
--- a/prism/src/explicit/DTMCModelChecker.java
+++ b/prism/src/explicit/DTMCModelChecker.java
@@ -196,6 +196,12 @@ public class DTMCModelChecker extends ProbModelChecker
 			case ExpressionTemporal.R_F:
 				rewards = checkRewardReach(model, modelRewards, exprTemp);
 				break;
+			case ExpressionTemporal.R_I:
+				rewards = checkRewardInstantaneous(model, modelRewards, exprTemp);
+				break;
+			case ExpressionTemporal.R_C:
+				rewards = checkRewardCumulative(model, modelRewards, exprTemp);
+				break;
 			default:
 				throw new PrismException("Explicit engine does not yet handle the " + exprTemp.getOperatorSymbol() + " operator in the R operator");
 			}
@@ -229,6 +235,145 @@ public class DTMCModelChecker extends ProbModelChecker
 		return rewards;
 	}
 
+	/**
+	 * Compute rewards for an instantaneous reward operator.
+	 */
+	protected StateValues checkRewardInstantaneous(Model model, MCRewards modelRewards, ExpressionTemporal expr) throws PrismException
+	{
+		StateValues rewards = null;
+		ModelCheckerResult res = null;
+
+		// get time bound
+		double t = expr.getUpperBound().evaluateDouble(constantValues);
+
+		// print out some info about num states
+		// mainLog.print("\nb = " + JDD.GetNumMintermsString(b1,
+		// allDDRowVars.n()));
+
+		res = computeInstantaneousRewards((DTMC) model, modelRewards, t);
+		rewards = StateValues.createFromDoubleArray(res.soln, model);
+
+		return rewards;
+	}
+
+	/**
+	 * Compute rewards for a cumulative reward operator.
+	 */
+	protected StateValues checkRewardCumulative(Model model, MCRewards modelRewards, ExpressionTemporal expr) throws PrismException
+	{
+		StateValues rewards = null;
+		ModelCheckerResult res = null;
+
+		// get time bound
+		double t = expr.getUpperBound().evaluateDouble(constantValues);
+
+		// print out some info about num states
+		// mainLog.print("\nb = " + JDD.GetNumMintermsString(b1,
+		// allDDRowVars.n()));
+
+		res = computeCumulativeRewards((DTMC) model, modelRewards, t);
+		rewards = StateValues.createFromDoubleArray(res.soln, model);
+
+		return rewards;
+	}
+
+	public ModelCheckerResult computeInstantaneousRewards(DTMC dtmc, MCRewards mcRewards, double t) throws PrismException
+	{
+		ModelCheckerResult res = null;
+		int i, n, iters;
+		double soln[], soln2[], tmpsoln[];
+		long timer;
+		int right = (int) t;
+
+		// Store num states
+		n = dtmc.getNumStates();
+
+		// Start backwards transient computation
+		timer = System.currentTimeMillis();
+		mainLog.println("Starting backwards instantaneous rewards computation...");
+
+		// Create solution vector(s)
+		soln = new double[n];
+		soln2 = new double[n];
+
+		// Initialise solution vectors.
+		for (i = 0; i < n; i++)
+			soln[i] = mcRewards.getStateReward(i);
+
+		// Start iterations
+		for (iters = 0; iters < right; iters++) {
+			// Matrix-vector multiply
+			dtmc.mvMult(soln, soln2, null, false);
+			// Swap vectors for next iter
+			tmpsoln = soln;
+			soln = soln2;
+			soln2 = tmpsoln;
+			iters++;
+		}
+
+		// Finished backwards transient computation
+		timer = System.currentTimeMillis() - timer;
+		mainLog.print("Backwards transient instantaneous rewards computation");
+		mainLog.println(" took " + iters + " iters and " + timer / 1000.0 + " seconds.");
+
+		// Return results
+		res = new ModelCheckerResult();
+		res.soln = soln;
+		res.lastSoln = soln2;
+		res.numIters = iters;
+		res.timeTaken = timer / 1000.0;
+		res.timePre = 0.0;
+		return res;
+	}
+
+	public ModelCheckerResult computeCumulativeRewards(DTMC dtmc, MCRewards mcRewards, double t) throws PrismException
+	{
+		ModelCheckerResult res = null;
+		int i, n, iters;
+		double soln[], soln2[], tmpsoln[];
+		long timer;
+		int right = (int) t;
+
+		// Store num states
+		n = dtmc.getNumStates();
+
+		// Start backwards transient computation
+		timer = System.currentTimeMillis();
+		mainLog.println("Starting backwards cumulative rewards computation...");
+
+		// Create solution vector(s)
+		soln = new double[n];
+		soln2 = new double[n];
+
+		// Start iterations
+		for (iters = 0; iters < right; iters++) {
+			// Matrix-vector multiply plus adding rewards
+			dtmc.mvMult(soln, soln2, null, false);
+			for (i = 0; i < n; i++) {
+				soln2[i] += mcRewards.getStateReward(i);
+			}
+			// Swap vectors for next iter
+			tmpsoln = soln;
+			soln = soln2;
+			soln2 = tmpsoln;
+			iters++;
+		}
+
+		// Finished backwards transient computation
+		timer = System.currentTimeMillis() - timer;
+		mainLog.print("Backwards cumulative rewards computation");
+		mainLog.println(" took " + iters + " iters and " + timer / 1000.0 + " seconds.");
+
+		// Return results
+		res = new ModelCheckerResult();
+		res.soln = soln;
+		res.lastSoln = soln2;
+		res.numIters = iters;
+		res.timeTaken = timer / 1000.0;
+		res.timePre = 0.0;
+		return res;
+	}
+
 	/**
 	 * Compute steady-state probabilities for an S operator.
 	 */