Use BigDecimal for subtraction to resolve precision issues with double

arity/src/main/java/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;
+    }
+}

arity/src/main/java/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;

arity/src/main/java/org/javia/arity/Complex.java

@@ -176,8 +176,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;

arity/src/test/java/org/javia/arity/Eval.java

@@ -176,6 +176,8 @@ private static class EvalCase {
           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)