MATH-1678: Avoid Canberra distance overflow

commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/distance/CanberraDistance.java

@@ -31,9 +31,16 @@ public double compute(double[] a, double[] b) {
         MathArrays.checkEqualLength(a, b);
         double sum = 0;
         for (int i = 0; i < a.length; i++) {
+            final double absA = JdkMath.abs(a[i]);
+            final double absB = JdkMath.abs(b[i]);
             final double num = JdkMath.abs(a[i] - b[i]);
-            final double denom = JdkMath.abs(a[i]) + JdkMath.abs(b[i]);
-            sum += num == 0.0 && denom == 0.0 ? 0.0 : num / denom;
+            final double denom = absA + absB;
+            if (Double.isInfinite(denom) && Double.isFinite(a[i]) && Double.isFinite(b[i])) {
+                // Rescale large finite coordinates to avoid intermediate overflow.
+                sum += JdkMath.abs(a[i] * 0.5 - b[i] * 0.5) / (absA * 0.5 + absB * 0.5);
+            } else {
+                sum += denom == 0.0 ? 0.0 : num / denom;
+            }
         }
         return sum;
     }

commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/ml/distance/CanberraDistanceTest.java

@@ -16,6 +16,9 @@
  */
 package org.apache.commons.math4.legacy.ml.distance;
 
+import java.math.BigDecimal;
+import java.math.MathContext;
+
 import org.junit.Assert;
 import org.junit.Test;
 
@@ -45,4 +48,81 @@ public void test() {
         Assert.assertEquals(expected, distance.compute(a, b), 0d);
         Assert.assertEquals(expected, distance.compute(b, a), 0d);
     }
+
+    @Test
+    public void testExtremeOppositeSigns() {
+        // Both |MAX - (-MAX)| and |MAX| + |-MAX| overflow to infinity: the true ratio is 1.
+        final double max = Double.MAX_VALUE;
+        Assert.assertEquals(1, distance.compute(new double[] {max}, new double[] {-max}), 0d);
+        Assert.assertEquals(1, distance.compute(new double[] {-max}, new double[] {max}), 0d);
+        Assert.assertEquals(1, distance.compute(new double[] {max}, new double[] {-max / 2}), 0d);
+        Assert.assertEquals(1, distance.compute(new double[] {-max / 2}, new double[] {max}), 0d);
+    }
+
+    @Test
+    public void testExtremeSameSign() {
+        // |MAX - MAX/2| / (|MAX| + |MAX/2|) == 1/3; only the denominator overflows.
+        final double max = Double.MAX_VALUE;
+        // 1/3 is not representable, so allow one ulp instead of an exact match.
+        final double expected = 1.0 / 3.0;
+        final double tol = Math.ulp(expected);
+        Assert.assertEquals(expected, distance.compute(new double[] {max}, new double[] {max / 2}), tol);
+        Assert.assertEquals(expected, distance.compute(new double[] {max / 2}, new double[] {max}), tol);
+        Assert.assertEquals(expected, distance.compute(new double[] {-max}, new double[] {-max / 2}), tol);
+        Assert.assertEquals(expected, distance.compute(new double[] {-max / 2}, new double[] {-max}), tol);
+        // Identical extreme values have distance 0.
+        Assert.assertEquals(0, distance.compute(new double[] {max}, new double[] {max}), 0d);
+        Assert.assertEquals(0, distance.compute(new double[] {-max}, new double[] {-max}), 0d);
+    }
+
+    @Test
+    public void testExtremeVector() {
+        final double max = Double.MAX_VALUE;
+        final double[] a = { max, max / 2, max };
+        final double[] b = { -max, max, max };
+        // 1 + 1/3 + 0
+        Assert.assertEquals(4.0 / 3.0, distance.compute(a, b), 1e-15);
+        Assert.assertEquals(4.0 / 3.0, distance.compute(b, a), 1e-15);
+    }
+
+    @Test
+    public void testExtremeFiniteValuesAgainstBigDecimalReference() {
+        final double max = Double.MAX_VALUE;
+        final double[] a = {max, max, 0.75 * max, max, -max, max};
+        final double[] b = {max / 2, Math.nextDown(max), Math.nextDown(0.75 * max), -max, -max / 2, max / 16};
+        for (int i = 0; i < a.length; i++) {
+            final double expected = referenceCanberra(a[i], b[i]);
+            assertWithinOneUlp(expected, distance.compute(new double[] {a[i]}, new double[] {b[i]}));
+            assertWithinOneUlp(expected, distance.compute(new double[] {b[i]}, new double[] {a[i]}));
+        }
+    }
+
+    /**
+     * Computes a high-precision reference from the exact decimal representations of the input doubles.
+     */
+    private static double referenceCanberra(double a, double b) {
+        final BigDecimal bdA = new BigDecimal(a);
+        final BigDecimal bdB = new BigDecimal(b);
+        final BigDecimal num = bdA.subtract(bdB).abs();
+        final BigDecimal den = bdA.abs().add(bdB.abs());
+        return num.divide(den, MathContext.DECIMAL128).doubleValue();
+    }
+
+    private static void assertWithinOneUlp(double expected, double actual) {
+        Assert.assertEquals("expected " + expected + " but was " + actual,
+                            expected, actual, Math.ulp(expected));
+    }
+
+    @Test
+    public void testSpecialValuesPreserved() {
+        final double min = Double.MIN_VALUE;
+        // Subnormal ratios are unaffected by the extreme-value handling.
+        Assert.assertEquals(1.0 / 3.0, distance.compute(new double[] {min}, new double[] {2 * min}), 0d);
+        Assert.assertEquals(1, distance.compute(new double[] {min}, new double[] {-min}), 0d);
+        // Signed zero: the zero/zero convention contributes 0.
+        Assert.assertEquals(0, distance.compute(new double[] {0.0}, new double[] {-0.0}), 0d);
+        // Infinite and NaN inputs keep their historical NaN result.
+        Assert.assertTrue(Double.isNaN(distance.compute(new double[] {Double.POSITIVE_INFINITY}, new double[] {1})));
+        Assert.assertTrue(Double.isNaN(distance.compute(new double[] {Double.NaN}, new double[] {0})));
+    }
 }

src/changes/changes.xml

@@ -101,6 +101,9 @@ Caveat:
         Instances not cannot be constructed that would silently overwrite entries
         due to key duplication from interger overflow.
       </action>
+      <action type="fix" issue="MATH-1678">
+        "CanberraDistance": Avoid overflow for large finite coordinates.
+      </action>
       <action dev="aherbert" type="fix" due-to="Lyin-OTTF">
         "SynchronizedDescriptiveStatistics": Synchronize all mutators
         of the internal state. Updates previously unsynchronized methods