You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2015/05/07 15:31:10 UTC
[math] Fixed wrong splitting of huge number in extended accuracy
algorithms.
Repository: commons-math
Updated Branches:
refs/heads/master 6571233ed -> e4b3ac859
Fixed wrong splitting of huge number in extended accuracy algorithms.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/e4b3ac85
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/e4b3ac85
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/e4b3ac85
Branch: refs/heads/master
Commit: e4b3ac8597a186ec8aafeda561031b95fe969d57
Parents: 6571233
Author: Luc Maisonobe <lu...@apache.org>
Authored: Thu May 7 15:23:05 2015 +0200
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Thu May 7 15:23:05 2015 +0200
----------------------------------------------------------------------
src/changes/changes.xml | 3 +
.../org/apache/commons/math4/util/FastMath.java | 30 ++++---
.../apache/commons/math4/util/MathArrays.java | 93 +++++++-------------
.../euclidean/threed/FieldVector3DTest.java | 2 +-
.../geometry/euclidean/threed/Vector3DTest.java | 2 +-
.../apache/commons/math4/util/FastMathTest.java | 21 +++++
.../commons/math4/util/MathArraysTest.java | 33 +++++++
7 files changed, 112 insertions(+), 72 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/e4b3ac85/src/changes/changes.xml
----------------------------------------------------------------------
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index d39acc0..7967a96 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -54,6 +54,9 @@ If the output is not quite correct, check for invisible trailing spaces!
</release>
<release version="4.0" date="XXXX-XX-XX" description="">
+ <action dev="luc" type="fix" issue="MATH-1223">
+ Fixed wrong splitting of huge number in extended accuracy algorithms.
+ </action>
<action dev="luc" type="fix" issue="MATH-1143">
Added helper methods to FunctionUtils for univariate and multivariate differentiable functions conversion.
</action>
http://git-wip-us.apache.org/repos/asf/commons-math/blob/e4b3ac85/src/main/java/org/apache/commons/math4/util/FastMath.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/util/FastMath.java b/src/main/java/org/apache/commons/math4/util/FastMath.java
index 24bb857..bdb420e 100644
--- a/src/main/java/org/apache/commons/math4/util/FastMath.java
+++ b/src/main/java/org/apache/commons/math4/util/FastMath.java
@@ -1634,12 +1634,10 @@ public class FastMath {
d = 1.0 / d;
}
- // split d as two 26 bits numbers
+ // split d as one 26 bits number and one 27 bits number
// beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
- final int splitFactor = 0x8000001;
- final double cd = splitFactor * d;
- final double d1High = cd - (cd - d);
- final double d1Low = d - d1High;
+ final double d1High = Double.longBitsToDouble(Double.doubleToRawLongBits(d) & ((-1L) << 27));
+ final double d1Low = d - d1High;
// prepare result
double resultHigh = 1;
@@ -1656,8 +1654,7 @@ public class FastMath {
// accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
// beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
final double tmpHigh = resultHigh * d2p;
- final double cRH = splitFactor * resultHigh;
- final double rHH = cRH - (cRH - resultHigh);
+ final double rHH = Double.longBitsToDouble(Double.doubleToRawLongBits(resultHigh) & ((-1L) << 27));
final double rHL = resultHigh - rHH;
final double tmpLow = rHL * d2pLow - (((tmpHigh - rHH * d2pHigh) - rHL * d2pHigh) - rHH * d2pLow);
resultHigh = tmpHigh;
@@ -1667,12 +1664,11 @@ public class FastMath {
// accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
// beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
final double tmpHigh = d2pHigh * d2p;
- final double cD2pH = splitFactor * d2pHigh;
+ final double cD2pH = Double.longBitsToDouble(Double.doubleToRawLongBits(d2pHigh) & ((-1L) << 27));
final double d2pHH = cD2pH - (cD2pH - d2pHigh);
final double d2pHL = d2pHigh - d2pHH;
final double tmpLow = d2pHL * d2pLow - (((tmpHigh - d2pHH * d2pHigh) - d2pHL * d2pHigh) - d2pHH * d2pLow);
- final double cTmpH = splitFactor * tmpHigh;
- d2pHigh = cTmpH - (cTmpH - tmpHigh);
+ d2pHigh = Double.longBitsToDouble(Double.doubleToRawLongBits(tmpHigh) & ((-1L) << 27));
d2pLow = d2pLow * d2p + tmpLow + (tmpHigh - d2pHigh);
d2p = d2pHigh + d2pLow;
@@ -1680,7 +1676,19 @@ public class FastMath {
}
- return resultHigh + resultLow;
+ final double result = resultHigh + resultLow;
+
+ if (Double.isNaN(result)) {
+ if (Double.isNaN(d)) {
+ return Double.NaN;
+ } else {
+ // some intermediate numbers exceeded capacity,
+ // and the low order bits became NaN (because infinity - infinity = NaN)
+ return (d < 0 && (e & 0x1) == 1) ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
+ }
+ } else {
+ return result;
+ }
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/e4b3ac85/src/main/java/org/apache/commons/math4/util/MathArrays.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/util/MathArrays.java b/src/main/java/org/apache/commons/math4/util/MathArrays.java
index ea2e9f2..f4e7f4e 100644
--- a/src/main/java/org/apache/commons/math4/util/MathArrays.java
+++ b/src/main/java/org/apache/commons/math4/util/MathArrays.java
@@ -47,8 +47,6 @@ import org.apache.commons.math4.random.Well19937c;
* @since 3.0
*/
public class MathArrays {
- /** Factor used for splitting double numbers: n = 2^27 + 1 (i.e. {@value}). */
- private static final int SPLIT_FACTOR = 0x8000001;
/**
* Private constructor.
@@ -864,15 +862,13 @@ public class MathArrays {
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;
+ final double ai = a[i];
+ final double aHigh = Double.longBitsToDouble(Double.doubleToRawLongBits(ai) & ((-1L) << 27));
+ final double aLow = ai - aHigh;
+
+ final double bi = b[i];
+ final double bHigh = Double.longBitsToDouble(Double.doubleToRawLongBits(bi) & ((-1L) << 27));
+ final double bLow = bi - bHigh;
prodHigh[i] = ai * bi;
final double prodLow = aLow * bLow - (((prodHigh[i] -
aHigh * bHigh) -
@@ -938,7 +934,6 @@ public class MathArrays {
// the code below is split in many additions/subtractions that may
// appear redundant. However, they should NOT be simplified, as they
// use IEEE754 floating point arithmetic rounding properties.
- // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
// The variable naming conventions are that xyzHigh contains the most significant
// bits of xyz and xyzLow contains its least significant bits. So theoretically
// xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
@@ -946,24 +941,20 @@ public class MathArrays {
// to hold it as long as we can, combining the high and low order bits together
// only at the end, after cancellation may have occurred on high order bits
- // split a1 and b1 as two 26 bits numbers
- final double ca1 = SPLIT_FACTOR * a1;
- final double a1High = ca1 - (ca1 - a1);
+ // split a1 and b1 as one 26 bits number and one 27 bits number
+ final double a1High = Double.longBitsToDouble(Double.doubleToRawLongBits(a1) & ((-1L) << 27));
final double a1Low = a1 - a1High;
- final double cb1 = SPLIT_FACTOR * b1;
- final double b1High = cb1 - (cb1 - b1);
+ final double b1High = Double.longBitsToDouble(Double.doubleToRawLongBits(b1) & ((-1L) << 27));
final double b1Low = b1 - b1High;
// accurate multiplication a1 * b1
final double prod1High = a1 * b1;
final double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
- // split a2 and b2 as two 26 bits numbers
- final double ca2 = SPLIT_FACTOR * a2;
- final double a2High = ca2 - (ca2 - a2);
+ // split a2 and b2 as one 26 bits number and one 27 bits number
+ final double a2High = Double.longBitsToDouble(Double.doubleToRawLongBits(a2) & ((-1L) << 27));
final double a2Low = a2 - a2High;
- final double cb2 = SPLIT_FACTOR * b2;
- final double b2High = cb2 - (cb2 - b2);
+ final double b2High = Double.longBitsToDouble(Double.doubleToRawLongBits(b2) & ((-1L) << 27));
final double b2Low = b2 - b2High;
// accurate multiplication a2 * b2
@@ -1018,7 +1009,6 @@ public class MathArrays {
// the code below is split in many additions/subtractions that may
// appear redundant. However, they should NOT be simplified, as they
// do use IEEE754 floating point arithmetic rounding properties.
- // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
// The variables naming conventions are that xyzHigh contains the most significant
// bits of xyz and xyzLow contains its least significant bits. So theoretically
// xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
@@ -1026,36 +1016,30 @@ public class MathArrays {
// to hold it as long as we can, combining the high and low order bits together
// only at the end, after cancellation may have occurred on high order bits
- // split a1 and b1 as two 26 bits numbers
- final double ca1 = SPLIT_FACTOR * a1;
- final double a1High = ca1 - (ca1 - a1);
+ // split a1 and b1 as one 26 bits number and one 27 bits number
+ final double a1High = Double.longBitsToDouble(Double.doubleToRawLongBits(a1) & ((-1L) << 27));
final double a1Low = a1 - a1High;
- final double cb1 = SPLIT_FACTOR * b1;
- final double b1High = cb1 - (cb1 - b1);
+ final double b1High = Double.longBitsToDouble(Double.doubleToRawLongBits(b1) & ((-1L) << 27));
final double b1Low = b1 - b1High;
// accurate multiplication a1 * b1
final double prod1High = a1 * b1;
final double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
- // split a2 and b2 as two 26 bits numbers
- final double ca2 = SPLIT_FACTOR * a2;
- final double a2High = ca2 - (ca2 - a2);
+ // split a2 and b2 as one 26 bits number and one 27 bits number
+ final double a2High = Double.longBitsToDouble(Double.doubleToRawLongBits(a2) & ((-1L) << 27));
final double a2Low = a2 - a2High;
- final double cb2 = SPLIT_FACTOR * b2;
- final double b2High = cb2 - (cb2 - b2);
+ final double b2High = Double.longBitsToDouble(Double.doubleToRawLongBits(b2) & ((-1L) << 27));
final double b2Low = b2 - b2High;
// accurate multiplication a2 * b2
final double prod2High = a2 * b2;
final double prod2Low = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
- // split a3 and b3 as two 26 bits numbers
- final double ca3 = SPLIT_FACTOR * a3;
- final double a3High = ca3 - (ca3 - a3);
+ // split a3 and b3 as one 26 bits number and one 27 bits number
+ final double a3High = Double.longBitsToDouble(Double.doubleToRawLongBits(a3) & ((-1L) << 27));
final double a3Low = a3 - a3High;
- final double cb3 = SPLIT_FACTOR * b3;
- final double b3High = cb3 - (cb3 - b3);
+ final double b3High = Double.longBitsToDouble(Double.doubleToRawLongBits(b3) & ((-1L) << 27));
final double b3Low = b3 - b3High;
// accurate multiplication a3 * b3
@@ -1120,7 +1104,6 @@ public class MathArrays {
// the code below is split in many additions/subtractions that may
// appear redundant. However, they should NOT be simplified, as they
// do use IEEE754 floating point arithmetic rounding properties.
- // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
// The variables naming conventions are that xyzHigh contains the most significant
// bits of xyz and xyzLow contains its least significant bits. So theoretically
// xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
@@ -1128,48 +1111,40 @@ public class MathArrays {
// to hold it as long as we can, combining the high and low order bits together
// only at the end, after cancellation may have occurred on high order bits
- // split a1 and b1 as two 26 bits numbers
- final double ca1 = SPLIT_FACTOR * a1;
- final double a1High = ca1 - (ca1 - a1);
+ // split a1 and b1 as one 26 bits number and one 27 bits number
+ final double a1High = Double.longBitsToDouble(Double.doubleToRawLongBits(a1) & ((-1L) << 27));
final double a1Low = a1 - a1High;
- final double cb1 = SPLIT_FACTOR * b1;
- final double b1High = cb1 - (cb1 - b1);
+ final double b1High = Double.longBitsToDouble(Double.doubleToRawLongBits(b1) & ((-1L) << 27));
final double b1Low = b1 - b1High;
// accurate multiplication a1 * b1
final double prod1High = a1 * b1;
final double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
- // split a2 and b2 as two 26 bits numbers
- final double ca2 = SPLIT_FACTOR * a2;
- final double a2High = ca2 - (ca2 - a2);
+ // split a2 and b2 as one 26 bits number and one 27 bits number
+ final double a2High = Double.longBitsToDouble(Double.doubleToRawLongBits(a2) & ((-1L) << 27));
final double a2Low = a2 - a2High;
- final double cb2 = SPLIT_FACTOR * b2;
- final double b2High = cb2 - (cb2 - b2);
+ final double b2High = Double.longBitsToDouble(Double.doubleToRawLongBits(b2) & ((-1L) << 27));
final double b2Low = b2 - b2High;
// accurate multiplication a2 * b2
final double prod2High = a2 * b2;
final double prod2Low = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
- // split a3 and b3 as two 26 bits numbers
- final double ca3 = SPLIT_FACTOR * a3;
- final double a3High = ca3 - (ca3 - a3);
+ // split a3 and b3 as one 26 bits number and one 27 bits number
+ final double a3High = Double.longBitsToDouble(Double.doubleToRawLongBits(a3) & ((-1L) << 27));
final double a3Low = a3 - a3High;
- final double cb3 = SPLIT_FACTOR * b3;
- final double b3High = cb3 - (cb3 - b3);
+ final double b3High = Double.longBitsToDouble(Double.doubleToRawLongBits(b3) & ((-1L) << 27));
final double b3Low = b3 - b3High;
// accurate multiplication a3 * b3
final double prod3High = a3 * b3;
final double prod3Low = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low);
- // split a4 and b4 as two 26 bits numbers
- final double ca4 = SPLIT_FACTOR * a4;
- final double a4High = ca4 - (ca4 - a4);
+ // split a4 and b4 as one 26 bits number and one 27 bits number
+ final double a4High = Double.longBitsToDouble(Double.doubleToRawLongBits(a4) & ((-1L) << 27));
final double a4Low = a4 - a4High;
- final double cb4 = SPLIT_FACTOR * b4;
- final double b4High = cb4 - (cb4 - b4);
+ final double b4High = Double.longBitsToDouble(Double.doubleToRawLongBits(b4) & ((-1L) << 27));
final double b4Low = b4 - b4High;
// accurate multiplication a4 * b4
http://git-wip-us.apache.org/repos/asf/commons-math/blob/e4b3ac85/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/FieldVector3DTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/FieldVector3DTest.java b/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/FieldVector3DTest.java
index 48e41e2..d59850f 100644
--- a/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/FieldVector3DTest.java
+++ b/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/FieldVector3DTest.java
@@ -585,7 +585,7 @@ public class FieldVector3DTest {
DerivativeStructure sNaive = u1.getX().multiply(u2.getX()).add(u1.getY().multiply(u2.getY())).add(u1.getZ().multiply(u2.getZ()));
DerivativeStructure sAccurate = FieldVector3D.dotProduct(u1, u2);
Assert.assertEquals(0.0, sNaive.getReal(), 1.0e-30);
- Assert.assertEquals(-2088690039198397.0 / 1125899906842624.0, sAccurate.getReal(), 1.0e-16);
+ Assert.assertEquals(-2088690039198397.0 / 1125899906842624.0, sAccurate.getReal(), 1.0e-15);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/e4b3ac85/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/Vector3DTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/Vector3DTest.java b/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/Vector3DTest.java
index c3bdf9f..d53a0f5 100644
--- a/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/Vector3DTest.java
+++ b/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/Vector3DTest.java
@@ -343,7 +343,7 @@ public class Vector3DTest {
double sNaive = u1.getX() * u2.getX() + u1.getY() * u2.getY() + u1.getZ() * u2.getZ();
double sAccurate = u1.dotProduct(u2);
Assert.assertEquals(0.0, sNaive, 1.0e-30);
- Assert.assertEquals(-2088690039198397.0 / 1125899906842624.0, sAccurate, 1.0e-16);
+ Assert.assertEquals(-2088690039198397.0 / 1125899906842624.0, sAccurate, 1.0e-15);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/e4b3ac85/src/test/java/org/apache/commons/math4/util/FastMathTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/util/FastMathTest.java b/src/test/java/org/apache/commons/math4/util/FastMathTest.java
index 5d36fea..1a2bca4 100644
--- a/src/test/java/org/apache/commons/math4/util/FastMathTest.java
+++ b/src/test/java/org/apache/commons/math4/util/FastMathTest.java
@@ -383,6 +383,12 @@ public class FastMathTest {
Assert.assertTrue("pow(-0.0, -3.0) should be -Inf", Double.isInfinite(FastMath.pow(-0.0, -3.0)));
+ Assert.assertEquals("pow(-0.0, Infinity) should be 0.0", 0.0, FastMath.pow(-0.0, Double.POSITIVE_INFINITY), Precision.EPSILON);
+
+ Assert.assertTrue("pow(-0.0, NaN) should be NaN", Double.isNaN(FastMath.pow(-0.0, Double.NaN)));
+
+ Assert.assertTrue("pow(-0.0, -tiny) should be Infinity", Double.isInfinite(FastMath.pow(-0.0, -Double.MIN_VALUE)));
+
Assert.assertTrue("pow(-Inf, -3.0) should be -Inf", Double.isInfinite(FastMath.pow(Double.NEGATIVE_INFINITY, 3.0)));
Assert.assertTrue("pow(-0.0, -3.5) should be Inf", Double.isInfinite(FastMath.pow(-0.0, -3.5)));
@@ -393,6 +399,16 @@ public class FastMathTest {
Assert.assertTrue("pow(-2.0, 3.5) should be NaN", Double.isNaN(FastMath.pow(-2.0, 3.5)));
+ Assert.assertTrue("pow(NaN, -Infinity) should be NaN", Double.isNaN(FastMath.pow(Double.NaN, Double.NEGATIVE_INFINITY)));
+
+ Assert.assertEquals("pow(NaN, 0.0) should be 1.0", 1.0, FastMath.pow(Double.NaN, 0.0), Precision.EPSILON);
+
+ Assert.assertEquals("pow(-Infinity, -Infinity) should be 0.0", 0.0, FastMath.pow(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY), Precision.EPSILON);
+
+ Assert.assertEquals("pow(-huge, -huge) should be 0.0", 0.0, FastMath.pow(-Double.MAX_VALUE, -Double.MAX_VALUE), Precision.EPSILON);
+
+ Assert.assertTrue("pow(-huge, huge) should be +Inf", Double.isInfinite(FastMath.pow(-Double.MAX_VALUE, Double.MAX_VALUE)));
+
// Added tests for a 100% coverage
Assert.assertTrue("pow(+Inf, NaN) should be NaN", Double.isNaN(FastMath.pow(Double.POSITIVE_INFINITY, Double.NaN)));
@@ -1188,6 +1204,11 @@ public class FastMathTest {
}
@Test
+ public void testIntPowHuge() {
+ Assert.assertTrue(Double.isInfinite(FastMath.pow(FastMath.scalb(1.0, 500), 4)));
+ }
+
+ @Test
public void testIncrementExactInt() {
int[] specialValues = new int[] {
Integer.MIN_VALUE, Integer.MIN_VALUE + 1, Integer.MIN_VALUE + 2,
http://git-wip-us.apache.org/repos/asf/commons-math/blob/e4b3ac85/src/test/java/org/apache/commons/math4/util/MathArraysTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/util/MathArraysTest.java b/src/test/java/org/apache/commons/math4/util/MathArraysTest.java
index 16e6a52..3abd0cd 100644
--- a/src/test/java/org/apache/commons/math4/util/MathArraysTest.java
+++ b/src/test/java/org/apache/commons/math4/util/MathArraysTest.java
@@ -712,6 +712,39 @@ public class MathArraysTest {
}
@Test
+ public void testLinearCombinationHuge() {
+ int scale = 971;
+ 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
+ };
+
+ double[] scaledA = new double[a.length];
+ double[] scaledB = new double[b.length];
+ for (int i = 0; i < scaledA.length; ++i) {
+ scaledA[i] = FastMath.scalb(a[i], -scale);
+ scaledB[i] = FastMath.scalb(b[i], +scale);
+ }
+ final double abSumInline = MathArrays.linearCombination(scaledA[0], scaledB[0],
+ scaledA[1], scaledB[1],
+ scaledA[2], scaledB[2]);
+ final double abSumArray = MathArrays.linearCombination(scaledA, scaledB);
+
+ Assert.assertEquals(abSumInline, abSumArray, 0);
+ Assert.assertEquals(-1.8551294182586248737720779899, abSumInline, 1.0e-15);
+
+ final double naive = scaledA[0] * scaledB[0] + scaledA[1] * scaledB[1] + scaledA[2] * scaledB[2];
+ Assert.assertTrue(FastMath.abs(naive - abSumInline) > 1.5);
+
+ }
+
+ @Test
public void testLinearCombinationInfinite() {
final double[][] a = new double[][] {
{ 1, 2, 3, 4},