Repository: commons-numbers Updated Branches: refs/heads/master 7a08b2a43 -> 9f1d9ef6f
NUMBERS-30: "linear combination" ported from "Commons Math". Project: http://git-wip-us.apache.org/repos/asf/commons-numbers/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-numbers/commit/7a896524 Tree: http://git-wip-us.apache.org/repos/asf/commons-numbers/tree/7a896524 Diff: http://git-wip-us.apache.org/repos/asf/commons-numbers/diff/7a896524 Branch: refs/heads/master Commit: 7a8965240785a49a5c78c238f63582c364a98160 Parents: 7a08b2a Author: Gilles Sadowski <[email protected]> Authored: Wed May 24 17:39:47 2017 +0200 Committer: Gilles Sadowski <[email protected]> Committed: Wed May 24 17:39:47 2017 +0200 ---------------------------------------------------------------------- commons-numbers-core/pom.xml | 16 + .../commons/numbers/core/LinearCombination.java | 341 +++++++++++++++++++ .../numbers/core/LinearCombinationTest.java | 297 ++++++++++++++++ 3 files changed, 654 insertions(+) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-numbers/blob/7a896524/commons-numbers-core/pom.xml ---------------------------------------------------------------------- diff --git a/commons-numbers-core/pom.xml b/commons-numbers-core/pom.xml index 76c2512..181bab7 100644 --- a/commons-numbers-core/pom.xml +++ b/commons-numbers-core/pom.xml @@ -42,6 +42,22 @@ <numbers.parent.dir>${basedir}/..</numbers.parent.dir> </properties> + <dependencies> + <dependency> + <groupId>org.apache.commons</groupId> + <artifactId>commons-rng-simple</artifactId> + <version>1.0</version> + <scope>test</scope> + </dependency> + + <dependency> + <groupId>org.apache.commons</groupId> + <artifactId>commons-numbers-fraction</artifactId> + <version>1.0-SNAPSHOT</version> + <scope>test</scope> + </dependency> + </dependencies> + <build> <plugins> <plugin> http://git-wip-us.apache.org/repos/asf/commons-numbers/blob/7a896524/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/LinearCombination.java ---------------------------------------------------------------------- diff --git a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/LinearCombination.java b/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/LinearCombination.java new file mode 100644 index 0000000..147ff5b --- /dev/null +++ b/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/LinearCombination.java @@ -0,0 +1,341 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.numbers.core; + +/** + * Computes linear combinations 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. + * + * 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 <em>SIAM J. Sci. Comput</em>. + */ +public class LinearCombination { + /* + * Caveat: + * + * 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. + * 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 + * be represented in only one double precision number so we preserve two numbers + * 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 + */ + + /** + * @param a Factors. + * @param b Factors. + * @return \( \Sum_i a_i b_i \). + * @throws IllegalArgumentException if the sizes of the arrays are different. + */ + public static double value(double[] a, + double[] b) { + if (a.length != b.length) { + throw new IllegalArgumentException("Dimension mismatch: " + a.length + " != " + b.length); + } + + final int len = a.length; + + if (len == 1) { + // Revert to scalar multiplication. + return a[0] * b[0]; + } + + final double[] prodHigh = new double[len]; + double prodLowSum = 0; + + for (int i = 0; i < len; i++) { + final double ai = a[i]; + final double aHigh = highPart(ai); + final double aLow = ai - aHigh; + + final double bi = b[i]; + final double bHigh = highPart(bi); + final double bLow = bi - bHigh; + prodHigh[i] = ai * bi; + final double prodLow = prodLow(aLow, bLow, prodHigh[i], aHigh, bHigh); + 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; + for (int i = 1; i < lenMinusOne; i++) { + prodHighNext = prodHigh[i + 1]; + final double sHighCur = sHighPrev + prodHighNext; + sPrime = sHighCur - prodHighNext; + sLowSum += (prodHighNext - (sHighCur - sPrime)) + (sHighPrev - sPrime); + sHighPrev = sHighCur; + } + + double result = sHighPrev + (prodLowSum + sLowSum); + + if (Double.isNaN(result)) { + // either we have split infinite numbers or some coefficients were NaNs, + // just rely on the naive implementation and let IEEE754 handle this + result = 0; + for (int i = 0; i < len; ++i) { + result += a[i] * b[i]; + } + } + + return result; + } + + /** + * @param a1 First factor of the first term. + * @param b1 Second factor of the first term. + * @param a2 First factor of the second term. + * @param b2 Second factor of the second term. + * @return \( a_1 b_1 + a_2 b_2 \) + * + * @see #value(double, double, double, double, double, double) + * @see #value(double, double, double, double, double, double, double, double) + * @see #value(double[], double[]) + */ + public static double value(double a1, double b1, + double a2, double b2) { + // split a1 and b1 as one 26 bits number and one 27 bits number + final double a1High = highPart(a1); + final double a1Low = a1 - a1High; + final double b1High = highPart(b1); + final double b1Low = b1 - b1High; + + // accurate multiplication a1 * b1 + final double prod1High = a1 * b1; + final double prod1Low = prodLow(a1Low, b1Low, prod1High, a1High, b1High); + + // split a2 and b2 as one 26 bits number and one 27 bits number + final double a2High = highPart(a2); + final double a2Low = a2 - a2High; + final double b2High = highPart(b2); + final double b2Low = b2 - b2High; + + // accurate multiplication a2 * b2 + final double prod2High = a2 * b2; + final double prod2Low = prodLow(a2Low, b2Low, prod2High, a2High, b2High); + + // accurate addition a1 * b1 + a2 * b2 + final double s12High = prod1High + prod2High; + final double s12Prime = s12High - prod2High; + final double s12Low = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime); + + // final rounding, s12 may have suffered many cancellations, we try + // to recover some bits from the extra words we have saved up to now + double result = s12High + (prod1Low + prod2Low + s12Low); + + if (Double.isNaN(result)) { + // either we have split infinite numbers or some coefficients were NaNs, + // just rely on the naive implementation and let IEEE754 handle this + result = a1 * b1 + a2 * b2; + } + + return result; + } + + /** + * @param a1 First factor of the first term. + * @param b1 Second factor of the first term. + * @param a2 First factor of the second term. + * @param b2 Second factor of the second term. + * @param a3 First factor of the third term. + * @param b3 Second factor of the third term. + * @return \( a_1 b_1 + a_2 b_2 + a_3 b_3 \) + * + * @see #value(double, double, double, double) + * @see #value(double, double, double, double, double, double, double, double) + * @see #value(double[], double[]) + */ + public static double value(double a1, double b1, + double a2, double b2, + double a3, double b3) { + // split a1 and b1 as one 26 bits number and one 27 bits number + final double a1High = highPart(a1); + final double a1Low = a1 - a1High; + final double b1High = highPart(b1); + final double b1Low = b1 - b1High; + + // accurate multiplication a1 * b1 + final double prod1High = a1 * b1; + final double prod1Low = prodLow(a1Low, b1Low, prod1High, a1High, b1High); + + // split a2 and b2 as one 26 bits number and one 27 bits number + final double a2High = highPart(a2); + final double a2Low = a2 - a2High; + final double b2High = highPart(b2); + final double b2Low = b2 - b2High; + + // accurate multiplication a2 * b2 + final double prod2High = a2 * b2; + final double prod2Low = prodLow(a2Low, b2Low, prod2High, a2High, b2High); + + // split a3 and b3 as one 26 bits number and one 27 bits number + final double a3High = highPart(a3); + final double a3Low = a3 - a3High; + final double b3High = highPart(b3); + final double b3Low = b3 - b3High; + + // accurate multiplication a3 * b3 + final double prod3High = a3 * b3; + final double prod3Low = prodLow(a3Low, b3Low, prod3High, a3High, b3High); + + // accurate addition a1 * b1 + a2 * b2 + final double s12High = prod1High + prod2High; + final double s12Prime = s12High - prod2High; + final double s12Low = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime); + + // accurate addition a1 * b1 + a2 * b2 + a3 * b3 + final double s123High = s12High + prod3High; + final double s123Prime = s123High - prod3High; + final double s123Low = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime); + + // final rounding, s123 may have suffered many cancellations, we try + // to recover some bits from the extra words we have saved up to now + double result = s123High + (prod1Low + prod2Low + prod3Low + s12Low + s123Low); + + if (Double.isNaN(result)) { + // either we have split infinite numbers or some coefficients were NaNs, + // just rely on the naive implementation and let IEEE754 handle this + result = a1 * b1 + a2 * b2 + a3 * b3; + } + + return result; + } + + /** + * @param a1 First factor of the first term. + * @param b1 Second factor of the first term. + * @param a2 First factor of the second term. + * @param b2 Second factor of the second term. + * @param a3 First factor of the third term. + * @param b3 Second factor of the third term. + * @param a4 First factor of the fourth term. + * @param b4 Second factor of the fourth term. + * @return \( a_1 b_1 + a_2 b_2 + a_3 b_3 + a_4 b_4 \) + * + * @see #value(double, double, double, double) + * @see #value(double, double, double, double, double, double) + * @see #value(double[], double[]) + */ + public static double value(double a1, double b1, + double a2, double b2, + double a3, double b3, + double a4, double b4) { + // split a1 and b1 as one 26 bits number and one 27 bits number + final double a1High = highPart(a1); + final double a1Low = a1 - a1High; + final double b1High = highPart(b1); + final double b1Low = b1 - b1High; + + // accurate multiplication a1 * b1 + final double prod1High = a1 * b1; + final double prod1Low = prodLow(a1Low, b1Low, prod1High, a1High, b1High); + + // split a2 and b2 as one 26 bits number and one 27 bits number + final double a2High = highPart(a2); + final double a2Low = a2 - a2High; + final double b2High = highPart(b2); + final double b2Low = b2 - b2High; + + // accurate multiplication a2 * b2 + final double prod2High = a2 * b2; + final double prod2Low = prodLow(a2Low, b2Low, prod2High, a2High, b2High); + + // split a3 and b3 as one 26 bits number and one 27 bits number + final double a3High = highPart(a3); + final double a3Low = a3 - a3High; + final double b3High = highPart(b3); + final double b3Low = b3 - b3High; + + // accurate multiplication a3 * b3 + final double prod3High = a3 * b3; + final double prod3Low = prodLow(a3Low, b3Low, prod3High, a3High, b3High); + + // split a4 and b4 as one 26 bits number and one 27 bits number + final double a4High = highPart(a4); + final double a4Low = a4 - a4High; + final double b4High = highPart(b4); + final double b4Low = b4 - b4High; + + // accurate multiplication a4 * b4 + final double prod4High = a4 * b4; + final double prod4Low = prodLow(a4Low, b4Low, prod4High, a4High, b4High); + + // accurate addition a1 * b1 + a2 * b2 + final double s12High = prod1High + prod2High; + final double s12Prime = s12High - prod2High; + final double s12Low = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime); + + // accurate addition a1 * b1 + a2 * b2 + a3 * b3 + final double s123High = s12High + prod3High; + final double s123Prime = s123High - prod3High; + final double s123Low = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime); + + // accurate addition a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4 + final double s1234High = s123High + prod4High; + final double s1234Prime = s1234High - prod4High; + final double s1234Low = (prod4High - (s1234High - s1234Prime)) + (s123High - s1234Prime); + + // final rounding, s1234 may have suffered many cancellations, we try + // to recover some bits from the extra words we have saved up to now + double result = s1234High + (prod1Low + prod2Low + prod3Low + prod4Low + s12Low + s123Low + s1234Low); + + if (Double.isNaN(result)) { + // either we have split infinite numbers or some coefficients were NaNs, + // just rely on the naive implementation and let IEEE754 handle this + result = a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4; + } + + return result; + } + + /** + * @param value Value. + * @return the high part of the value. + */ + private static double highPart(double value) { + return Double.longBitsToDouble(Double.doubleToRawLongBits(value) & ((-1L) << 27)); + } + + /** + * @param aLow Low part of first factor. + * @param bLow Low part of second factor. + * @param prodHigh Product of the factors. + * @param aHigh High part of first factor. + * @param bHigh High part of second factor. + * @return <code>aLow * bLow - (((prodHigh - aHigh * bHigh) - aLow * bHigh) - aHigh * bLow)</code> + */ + private static double prodLow(double aLow, + double bLow, + double prodHigh, + double aHigh, + double bHigh) { + return aLow * bLow - (((prodHigh - aHigh * bHigh) - aLow * bHigh) - aHigh * bLow); + } +} http://git-wip-us.apache.org/repos/asf/commons-numbers/blob/7a896524/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/LinearCombinationTest.java ---------------------------------------------------------------------- diff --git a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/LinearCombinationTest.java b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/LinearCombinationTest.java new file mode 100644 index 0000000..d8a00c9 --- /dev/null +++ b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/LinearCombinationTest.java @@ -0,0 +1,297 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with this + * work for additional information regarding copyright ownership. The ASF + * licenses this file to You under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law + * or agreed to in writing, software distributed under the License is + * distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ +package org.apache.commons.numbers.core; + +import org.junit.Assert; +import org.junit.Test; + +import org.apache.commons.rng.UniformRandomProvider; +import org.apache.commons.rng.simple.RandomSource; +import org.apache.commons.numbers.fraction.BigFraction; + +/** + * Test cases for the {@link LinearCombination} class. + */ +public class LinearCombinationTest { + // MATH-1005 + @Test + public void testSingleElementArray() { + final double[] a = { 1.23456789 }; + final double[] b = { 98765432.1 }; + + Assert.assertEquals(a[0] * b[0], LinearCombination.value(a, b), 0d); + } + + @Test + public void testTwoSums() { + final BigFraction[] aF = new BigFraction[] { + new BigFraction(-1321008684645961L, 268435456L), + new BigFraction(-5774608829631843L, 268435456L), + new BigFraction(-7645843051051357L, 8589934592L) + }; + final BigFraction[] bF = new BigFraction[] { + new BigFraction(-5712344449280879L, 2097152L), + new BigFraction(-4550117129121957L, 2097152L), + new BigFraction(8846951984510141L, 131072L) + }; + + final int len = aF.length; + final double[] a = new double[len]; + final double[] b = new double[len]; + for (int i = 0; i < len; i++) { + a[i] = aF[i].getNumerator().doubleValue() / aF[i].getDenominator().doubleValue(); + b[i] = bF[i].getNumerator().doubleValue() / bF[i].getDenominator().doubleValue(); + } + + // Ensure "array" and "inline" implementations give the same result. + final double abSumInline = LinearCombination.value(a[0], b[0], + a[1], b[1], + a[2], b[2]); + final double abSumArray = LinearCombination.value(a, b); + Assert.assertEquals(abSumInline, abSumArray, 0); + + // Compare with arbitrary precision computation. + BigFraction result = BigFraction.ZERO; + for (int i = 0; i < a.length; i++) { + result = result.add(aF[i].multiply(bF[i])); + } + final double expected = result.doubleValue(); + Assert.assertEquals(expected, abSumInline, 1e-15); + + final double naive = a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; + Assert.assertTrue(Math.abs(naive - abSumInline) > 1.5); + } + + @Test + public void testArrayVsInline() { + final UniformRandomProvider rng = RandomSource.create(RandomSource.XOR_SHIFT_1024_S); + + double sInline; + double sArray; + final double scale = 1e17; + for (int i = 0; i < 10000; ++i) { + final double u1 = scale * rng.nextDouble(); + final double u2 = scale * rng.nextDouble(); + final double u3 = scale * rng.nextDouble(); + final double u4 = scale * rng.nextDouble(); + final double v1 = scale * rng.nextDouble(); + final double v2 = scale * rng.nextDouble(); + final double v3 = scale * rng.nextDouble(); + final double v4 = scale * rng.nextDouble(); + + // One sum. + sInline = LinearCombination.value(u1, v1, u2, v2); + sArray = LinearCombination.value(new double[] { u1, u2 }, + new double[] { v1, v2 }); + Assert.assertEquals(sInline, sArray, 0); + + // Two sums. + sInline = LinearCombination.value(u1, v1, u2, v2, u3, v3); + sArray = LinearCombination.value(new double[] { u1, u2, u3 }, + new double[] { v1, v2, v3 }); + Assert.assertEquals(sInline, sArray, 0); + + // Three sums. + sInline = LinearCombination.value(u1, v1, u2, v2, u3, v3, u4, v4); + sArray = LinearCombination.value(new double[] { u1, u2, u3, u4 }, + new double[] { v1, v2, v3, v4 }); + Assert.assertEquals(sInline, sArray, 0); + } + } + + @Test + public void testHuge() { + 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 + }; + + final int len = a.length; + final double[] scaledA = new double[len]; + final double[] scaledB = new double[len]; + for (int i = 0; i < len; ++i) { + scaledA[i] = Math.scalb(a[i], -scale); + scaledB[i] = Math.scalb(b[i], scale); + } + final double abSumInline = LinearCombination.value(scaledA[0], scaledB[0], + scaledA[1], scaledB[1], + scaledA[2], scaledB[2]); + final double abSumArray = LinearCombination.value(scaledA, scaledB); + + Assert.assertEquals(abSumInline, abSumArray, 0); + Assert.assertEquals(-1.8551294182586248737720779899, abSumInline, 1e-15); + + final double naive = scaledA[0] * scaledB[0] + scaledA[1] * scaledB[1] + scaledA[2] * scaledB[2]; + Assert.assertTrue(Math.abs(naive - abSumInline) > 1.5); + } + + @Test + public void testInfinite() { + final double[][] a = new double[][] { + { 1, 2, 3, 4 }, + { 1, Double.POSITIVE_INFINITY, 3, 4 }, + { 1, 2, Double.POSITIVE_INFINITY, 4 }, + { 1, Double.POSITIVE_INFINITY, 3, Double.NEGATIVE_INFINITY }, + { 1, 2, 3, 4 }, + { 1, 2, 3, 4 }, + { 1, 2, 3, 4 }, + { 1, 2, 3, 4 } + }; + final double[][] b = new double[][] { + { 1, -2, 3, 4 }, + { 1, -2, 3, 4 }, + { 1, -2, 3, 4 }, + { 1, -2, 3, 4 }, + { 1, Double.POSITIVE_INFINITY, 3, 4 }, + { 1, -2, Double.POSITIVE_INFINITY, 4 }, + { 1, Double.POSITIVE_INFINITY, 3, Double.NEGATIVE_INFINITY }, + { Double.NaN, -2, 3, 4 } + }; + + Assert.assertEquals(-3, + LinearCombination.value(a[0][0], b[0][0], + a[0][1], b[0][1]), + 1e-10); + Assert.assertEquals(6, + LinearCombination.value(a[0][0], b[0][0], + a[0][1], b[0][1], + a[0][2], b[0][2]), + 1e-10); + Assert.assertEquals(22, + LinearCombination.value(a[0][0], b[0][0], + a[0][1], b[0][1], + a[0][2], b[0][2], + a[0][3], b[0][3]), + 1e-10); + Assert.assertEquals(22, LinearCombination.value(a[0], b[0]), 1e-10); + + Assert.assertEquals(Double.NEGATIVE_INFINITY, + LinearCombination.value(a[1][0], b[1][0], + a[1][1], b[1][1]), + 1e-10); + Assert.assertEquals(Double.NEGATIVE_INFINITY, + LinearCombination.value(a[1][0], b[1][0], + a[1][1], b[1][1], + a[1][2], b[1][2]), + 1e-10); + Assert.assertEquals(Double.NEGATIVE_INFINITY, + LinearCombination.value(a[1][0], b[1][0], + a[1][1], b[1][1], + a[1][2], b[1][2], + a[1][3], b[1][3]), + 1e-10); + Assert.assertEquals(Double.NEGATIVE_INFINITY, LinearCombination.value(a[1], b[1]), 1e-10); + + Assert.assertEquals(-3, + LinearCombination.value(a[2][0], b[2][0], + a[2][1], b[2][1]), + 1e-10); + Assert.assertEquals(Double.POSITIVE_INFINITY, + LinearCombination.value(a[2][0], b[2][0], + a[2][1], b[2][1], + a[2][2], b[2][2]), + 1e-10); + Assert.assertEquals(Double.POSITIVE_INFINITY, + LinearCombination.value(a[2][0], b[2][0], + a[2][1], b[2][1], + a[2][2], b[2][2], + a[2][3], b[2][3]), + 1e-10); + Assert.assertEquals(Double.POSITIVE_INFINITY, LinearCombination.value(a[2], b[2]), 1e-10); + + Assert.assertEquals(Double.NEGATIVE_INFINITY, + LinearCombination.value(a[3][0], b[3][0], + a[3][1], b[3][1]), + 1e-10); + Assert.assertEquals(Double.NEGATIVE_INFINITY, + LinearCombination.value(a[3][0], b[3][0], + a[3][1], b[3][1], + a[3][2], b[3][2]), + 1e-10); + Assert.assertEquals(Double.NEGATIVE_INFINITY, + LinearCombination.value(a[3][0], b[3][0], + a[3][1], b[3][1], + a[3][2], b[3][2], + a[3][3], b[3][3]), + 1e-10); + Assert.assertEquals(Double.NEGATIVE_INFINITY, LinearCombination.value(a[3], b[3]), 1e-10); + + Assert.assertEquals(Double.POSITIVE_INFINITY, + LinearCombination.value(a[4][0], b[4][0], + a[4][1], b[4][1]), + 1e-10); + Assert.assertEquals(Double.POSITIVE_INFINITY, + LinearCombination.value(a[4][0], b[4][0], + a[4][1], b[4][1], + a[4][2], b[4][2]), + 1e-10); + Assert.assertEquals(Double.POSITIVE_INFINITY, + LinearCombination.value(a[4][0], b[4][0], + a[4][1], b[4][1], + a[4][2], b[4][2], + a[4][3], b[4][3]), + 1e-10); + Assert.assertEquals(Double.POSITIVE_INFINITY, LinearCombination.value(a[4], b[4]), 1e-10); + + Assert.assertEquals(-3, + LinearCombination.value(a[5][0], b[5][0], + a[5][1], b[5][1]), + 1e-10); + Assert.assertEquals(Double.POSITIVE_INFINITY, + LinearCombination.value(a[5][0], b[5][0], + a[5][1], b[5][1], + a[5][2], b[5][2]), + 1e-10); + Assert.assertEquals(Double.POSITIVE_INFINITY, + LinearCombination.value(a[5][0], b[5][0], + a[5][1], b[5][1], + a[5][2], b[5][2], + a[5][3], b[5][3]), + 1e-10); + Assert.assertEquals(Double.POSITIVE_INFINITY, LinearCombination.value(a[5], b[5]), 1e-10); + + Assert.assertEquals(Double.POSITIVE_INFINITY, + LinearCombination.value(a[6][0], b[6][0], + a[6][1], b[6][1]), + 1e-10); + Assert.assertEquals(Double.POSITIVE_INFINITY, + LinearCombination.value(a[6][0], b[6][0], + a[6][1], b[6][1], + a[6][2], b[6][2]), + 1e-10); + Assert.assertTrue(Double.isNaN(LinearCombination.value(a[6][0], b[6][0], + a[6][1], b[6][1], + a[6][2], b[6][2], + a[6][3], b[6][3]))); + Assert.assertTrue(Double.isNaN(LinearCombination.value(a[6], b[6]))); + + Assert.assertTrue(Double.isNaN(LinearCombination.value(a[7][0], b[7][0], + a[7][1], b[7][1]))); + Assert.assertTrue(Double.isNaN(LinearCombination.value(a[7][0], b[7][0], + a[7][1], b[7][1], + a[7][2], b[7][2]))); + Assert.assertTrue(Double.isNaN(LinearCombination.value(a[7][0], b[7][0], + a[7][1], b[7][1], + a[7][2], b[7][2], + a[7][3], b[7][3]))); + Assert.assertTrue(Double.isNaN(LinearCombination.value(a[7], b[7]))); + } +}
