From a94f7bc7a5804dd8de478aa577b257c4de36ed56 Mon Sep 17 00:00:00 2001
From: Dave Parker <dave.x.parker@gmail.com>
Date: Wed, 15 Apr 2020 11:07:39 +0100
Subject: [PATCH] Make C++/Java relative error/convergence checks consistent.

In the C++ code, use |v1-v2/v1| rather than |v1-v2/v2|.
In the Java code, remove the special case where v1 and v2 are < 1e-12.

Should have minimal effect on convergence in practice,
but makes debugging easier to have engines consistent.
---
 prism/include/Measures.h        | 12 ++++++------
 prism/src/dv/dv.cc              |  7 ++++---
 prism/src/prism/PrismUtils.java | 11 +++++------
 3 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/prism/include/Measures.h b/prism/include/Measures.h
index 91d6487c..5d3055c6 100644
--- a/prism/include/Measures.h
+++ b/prism/include/Measures.h
@@ -61,7 +61,7 @@ public:
 
 	/**
 	 * Do the measurement for a single pair of values.
-	 * For relative mode, the second value is used as the divisor.
+	 * For relative mode, the first value is used as the divisor.
 	 */
 	inline void measure(double v1, double v2) {
 		double x;
@@ -70,14 +70,14 @@ public:
 		x = fabs(v2 - v1);
 		if (relative) {
 			// for relative mode: divide by second value
-			// We take the absolute value of v2 to ensure that
-			// x remains non-negative. v2 can become negative e.g.
+			// We take the absolute value of v1 to ensure that
+			// x remains non-negative. v1 can become negative e.g.
 			// during iterations with over-relaxation with large
 			// omega values.
-			// Note: if v2 is 0, then x will become +inf for x>0 and NaN for x=0, i.e., v1=v2=0
+			// Note: if v1 is 0, then x will become +inf for x>0 and NaN for x=0, i.e., v1=v2=0
 			// In the later case, the max computation below will ignore the NaN,
 			// as NaN > y is false for all y
-			x /= fabs(v2);
+			x /= fabs(v1);
 		}
 
 		// sup_norm = max { x, sup_norm }
@@ -88,7 +88,7 @@ public:
 
 	/**
 	 * Do the measurement for two value arrays of size n.
-	 * For relative mode, the values of the second array are used as the divisors.
+	 * For relative mode, the values of the first array are used as the divisors.
 	 */
 	inline void measure(const double *soln, const double *soln2, std::size_t n) {
 		for (std::size_t i = 0; i < n; i++) {
diff --git a/prism/src/dv/dv.cc b/prism/src/dv/dv.cc
index 1bf2ee2d..0a22c3b1 100644
--- a/prism/src/dv/dv.cc
+++ b/prism/src/dv/dv.cc
@@ -639,9 +639,10 @@ EXPORT bool doubles_are_close_rel(double d1, double d2, double epsilon)
 	// Compute/check error
 	d1 = fabs(d1);
 	d2 = fabs(d2);
-	// For two (near) zero values, return true, for just one, return false
-	if (d1 < epsilon_double)
-		return (d2 < epsilon_double);
+	if (d2 == 0) {
+		   // If both are zero, error is 0; otherwise +inf
+		   return (d1 == 0);
+	}
 	return (fabs(d1 - d2) / d1 < epsilon);
 }
 
diff --git a/prism/src/prism/PrismUtils.java b/prism/src/prism/PrismUtils.java
index c2aaf652..aa88a739 100644
--- a/prism/src/prism/PrismUtils.java
+++ b/prism/src/prism/PrismUtils.java
@@ -94,12 +94,11 @@ public class PrismUtils
 			return false;
 		}
 		// Compute/check error
-		d1 = Math.abs(d1);
-		d2 = Math.abs(d2);
-		// For two (near) zero values, return true, for just one, return false
-		if (d1 < epsilonDouble)
-			return (d2 < epsilonDouble);
-		return (Math.abs(d1 - d2) / d1 < epsilon);
+		if (d2 == 0) {
+			// If both are zero, error is 0; otherwise +inf
+			return d1 == 0;
+		}
+		return Math.abs((d1 - d2) / d1) < epsilon;
 	}
 
 	/**