diff --git a/src/main/java/org/apache/commons/math/util/MathUtils.java b/src/main/java/org/apache/commons/math/util/MathUtils.java index e54b04c5a..c7867c11e 100644 --- a/src/main/java/org/apache/commons/math/util/MathUtils.java +++ b/src/main/java/org/apache/commons/math/util/MathUtils.java @@ -2651,21 +2651,22 @@ public final class MathUtils { prodLowSum += prodLow; } + + final double prodHighCur = prodHigh[0]; + double prodHighNext = prodHigh[1]; + double sHighPrev = prodHighCur + prodHighNext; + double sPrime = sHighPrev - prodHighNext; + double sLowSum = (prodHighNext - (sHighPrev - sPrime)) + (prodHighCur - sPrime); + 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); + prodHighNext = prodHigh[i + 1]; + final double sHighCur = sHighPrev + prodHighNext; + sPrime = sHighCur - prodHighNext; + sLowSum += (prodHighNext - (sHighCur - sPrime)) + (sHighPrev - sPrime); + sHighPrev = sHighCur; } - return sHigh[lenMinusOne - 1] + (prodLowSum + sLowSum); + return sHighPrev + (prodLowSum + sLowSum); } } diff --git a/src/test/java/org/apache/commons/math/util/MathUtilsTest.java b/src/test/java/org/apache/commons/math/util/MathUtilsTest.java index cd5971291..7b6ad8e8e 100644 --- a/src/test/java/org/apache/commons/math/util/MathUtilsTest.java +++ b/src/test/java/org/apache/commons/math/util/MathUtilsTest.java @@ -29,6 +29,7 @@ import org.apache.commons.math.exception.NotFiniteNumberException; import org.apache.commons.math.exception.NullArgumentException; import org.apache.commons.math.exception.util.LocalizedFormats; import org.apache.commons.math.random.RandomDataImpl; +import org.apache.commons.math.random.Well1024a; import org.junit.Assert; import org.junit.Test; @@ -1862,4 +1863,26 @@ public final class MathUtilsTest { Assert.assertEquals(abSumInline, abSumArray, 0); } + + @Test + public void testLinearCombination2() { + // we compare accurate versus naive dot product implementations + // on regular vectors (i.e. not extreme cases like in the previous test) + Well1024a random = new Well1024a(553267312521321234l); + + for (int i = 0; i < 10000; ++i) { + final double ux = 1e17 * random.nextDouble(); + final double uy = 1e17 * random.nextDouble(); + final double uz = 1e17 * random.nextDouble(); + final double vx = 1e17 * random.nextDouble(); + final double vy = 1e17 * random.nextDouble(); + final double vz = 1e17 * random.nextDouble(); + final double sInline = MathUtils.linearCombination(ux, vx, + uy, vy, + uz, vz); + final double sArray = MathUtils.linearCombination(new double[] {ux, uy, uz}, + new double[] {vx, vy, vz}); + Assert.assertEquals(sInline, sArray, 0); + } + } }