Repository: Hipparchus-Math/hipparchus (java)

Issue:
Add Møller's 2Sum algorithm
Add Møller's 2Sum algorithm based on discussion [here](https://gitlab.orekit.org/orekit/orekit/-/merge_requests/149).

Pull request:
Added Møller's 2Sum algorithm to MathUtils.
Closes #116.

Patch:
diff --git a/hipparchus-core/src/changes/changes.xml b/hipparchus-core/src/changes/changes.xml
index c745a2c48..731bcd23a 100644
--- a/hipparchus-core/src/changes/changes.xml
+++ b/hipparchus-core/src/changes/changes.xml
@@ -50,6 +50,9 @@ If the output is not quite correct, check for invisible trailing spaces!
   </properties>
   <body>
     <release version="X.Y" date="TBD" description="TBD">
+      <action dev="andrewsgoetz" type="add" issue="issues/116">
+        Added Møller's 2Sum algorithm to MathUtils.
+      </action>
       <action dev="luc" type="fix" issue="issues/111">
         Allow selecting customized threshold for tests in OrderedComplexEigenDecomposition.
       </action>
diff --git a/hipparchus-core/src/main/java/org/hipparchus/util/MathUtils.java b/hipparchus-core/src/main/java/org/hipparchus/util/MathUtils.java
index 8ddfa1185..6433348d3 100644
--- a/hipparchus-core/src/main/java/org/hipparchus/util/MathUtils.java
+++ b/hipparchus-core/src/main/java/org/hipparchus/util/MathUtils.java
@@ -45,8 +45,6 @@ public final class MathUtils {
     /** \(\pi^2\) */
     public static final double PI_SQUARED = FastMath.PI * FastMath.PI;
 
-    /** \(\pi/2\). */
-    public static final double SEMI_PI = 0.5 * FastMath.PI;
 
     /**
      * Class contains only static methods.
@@ -370,4 +368,139 @@ public final class MathUtils {
         }
     }
 
+    /**
+     * Sums {@code a} and {@code b} using Møller's 2Sum algorithm.
+     * <p>
+     * References:
+     * <ul>
+     * <li>Møller, Ole. "Quasi double-precision in floating point addition." BIT
+     * 5, 37–50 (1965).</li>
+     * <li>Shewchuk, Richard J. "Adaptive Precision Floating-Point Arithmetic
+     * and Fast Robust Geometric Predicates." Discrete & Computational Geometry
+     * 18, 305–363 (1997).</li>
+     * <li><a href=
+     * "https://en.wikipedia.org/wiki/2Sum">https://en.wikipedia.org/wiki/2Sum</a></li>
+     * </ul>
+     * @param a first summand
+     * @param b second summand
+     * @return sum and residual error in the sum
+     */
+    public static SumAndResidual twoSum(final double a, final double b) {
+        final double s = a + b;
+        final double aPrime = s - b;
+        final double bPrime = s - aPrime;
+        final double deltaA = a - aPrime;
+        final double deltaB = b - bPrime;
+        final double t = deltaA + deltaB;
+        return new SumAndResidual(s, t);
+    }
+
+    /**
+     * Sums {@code a} and {@code b} using Møller's 2Sum algorithm.
+     * <p>
+     * References:
+     * <ul>
+     * <li>Møller, Ole. "Quasi double-precision in floating point addition." BIT
+     * 5, 37–50 (1965).</li>
+     * <li>Shewchuk, Richard J. "Adaptive Precision Floating-Point Arithmetic
+     * and Fast Robust Geometric Predicates." Discrete & Computational Geometry
+     * 18, 305–363 (1997).</li>
+     * <li><a href=
+     * "https://en.wikipedia.org/wiki/2Sum">https://en.wikipedia.org/wiki/2Sum</a></li>
+     * </ul>
+     * @param <T> field element type
+     * @param a first summand
+     * @param b second summand
+     * @return sum and residual error in the sum
+     */
+    public static <T extends RealFieldElement<T>> FieldSumAndResidual<T> twoSum(final T a, final T b) {
+        final T s = a.add(b);
+        final T aPrime = s.subtract(b);
+        final T bPrime = s.subtract(aPrime);
+        final T deltaA = a.subtract(aPrime);
+        final T deltaB = b.subtract(bPrime);
+        final T t = deltaA.add(deltaB);
+        return new FieldSumAndResidual<>(s, t);
+    }
+
+    /**
+     * Result class for {@link MathUtils#twoSum(double, double)} containing the
+     * sum and the residual error in the sum.
+     */
+    public static final class SumAndResidual {
+
+        /** Sum. */
+        private final double sum;
+        /** Residual error in the sum. */
+        private final double residual;
+
+        /**
+         * Constructs a {@link SumAndResidual} instance.
+         * @param sum sum
+         * @param residual residual error in the sum
+         */
+        private SumAndResidual(final double sum, final double residual) {
+            this.sum = sum;
+            this.residual = residual;
+        }
+
+        /**
+         * Returns the sum.
+         * @return sum
+         */
+        public double getSum() {
+            return sum;
+        }
+
+        /**
+         * Returns the residual error in the sum.
+         * @return residual error in the sum
+         */
+        public double getResidual() {
+            return residual;
+        }
+
+    }
+
+    /**
+     * Result class for
+     * {@link MathUtils#twoSum(RealFieldElement, RealFieldElement)} containing
+     * the sum and the residual error in the sum.
+     * @param <T> field element type
+     */
+    public static final class FieldSumAndResidual<T extends RealFieldElement<T>> {
+
+        /** Sum. */
+        private final T sum;
+        /** Residual error in the sum. */
+        private final T residual;
+
+        /**
+         * Constructs a {@link FieldSumAndResidual} instance.
+         * @param sum sum
+         * @param residual residual error in the sum
+         */
+        private FieldSumAndResidual(final T sum, final T residual) {
+            this.sum = sum;
+            this.residual = residual;
+        }
+
+        /**
+         * Returns the sum.
+         * @return sum
+         */
+        public T getSum() {
+            return sum;
+        }
+
+        /**
+         * Returns the residual error in the sum.
+         * @return residual error in the sum
+         */
+        public T getResidual() {
+            return residual;
+        }
+
+    }
+
 }


Test patch:
diff --git a/hipparchus-core/src/test/java/org/hipparchus/util/MathUtilsTest.java b/hipparchus-core/src/test/java/org/hipparchus/util/MathUtilsTest.java
index 80185df91..aee3051b2 100644
--- a/hipparchus-core/src/test/java/org/hipparchus/util/MathUtilsTest.java
+++ b/hipparchus-core/src/test/java/org/hipparchus/util/MathUtilsTest.java
@@ -20,6 +20,8 @@ import org.hipparchus.exception.MathIllegalArgumentException;
 import org.hipparchus.exception.MathRuntimeException;
 import org.hipparchus.exception.NullArgumentException;
 import org.hipparchus.random.RandomDataGenerator;
+import org.hipparchus.util.MathUtils.FieldSumAndResidual;
+import org.hipparchus.util.MathUtils.SumAndResidual;
 import org.junit.Assert;
 import org.junit.Test;
 
@@ -367,4 +369,63 @@ public final class MathUtilsTest {
     public void testCopySignByte2() {
         MathUtils.copySign(Byte.MIN_VALUE, (byte) 1);
     }
+
+    /**
+     * Tests {@link TwoSum#twoSum(double, double)}.
+     */
+    @Test
+    public void testTwoSum() {
+        //      sum = 0.30000000000000004
+        // residual = -2.7755575615628914E-17
+        final double a1 = 0.1;
+        final double b1 = 0.2;
+        final SumAndResidual result1 = MathUtils.twoSum(a1, b1);
+        Assert.assertEquals(a1 + b1, result1.getSum(), 0.);
+        Assert.assertEquals(a1 + b1, result1.getSum() + result1.getResidual(), 0.);
+        Assert.assertNotEquals(0., result1.getResidual(), 0.);
+
+        //      sum = -1580.3205849419005
+        // residual = -1.1368683772161603E-13
+        final double a2 = -615.7212034581913;
+        final double b2 = -964.5993814837093;
+        final SumAndResidual result2 = MathUtils.twoSum(a2, b2);
+        Assert.assertEquals(a2 + b2, result2.getSum(), 0.);
+        Assert.assertEquals(a2 + b2, result2.getSum() + result2.getResidual(), 0.);
+        Assert.assertNotEquals(0., result2.getResidual(), 0.);
+
+        //      sum = 251.8625825973395
+        // residual = 1.4210854715202004E-14
+        final double a3 = 60.348375484313706;
+        final double b3 = 191.5142071130258;
+        final SumAndResidual result3 = MathUtils.twoSum(a3, b3);
+        Assert.assertEquals(a3 + b3, result3.getSum(), 0.);
+        Assert.assertEquals(a3 + b3, result3.getSum() + result3.getResidual(), 0.);
+        Assert.assertNotEquals(0., result3.getResidual(), 0.);
+
+        //      sum = 622.8319023175123
+        // residual = -4.3315557304163255E-14
+        final double a4 = 622.8314146170453;
+        final double b4 = 0.0004877004669900762;
+        final SumAndResidual result4 = MathUtils.twoSum(a4, b4);
+        Assert.assertEquals(a4 + b4, result4.getSum(), 0.);
+        Assert.assertEquals(a4 + b4, result4.getSum() + result4.getResidual(), 0.);
+        Assert.assertNotEquals(0., result4.getResidual(), 0.);
+    }
+
+    /**
+     * Tests {@link TwoSum#twoSum(RealFieldElement, RealFieldElement)}.
+     */
+    @Test
+    public void testTwoSumField() {
+        final Tuple a = new Tuple(0.1, -615.7212034581913, 60.348375484313706, 622.8314146170453);
+        final Tuple b = new Tuple(0.2, -964.5993814837093, 191.5142071130258, 0.0004877004669900762);
+        final FieldSumAndResidual<Tuple> result = MathUtils.twoSum(a, b);
+        for (int i = 0; i < a.getDimension(); ++i) {
+            Assert.assertEquals(a.getComponent(i) + b.getComponent(i), result.getSum().getComponent(i), 0.);
+            Assert.assertEquals(a.getComponent(i) + b.getComponent(i),
+                    result.getSum().getComponent(i) + result.getResidual().getComponent(i), 0.);
+            Assert.assertNotEquals(0., result.getResidual().getComponent(i), 0.);
+        }
+    }
+
 }