Prevent intermediate overflow in Akima spline weighted derivatives

commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/AkimaSplineInterpolator.java

@@ -144,7 +144,23 @@ public PolynomialSplineFunction interpolate(double[] xvals,
                 final double xvM = xvals[i - 1];
                 firstDerivatives[i] = (((xvP - xv) * differences[i - 1]) + ((xv - xvM) * differences[i])) / (xvP - xvM);
             } else {
-                firstDerivatives[i] = ((wP * differences[i - 1]) + (wM * differences[i])) / (wP + wM);
+                final double weightedP = wP * differences[i - 1];
+                final double weightedM = wM * differences[i];
+                final double weightedSum = weightedP + weightedM;
+                final double weightSum = wP + wM;
+                if (Double.isFinite(weightedSum) && Double.isFinite(weightSum) &&
+                    (weightedP != 0.0 || weightedM != 0.0)) {
+                    // Preserve small contributions when normalizing would underflow a weight.
+                    firstDerivatives[i] = weightedSum / weightSum;
+                } else {
+                    // Avoid overflow or underflow of the intermediate weighted products.
+                    final double weightScale = JdkMath.max(wP, wM);
+                    final double scaledP = wP / weightScale;
+                    final double scaledM = wM / weightScale;
+                    final double scaledSum = scaledP + scaledM;
+                    firstDerivatives[i] = (scaledP / scaledSum) * differences[i - 1]
+                                        + (scaledM / scaledSum) * differences[i];
+                }
             }
         }
 

commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/analysis/interpolation/AkimaSplineInterpolatorTest.java

@@ -36,6 +36,65 @@
 import org.junit.Ignore;
 
 public class AkimaSplineInterpolatorTest {
+    @Test
+    public void testUnbalancedWeights() {
+        final double[] x = {0, 2, 4, 6, 8};
+        final double[] y = {1e-200, 0, 0, 1e200, 1e200};
+        for (boolean useModified : new boolean[] {false, true}) {
+            final UnivariateFunction interpolation =
+                new AkimaSplineInterpolator(useModified).interpolate(x, y);
+            assertEquals(-1.875e-201, interpolation.value(3), 1e-215);
+        }
+    }
+
+    @Test
+    public void testLargeFiniteValues() {
+        checkLargeFiniteValues(false);
+    }
+
+    @Test
+    public void testLargeFiniteValuesModified() {
+        checkLargeFiniteValues(true);
+    }
+
+    private void checkLargeFiniteValues(boolean useModified) {
+        final double[] x = {0, 2, 4, 6, 8};
+        final double[] y = {0, 1e200, 0, 1e200, 0};
+        final double[] expected = {0, 0.75, 1, 0.5, 0, 0.5, 1, 0.75, 0};
+        final UnivariateFunction interpolation = new AkimaSplineInterpolator(useModified).interpolate(x, y);
+        for (int i = 0; i < expected.length; i++) {
+            assertEquals(expected[i], interpolation.value(i) / 1e200, 1e-14);
+        }
+    }
+
+    @Test
+    public void testScaleInvariance() {
+        checkScaleInvariance(false);
+    }
+
+    @Test
+    public void testScaleInvarianceModified() {
+        checkScaleInvariance(true);
+    }
+
+    private void checkScaleInvariance(boolean useModified) {
+        final double[] x = {0, 2, 4, 6, 8};
+        final double[] y = {0, 1, 3, 2, 5};
+        final UnivariateFunction reference = new AkimaSplineInterpolator(useModified).interpolate(x, y);
+        for (double scale : new double[] {1e200, -1e200, 1e-200}) {
+            final double[] scaledY = new double[y.length];
+            for (int i = 0; i < y.length; i++) {
+                scaledY[i] = y[i] * scale;
+            }
+            final UnivariateFunction interpolation =
+                new AkimaSplineInterpolator(useModified).interpolate(x, scaledY);
+            for (int i = 0; i <= 16; i++) {
+                final double point = i * 0.5;
+                assertEquals(reference.value(point), interpolation.value(point) / scale, 1e-14);
+            }
+        }
+    }
+
     @Test
     public void testIllegalArguments() {
         // Data set arrays of different size.

src/changes/changes.xml

@@ -96,6 +96,10 @@ Caveat:
  to support the whole codebase (it was one of the main reasons for
  creating more focused components).
 ">
+      <action type="fix" due-to="JackieTien97">
+        "AkimaSplineInterpolator": Normalize derivative weights to avoid
+        overflow and underflow in intermediate weighted products.
+      </action>
       <action dev="aherbert" type="fix" due-to="simpleqt">
         "MathArrays": Fix int overflow in the int[] L1 and Linf distance.
       </action>