Fox-Glynn computation: Apply accuracy check also for the 'naive' case (avoid infinite loop)

(addresses #72)

If the desired accuracy is to small, the naive weight computation can run
into an infinite loop, as floating point precision / rounding leads to
non-termination of the loop.

To partially address this, we move the accuracy check up, so that it
is also applied in the case of the 'naive' computation. This should make
it much less likely to run into this in practice.

For a full fix, we'd either need to check for non-progress in the loop
or do an analysis that the floating-point precision of double always
suffices for the remaining allowed input values.

prism/src/explicit/FoxGlynn.java

@@ -74,7 +74,12 @@ public final class FoxGlynn
 		if (q_tmax == 0.0) {
 			throw new PrismException("Overflow: TA parameter qtmax = time * maxExitRate = 0.");
 		}
-		else if (q_tmax < 400)
+
+		if (accuracy < 1e-10) {
+			throw new PrismException("Overflow: Accuracy is smaller than Fox Glynn can handle (must be at least 1e-10).");
+		}
+
+		if (q_tmax < 400)
 		{ //here naive approach should have better performance than Fox Glynn
 			final double expcoef = Math.exp(-q_tmax); //the "e^-lambda" part of p.m.f. of Poisson dist.
 			int k; //denotes that we work with event "k steps occur"
@@ -91,6 +96,9 @@ public final class FoxGlynn
 			//add further steps until you have accumulated enough
 			k = 1;
 			do {
+				// TODO: catch case where lastval gets so small that
+				// accnum never reaches desval due to rounding/floating point precision errors (infinite loop)
+
 				lastval *= q_tmax / k; // invariant: lastval = q_tmax^k / k!
 				accum += lastval;
 				w.add(lastval * expcoef);
@@ -112,9 +120,6 @@ public final class FoxGlynn
 		}
 		else
 		{ //use actual Fox Glynn for q_tmax>400
-			if (accuracy < 1e-10) {
-				throw new PrismException("Overflow: Accuracy is smaller than Fox Glynn can handle (must be at least 1e-10).");
-			}
 			final double factor = 1e+10; //factor from the paper, it has no real explanation there
 			final int m = (int) q_tmax; //mode
 			//run FINDER to get left, right and weight[m]

prism/src/prism/prism.cc

@@ -116,7 +116,14 @@ EXPORT FoxGlynnWeights fox_glynn(double q_tmax, double underflow, double overflo
 		fgw.right = -1;
 		return fgw;
 	}
-	else if (q_tmax < 400)
+
+	if (accuracy < 1e-10) {
+		printf("Overflow: Accuracy is smaller than Fox Glynn can handle (must be at least 1e-10).");
+		fgw.right = -1;
+		return fgw;
+	}
+
+	if (q_tmax < 400)
 	{ //here naive approach should have better performance than Fox Glynn
 		const double expcoef = exp(-q_tmax); //the "e^-lambda" part of p.m.f. of Poisson dist.
 		int k; //denotes that we work with event "k steps occur"
@@ -133,6 +140,9 @@ EXPORT FoxGlynnWeights fox_glynn(double q_tmax, double underflow, double overflo
 		//add further steps until you have accumulated enough
 		k = 1;
 		do {
+			// TODO: catch case where lastval gets so small that
+			// accnum never reaches desval due to rounding/floating point precision errors (infinite loop)
+
 			lastval *= q_tmax / k; // invariant: lastval = q_tmax^k / k!
 			accum += lastval;
 			w.push_back(lastval * expcoef);
@@ -154,11 +164,6 @@ EXPORT FoxGlynnWeights fox_glynn(double q_tmax, double underflow, double overflo
 	}
 	else
 	{ //use actual Fox Glynn for q_tmax>400
-		if (accuracy < 1e-10) {
-			printf("Overflow: Accuracy is smaller than Fox Glynn can handle (must be at least 1e-10).");
-			fgw.right = -1;
-			return fgw;
-		}
 		const double factor = 1e+10;
 
 		const long m = (long) q_tmax; //mode