Repository: commons-numbers Updated Branches: refs/heads/master d195a7bbb -> e4c21beac
http://git-wip-us.apache.org/repos/asf/commons-numbers/blob/e4c21bea/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 deleted file mode 100644 index 147ff5b..0000000 --- a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/LinearCombination.java +++ /dev/null @@ -1,341 +0,0 @@ -/* - * 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/e4c21bea/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/SafeNorm.java ---------------------------------------------------------------------- diff --git a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/SafeNorm.java b/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/SafeNorm.java deleted file mode 100644 index 5d077bc..0000000 --- a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/SafeNorm.java +++ /dev/null @@ -1,86 +0,0 @@ -/* - * 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 the Cartesian norm (2-norm), handling both overflow and underflow. - * Translation of the <a href="http://www.netlib.org/minpack";>minpack</a> - * "enorm" subroutine. - */ -public class SafeNorm { - /** Constant. */ - private static final double R_DWARF = 3.834e-20; - /** Constant. */ - private static final double R_GIANT = 1.304e+19; - - /** - * @param v Cartesian coordinates. - * @return the 2-norm of the vector. - */ - public static double value(double[] v) { - double s1 = 0; - double s2 = 0; - double s3 = 0; - double x1max = 0; - double x3max = 0; - double floatn = v.length; - double agiant = R_GIANT / floatn; - for (int i = 0; i < v.length; i++) { - double xabs = Math.abs(v[i]); - if (xabs < R_DWARF || xabs > agiant) { - if (xabs > R_DWARF) { - if (xabs > x1max) { - double r = x1max / xabs; - s1 = 1 + s1 * r * r; - x1max = xabs; - } else { - double r = xabs / x1max; - s1 += r * r; - } - } else { - if (xabs > x3max) { - double r = x3max / xabs; - s3 = 1 + s3 * r * r; - x3max = xabs; - } else { - if (xabs != 0) { - double r = xabs / x3max; - s3 += r * r; - } - } - } - } else { - s2 += xabs * xabs; - } - } - double norm; - if (s1 != 0) { - norm = x1max * Math.sqrt(s1 + (s2 / x1max) / x1max); - } else { - if (s2 == 0) { - norm = x3max * Math.sqrt(s3); - } else { - if (s2 >= x3max) { - norm = Math.sqrt(s2 * (1 + (x3max / s2) * (x3max * s3))); - } else { - norm = Math.sqrt(x3max * ((s2 / x3max) + (x3max * s3))); - } - } - } - return norm; - } -} http://git-wip-us.apache.org/repos/asf/commons-numbers/blob/e4c21bea/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/CosAngleTest.java ---------------------------------------------------------------------- diff --git a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/CosAngleTest.java b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/CosAngleTest.java deleted file mode 100644 index 63eeca0..0000000 --- a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/CosAngleTest.java +++ /dev/null @@ -1,77 +0,0 @@ -/* - * 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; - -/** - * Test cases for the {@link CosAngle} class. - */ -public class CosAngleTest { - @Test - public void testCosAngle2D() { - double expected; - - final double[] v1 = { 1, 0 }; - expected = 1; - Assert.assertEquals(expected, CosAngle.value(v1, v1), 0d); - - final double[] v2 = { 0, 1 }; - expected = 0; - Assert.assertEquals(expected, CosAngle.value(v1, v2), 0d); - - final double[] v3 = { 7, 7 }; - expected = Math.sqrt(2) / 2; - Assert.assertEquals(expected, CosAngle.value(v1, v3), 1e-15); - Assert.assertEquals(expected, CosAngle.value(v3, v2), 1e-15); - - final double[] v4 = { -5, 0 }; - expected = -1; - Assert.assertEquals(expected, CosAngle.value(v1, v4), 0); - - final double[] v5 = { -100, 100 }; - expected = 0; - Assert.assertEquals(expected, CosAngle.value(v3, v5), 0); - } - - @Test - public void testCosAngle3D() { - double expected; - - final double[] v1 = { 1, 1, 0 }; - expected = 1; - Assert.assertEquals(expected, CosAngle.value(v1, v1), 1e-15); - - final double[] v2 = { 1, 1, 1 }; - expected = Math.sqrt(2) / Math.sqrt(3); - Assert.assertEquals(expected, CosAngle.value(v1, v2), 1e-15); - } - - @Test - public void testCosAngleExtreme() { - double expected; - - final double tiny = 1e-200; - final double[] v1 = { tiny, tiny }; - final double big = 1e200; - final double[] v2 = { -big, -big }; - expected = -1; - Assert.assertEquals(expected, CosAngle.value(v1, v2), 1e-15); - - final double[] v3 = { big, -big }; - expected = 0; - Assert.assertEquals(expected, CosAngle.value(v1, v3), 1e-15); - } -} http://git-wip-us.apache.org/repos/asf/commons-numbers/blob/e4c21bea/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 deleted file mode 100644 index d8a00c9..0000000 --- a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/LinearCombinationTest.java +++ /dev/null @@ -1,297 +0,0 @@ -/* - * 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]))); - } -} http://git-wip-us.apache.org/repos/asf/commons-numbers/blob/e4c21bea/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SafeNormTest.java ---------------------------------------------------------------------- diff --git a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SafeNormTest.java b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SafeNormTest.java deleted file mode 100644 index 0c4e9bc..0000000 --- a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SafeNormTest.java +++ /dev/null @@ -1,62 +0,0 @@ -/* - * 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; - -/** - * Test cases for the {@link SafeNorm} class. - */ -public class SafeNormTest { - - @Test - public void testTiny() { - final double s = 1e-320; - final double[] v = new double[] { s, s }; - Assert.assertEquals(Math.sqrt(2) * s, SafeNorm.value(v), 0d); - } - - @Test - public void testBig() { - final double s = 1e300; - final double[] v = new double[] { s, s }; - Assert.assertEquals(Math.sqrt(2) * s, SafeNorm.value(v), 0d); - } - - @Test - public void testOne3D() { - final double s = 1; - final double[] v = new double[] { s, s, s }; - Assert.assertEquals(Math.sqrt(3), SafeNorm.value(v), 0d); - } - - @Test - public void testUnit3D() { - Assert.assertEquals(1, SafeNorm.value(new double[] { 1, 0, 0 }), 0d); - Assert.assertEquals(1, SafeNorm.value(new double[] { 0, 1, 0 }), 0d); - Assert.assertEquals(1, SafeNorm.value(new double[] { 0, 0, 1 }), 0d); - } - - @Test - public void testSimple() { - final double[] v = new double[] { -0.9, 8.7, -6.5, -4.3, -2.1, 0, 1.2, 3.4, -5.6, 7.8, 9.0 }; - double n = 0; - for (int i = 0; i < v.length; i++) { - n += v[i] * v[i]; - } - final double expected = Math.sqrt(n); - Assert.assertEquals(expected, SafeNorm.value(v), 0d); - } -} http://git-wip-us.apache.org/repos/asf/commons-numbers/blob/e4c21bea/pom.xml ---------------------------------------------------------------------- diff --git a/pom.xml b/pom.xml index 248a42f..7f18506 100644 --- a/pom.xml +++ b/pom.xml @@ -600,6 +600,7 @@ <module>commons-numbers-angle</module> <module>commons-numbers-gamma</module> <module>commons-numbers-combinatorics</module> + <module>commons-numbers-arrays</module> </modules> </project>
