From a171143586c077e279d60c86f505fb87e3dba7a3 Mon Sep 17 00:00:00 2001
From: woheller69 <wolfgang.heller@web.de>
Date: Sun, 29 Jan 2023 18:21:17 +0100
Subject: [PATCH] Use BigDecimal for subtraction to resolve precision issues
 with double

Patch from:
https://github.com/piusvelte/arithmetic/commit/5fb74edd2ca32401a93dcbae215a155c1693bd27
---
 .../java/org/javia/arity/BigDecimalUtils.java | 27 +++++++++++++++++++
 .../org/javia/arity/CompiledFunction.java     |  2 +-
 .../main/java/org/javia/arity/Complex.java    |  4 +--
 arity/src/test/java/org/javia/arity/Eval.java |  2 ++
 4 files changed, 32 insertions(+), 3 deletions(-)
 create mode 100644 arity/src/main/java/org/javia/arity/BigDecimalUtils.java

diff --git a/arity/src/main/java/org/javia/arity/BigDecimalUtils.java b/arity/src/main/java/org/javia/arity/BigDecimalUtils.java
new file mode 100644
index 0000000..795ccc8
--- /dev/null
+++ b/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;
+    }
+}
diff --git a/arity/src/main/java/org/javia/arity/CompiledFunction.java b/arity/src/main/java/org/javia/arity/CompiledFunction.java
index c2ae621..ced3328 100644
--- a/arity/src/main/java/org/javia/arity/CompiledFunction.java
+++ 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;
diff --git a/arity/src/main/java/org/javia/arity/Complex.java b/arity/src/main/java/org/javia/arity/Complex.java
index 8140972..f89bdae 100644
--- a/arity/src/main/java/org/javia/arity/Complex.java
+++ b/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;
diff --git a/arity/src/test/java/org/javia/arity/Eval.java b/arity/src/test/java/org/javia/arity/Eval.java
index 47e91aa..40376c2 100644
--- a/arity/src/test/java/org/javia/arity/Eval.java
+++ b/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)