Array version of "linearCombination".

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1154416 13f79535-47bb-0310-9956-ffa450edef68

src/main/java/org/apache/commons/math/util/MathUtils.java

@@ -2608,4 +2608,64 @@ public final class MathUtils {
 
     }
 
+    /**
+     * Compute a linear combination accurately.
+     * This method computes the sum of the products
+     * <code>a<sub>i</sub> b<sub>i</sub></code> to high accuracy.
+     * It does so by using specific multiplication and addition algorithms to
+     * preserve accuracy and reduce cancellation effects.
+     * <br/>
+     * It is based on the 2005 paper
+     * <a href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
+     * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump,
+     * and Shin'ichi Oishi published in SIAM J. Sci. Comput.
+     *
+     * @param a Factors.
+     * @param b Factors.
+     * @return <code>&Sigma;<sub>i</sub> a<sub>i</sub> b<sub>i</sub></code>.
+     */
+    public static double linearCombination(final double[] a, final double[] b) {
+        final int len = a.length;
+        if (len != b.length) {
+            throw new DimensionMismatchException(len, b.length);
+        }
+
+        final double[] prodHigh = new double[len];
+        double prodLowSum = 0;
+
+        for (int i = 0; i < len; i++) {
+            final double ai = a[i];
+            final double ca = SPLIT_FACTOR * ai;
+            final double aHigh = ca - (ca - ai);
+            final double aLow = ai - aHigh;
+
+            final double bi = b[i];
+            final double cb = SPLIT_FACTOR * bi;
+            final double bHigh = cb - (cb - bi);
+            final double bLow = bi - bHigh;
+            prodHigh[i] = ai * bi;
+            final double prodLow = aLow * bLow - (((prodHigh[i] -
+                                                    aHigh * bHigh) -
+                                                   aLow * bHigh) -
+                                                  aHigh * bLow);
+            prodLowSum += prodLow;
+        }
+
+        final int lenMinusOne = len - 1;
+        final double[] sHigh = new double[lenMinusOne];
+
+        sHigh[0] = prodHigh[0] + prodHigh[1];
+        double sPrime = sHigh[0] - prodHigh[1];
+        double sLowSum = (prodHigh[1] - (sHigh[0] - sPrime)) + (prodHigh[0] - sPrime);
+
+        for (int i = 1; i < lenMinusOne; i++) {
+            final int prev = i - 1;
+            final int next = i + 1;
+            sHigh[i] = sHigh[prev] + prodHigh[next];
+            sPrime = sHigh[i] - prodHigh[next];
+            sLowSum += (prodHigh[next] - (sHigh[i] - sPrime)) + (sHigh[prev] - sPrime);
+        }
+
+        return sHigh[lenMinusOne - 1] + (prodLowSum + sLowSum);
+    }
 }

src/test/java/org/apache/commons/math/util/MathUtilsTest.java

@@ -1841,4 +1841,25 @@ public final class MathUtilsTest {
             // Expected.
         }
     }
+
+    @Test
+    public void testLinearCombination1() {
+        final double[] a = new double[] {
+            -1321008684645961.0 / 268435456.0,
+            -5774608829631843.0 / 268435456.0,
+            -7645843051051357.0 / 8589934592.0
+        };
+        final double[] b = new double[] {
+            -5712344449280879.0 / 2097152.0,
+            -4550117129121957.0 / 2097152.0,
+            8846951984510141.0 / 131072.0
+        };
+
+        final double abSumInline = MathUtils.linearCombination(a[0], b[0],
+                                                               a[1], b[1],
+                                                               a[2], b[2]);
+        final double abSumArray = MathUtils.linearCombination(a, b);
+
+        Assert.assertEquals(abSumInline, abSumArray, 0);
+    }
 }