Use BigDecimal for subtraction to resolve precision issues with double

src/org/javia/arity/BigDecimalUtils.java

@@ -0,0 +1,27 @@
+package org.javia.arity;
+
+import java.math.BigDecimal;
+
+/**
+ * Utilities for Math using BigDecimal
+ */
+public class BigDecimalUtils {
+
+    private BigDecimalUtils() {}
+
+    private static boolean isSupported(double value) {
+        return !Double.isInfinite(value) && !Double.isNaN(value);
+    }
+
+    private static BigDecimal getBigDecimalFrom(double value) {
+        return new BigDecimal(String.valueOf(value));
+    }
+
+    public static double substract(double a, double b) {
+        if (isSupported(a) && isSupported(b)) {
+            return getBigDecimalFrom(a).subtract(getBigDecimalFrom(b)).doubleValue();
+        }
+
+        return a - b;
+    }
+}

src/org/javia/arity/CompiledFunction.java

@@ -233,7 +233,7 @@ int execWithoutCheck(EvalContext context, int p) throws IsComplexException {
 
             case VM.SUB: { 
                 final double a = s[--p];
-                double res = a - (percentPC == pc-1 ? s[p] * s[p+1] : s[p+1]);
+                double res = BigDecimalUtils.substract(a, (percentPC == pc-1 ? s[p] * s[p+1] : s[p+1]));
                 if (Math.abs(res) < Math.ulp(a) * 1024) {
                     // hack for "1.1-1-.1"
                     res = 0;

src/org/javia/arity/Complex.java

@@ -150,8 +150,8 @@ public final Complex add(Complex o) {
      */
     public final Complex sub(Complex o) {
         final double ulp = Math.ulp(re);
-        re -= o.re;
-        im -= o.im;
+        re = BigDecimalUtils.substract(re, o.re);
+        im = BigDecimalUtils.substract(im, o.im);
         // hack for "1.1-1-.1"
         if (Math.abs(re) < ulp * 1024) {
             re = 0;

src/org/javia/arity/UnitTest.java

@@ -188,6 +188,8 @@ class TestEval {
         new EvalCase("imag(8.123)", 0),
         new EvalCase("im(sqrt(-1))", 1),
         new EvalCase("im(nan)", Double.NaN),
+
+        new EvalCase("56.3 - 55.7", .6),
     };
 
     private static final double ONE_SQRT2 = 0.7071067811865475; // sin(pi/4)