You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2021/06/20 02:30:25 UTC
[commons-numbers] 01/02: NUMBERS-163: replacing LinearCombination
and Summation with Sum class
This is an automated email from the ASF dual-hosted git repository.
erans pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit 06a613d29c950acff30feb06ac95d54c0d5151ea
Author: Matt Juntunen <ma...@apache.org>
AuthorDate: Sat Jun 19 08:33:35 2021 -0400
NUMBERS-163: replacing LinearCombination and Summation with Sum class
---
.../org/apache/commons/numbers/angle/CosAngle.java | 5 +-
.../commons/numbers/core/LinearCombination.java | 222 --------
.../java/org/apache/commons/numbers/core/Norm.java | 16 +-
.../java/org/apache/commons/numbers/core/Sum.java | 258 ++++++++++
.../org/apache/commons/numbers/core/Summation.java | 179 -------
.../numbers/core/LinearCombinationTest.java | 416 ---------------
.../org/apache/commons/numbers/core/SumTest.java | 568 +++++++++++++++++++++
.../apache/commons/numbers/core/SummationTest.java | 283 ----------
.../examples/jmh/core/EuclideanNormAlgorithms.java | 4 +-
.../jmh/core/LinearCombinationPerformance.java | 15 +-
.../numbers/examples/jmh/core/NormPerformance.java | 6 +-
.../numbers/examples/jmh/core/SumPerformance.java | 183 +++++++
12 files changed, 1033 insertions(+), 1122 deletions(-)
diff --git a/commons-numbers-angle/src/main/java/org/apache/commons/numbers/angle/CosAngle.java b/commons-numbers-angle/src/main/java/org/apache/commons/numbers/angle/CosAngle.java
index d9d7dd4..fcbabba 100644
--- a/commons-numbers-angle/src/main/java/org/apache/commons/numbers/angle/CosAngle.java
+++ b/commons-numbers-angle/src/main/java/org/apache/commons/numbers/angle/CosAngle.java
@@ -16,8 +16,8 @@
*/
package org.apache.commons.numbers.angle;
-import org.apache.commons.numbers.core.LinearCombination;
import org.apache.commons.numbers.core.Norm;
+import org.apache.commons.numbers.core.Sum;
/**
* Computes the cosine of the angle between two vectors.
@@ -39,7 +39,6 @@ public final class CosAngle {
*/
public static double value(double[] v1,
double[] v2) {
- return LinearCombination.value(v1, v2) / Norm.L2.of(v1) / Norm.L2.of(v2);
+ return Sum.ofProducts(v1, v2).getAsDouble() / Norm.L2.of(v1) / Norm.L2.of(v2);
}
}
-
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 09cc821..0000000
--- a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/LinearCombination.java
+++ /dev/null
@@ -1,222 +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.
- *
- * <p>It is based on the 2005 paper
- * <a href="https://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 final class LinearCombination {
- /*
- * Caveat:
- *
- * The code below uses many additions/subtractions that may
- * appear redundant. However, they should NOT be simplified, as they
- * do use IEEE754 floating point arithmetic rounding properties.
- *
- * Algorithms are based on computing the product or sum of two values x and y in
- * extended precision. The standard result is stored using a double (high part z) and
- * the round-off error (or low part zz) is stored in a second double, e.g:
- * x * y = (z, zz); z + zz = x * y
- * x + y = (z, zz); z + zz = x + y
- *
- * To sum multiple (z, zz) results ideally the parts are sorted in order of
- * non-decreasing magnitude and summed. This is exact if each number's most significant
- * bit is below the least significant bit of the next (i.e. does not
- * overlap). Creating non-overlapping parts requires a rebalancing
- * of adjacent pairs using a summation z + zz = (z1, zz1) iteratively through the parts
- * (see Shewchuk (1997) Grow-Expansion and Expansion-Sum [1]).
- *
- * In this class the sum of the low parts in computed separately from the sum of the
- * high parts for an approximate 2-fold increase in precision in the event of cancellation
- * (sum positives and negatives to a result of much smaller magnitude than the parts).
- * Uses the dot2s algorithm of Ogita to avoid allocation of an array to store intermediates.
- *
- * [1] Shewchuk (1997): Arbitrary Precision Floating-Point Arithmetic
- * http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps
- */
-
- /** Private constructor. */
- private LinearCombination() {
- // intentionally empty.
- }
-
- /**
- * @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);
- }
-
- // Implement dot2s (Algorithm 5.4) from Ogita et al (2005).
- final int len = a.length;
-
- // p is the standard scalar product sum.
- // s is the sum of round-off parts.
- double p = a[0] * b[0];
- double s = ExtendedPrecision.productLow(a[0], b[0], p);
-
- // Remaining split products added to the current sum and round-off sum.
- for (int i = 1; i < len; i++) {
- final double h = a[i] * b[i];
- final double r = ExtendedPrecision.productLow(a[i], b[i], h);
-
- final double x = p + h;
- // s_i = s_(i-1) + (q_i + r_i)
- s += ExtendedPrecision.twoSumLow(p, h, x) + r;
- p = x;
- }
-
- return getSum(p, p + s);
- }
-
- /**
- * @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) {
- // p/pn are the standard scalar product old/new sum.
- // s is the sum of round-off parts.
- final double p = a1 * b1;
- double s = ExtendedPrecision.productLow(a1, b1, p);
- final double h = a2 * b2;
- final double r = ExtendedPrecision.productLow(a2, b2, h);
- final double pn = p + h;
- s += ExtendedPrecision.twoSumLow(p, h, pn) + r;
-
- // Final summation
- return getSum(pn, pn + s);
- }
-
- /**
- * @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) {
- // p/q are the standard scalar product old/new sum (alternating).
- // s is the sum of round-off parts.
- // pn is the final scalar product sum.
- final double p = a1 * b1;
- double s = ExtendedPrecision.productLow(a1, b1, p);
- double h = a2 * b2;
- double r = ExtendedPrecision.productLow(a2, b2, h);
- final double q = p + h;
- s += r + ExtendedPrecision.twoSumLow(p, h, q);
- h = a3 * b3;
- r = ExtendedPrecision.productLow(a3, b3, h);
- final double pn = q + h;
- s += r + ExtendedPrecision.twoSumLow(q, h, pn);
-
- // Final summation
- return getSum(pn, pn + s);
- }
-
- /**
- * @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) {
- // p/q are the standard scalar product old/new sum (alternating).
- // s is the sum of round-off parts.
- // pn is the final scalar product sum.
- double p = a1 * b1;
- double s = ExtendedPrecision.productLow(a1, b1, p);
- double h = a2 * b2;
- double r = ExtendedPrecision.productLow(a2, b2, h);
- final double q = p + h;
- s += ExtendedPrecision.twoSumLow(p, h, q) + r;
- h = a3 * b3;
- r = ExtendedPrecision.productLow(a3, b3, h);
- p = q + h;
- s += ExtendedPrecision.twoSumLow(q, h, p) + r;
- h = a4 * b4;
- r = ExtendedPrecision.productLow(a4, b4, h);
- final double pn = p + h;
- s += ExtendedPrecision.twoSumLow(p, h, pn) + r;
-
- // Final summation
- return getSum(pn, pn + s);
- }
-
- /**
- * Gets the final sum. This checks the high precision sum is finite, otherwise
- * returns the standard precision sum for the IEEE754 result.
- *
- * <p>The high precision sum may be non-finite due to input infinite
- * or NaN numbers or overflow in the summation. In all cases returning the
- * standard sum ensures the IEEE754 result.
- *
- * @param sum Standard sum.
- * @param hpSum High precision sum.
- * @return the sum
- */
- private static double getSum(double sum, double hpSum) {
- if (!Double.isFinite(hpSum)) {
- // Either we have split infinite numbers, some coefficients were NaNs,
- // or the sum overflowed.
- // Return the naive implementation for the IEEE754 result.
- return sum;
- }
- return hpSum;
- }
-}
diff --git a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Norm.java b/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Norm.java
index 65859b9..255215a 100644
--- a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Norm.java
+++ b/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Norm.java
@@ -219,9 +219,9 @@ public enum Norm {
private static double manhattan(final double x,
final double y,
final double z) {
- return Summation.value(Math.abs(x),
- Math.abs(y),
- Math.abs(z));
+ return Sum.of(Math.abs(x),
+ Math.abs(y),
+ Math.abs(z)).getAsDouble();
}
/** Computes the Manhattan norm.
@@ -234,17 +234,13 @@ public enum Norm {
* @see #of(double[])
*/
private static double manhattan(final double[] v) {
- double sum = 0d;
- double comp = 0d;
+ final Sum sum = Sum.create();
for (int i = 0; i < v.length; ++i) {
- final double x = Math.abs(v[i]);
- final double sx = sum + x;
- comp += ExtendedPrecision.twoSumLow(sum, x, sx);
- sum = sx;
+ sum.add(Math.abs(v[i]));
}
- return Summation.summationResult(sum, comp);
+ return sum.getAsDouble();
}
/** Computes the Euclidean norm.
diff --git a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Sum.java b/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Sum.java
new file mode 100644
index 0000000..9939e62
--- /dev/null
+++ b/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Sum.java
@@ -0,0 +1,258 @@
+/*
+ * 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 java.util.function.DoubleConsumer;
+import java.util.function.DoubleSupplier;
+
+/** Class providing accurate floating-point sums and linear combinations. The methods
+ * provided use compensated summation and multiplication techniques to reduce numerical errors.
+ * The approach is based on the <em>Sum2S</em> and <em>Dot2S</em> algorithms described in the
+ * 2005 paper <a href="https://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>.
+ *
+ * <p>Method results follow the standard rules for IEEE 754 addition. For example,
+ * if any input value is NaN, the result is NaN.
+ */
+public final class Sum implements DoubleSupplier, DoubleConsumer {
+
+ /** Standard sum. */
+ private double sum;
+
+ /** Compensation value. */
+ private double comp;
+
+ /** Construct a new instance with an initial value of zero.
+ */
+ private Sum() {
+ this(0d);
+ }
+
+ /** Construct a new instance with the given initial value.
+ * @param initialValue initial value
+ */
+ private Sum(final double initialValue) {
+ this.sum = initialValue;
+ }
+
+ /** Add a single term to this sum.
+ * @param t value to add
+ * @return this instance
+ */
+ public Sum add(final double t) {
+ final double newSum = sum + t;
+ comp += ExtendedPrecision.twoSumLow(sum, t, newSum);
+ sum = newSum;
+
+ return this;
+ }
+
+ /** Add an array of value to the sum.
+ * @param terms terms to add
+ * @return this instance
+ */
+ public Sum add(final double[] terms) {
+ for (double t : terms) {
+ add(t);
+ }
+
+ return this;
+ }
+
+ /** Add the high-accuracy product \(a b\) to this sum.
+ * @param a first factor
+ * @param b second factor
+ * @return this instance
+ */
+ public Sum addProduct(final double a, final double b) {
+ final double ab = a * b;
+ final double pLow = ExtendedPrecision.productLow(a, b, ab);
+
+ final double newSum = sum + ab;
+ comp += ExtendedPrecision.twoSumLow(sum, ab, newSum) + pLow;
+ sum = newSum;
+
+ return this;
+ }
+
+ /** Add \( \sum_i a_i b_i \). In other words, multiply each element
+ * in {@code a} with its corresponding element in {@code b} and add the product
+ * to the sum.
+ * @param a factors
+ * @param b factors
+ * @return this instance
+ * @throws IllegalArgumentException if the arrays do not have the same length
+ */
+ public Sum addProducts(final double[] a, final double[] b) {
+ final int len = a.length;
+ if (len != b.length) {
+ throw new IllegalArgumentException("Dimension mismatch: " + a.length + " != " + b.length);
+ }
+
+ for (int i = 0; i < len; ++i) {
+ addProduct(a[i], b[i]);
+ }
+
+ return this;
+ }
+
+ /** Add another sum to this sum.
+ * @param other sum to add
+ * @return this instance
+ */
+ public Sum add(final Sum other) {
+ // pull both values first to ensure there are
+ // no issues when adding a sum to itself
+ final double s = other.sum;
+ final double c = other.comp;
+
+ return add(s)
+ .add(c);
+ }
+
+ /** Add a single term to this sum. This is equivalent to {@link #add(double)}.
+ * @param value value to add
+ * @see #add(double)
+ */
+ @Override
+ public void accept(final double value) {
+ add(value);
+ }
+
+ /** Get the sum value.
+ * @return sum value as a double
+ */
+ @Override
+ public double getAsDouble() {
+ // compute and return the high precision sum if it is finite; otherwise,
+ // return the standard IEEE754 result
+ final double hpsum = sum + comp;
+ return Double.isFinite(hpsum) ?
+ hpsum :
+ sum;
+ }
+
+ /** Create a new sum instance with an initial value of zero.
+ * @return new sum instance
+ */
+ public static Sum create() {
+ return new Sum();
+ }
+
+ /** Return a new sum instance containing a single value.
+ * @param a value
+ * @return new sum instance
+ * @see #add(double)
+ */
+ public static Sum of(final double a) {
+ return new Sum(a);
+ }
+
+ /** Return a new sum instance containing the value \(a + b\).
+ * @param a first term
+ * @param b second term
+ * @return new sum instance
+ * @see #add(double)
+ */
+ public static Sum of(final double a, final double b) {
+ return new Sum(a).add(b);
+ }
+
+ /** Return a new sum instance containing the value \(a + b + c\).
+ * @param a first term
+ * @param b second term
+ * @param c third term
+ * @return new sum instance
+ * @see #add(double)
+ */
+ public static Sum of(final double a, final double b, final double c) {
+ return new Sum(a)
+ .add(b)
+ .add(c);
+ }
+
+ /** Return a new sum instance containing the value \(a + b + c + d\).
+ * @param a first term
+ * @param b second term
+ * @param c third term
+ * @param d fourth term
+ * @return new sum instance
+ * @see #add(double)
+ */
+ public static Sum of(final double a, final double b, final double c, final double d) {
+ return new Sum(a)
+ .add(b)
+ .add(c)
+ .add(d);
+ }
+
+ /** Return a new sum instance containing the sum of the given values.
+ * @param values input values
+ * @return new sum instance
+ * @see #add(double[])
+ */
+ public static Sum of(final double[] values) {
+ return new Sum().add(values);
+ }
+
+ /** Return a new sum instance containing the linear combination
+ * \(a_1 b_1 + a_2 b_2\).
+ * @param a1 first factor of first term
+ * @param b1 second factor of first term
+ * @param a2 first factor of second term
+ * @param b2 second factor of second term
+ * @return new sum instance
+ * @see #addProduct(double, double)
+ */
+ public static Sum ofProducts(final double a1, final double b1,
+ final double a2, final double b2) {
+ return new Sum()
+ .addProduct(a1, b1)
+ .addProduct(a2, b2);
+ }
+
+ /** Return a new sum instance containing the linear combination
+ * \(a_1 b_1 + a_2 b_2 + a_3 b_3\).
+ * @param a1 first factor of first term
+ * @param b1 second factor of first term
+ * @param a2 first factor of second term
+ * @param b2 second factor of second term
+ * @param a3 first factor of third term
+ * @param b3 second factor of third term
+ * @return new sum instance
+ * @see #addProduct(double, double)
+ */
+ public static Sum ofProducts(final double a1, final double b1,
+ final double a2, final double b2,
+ final double a3, final double b3) {
+ return new Sum()
+ .addProduct(a1, b1)
+ .addProduct(a2, b2)
+ .addProduct(a3, b3);
+ }
+
+ /** Return a new sum instance containing \( \sum_i a_i b_i \).
+ * @param a first set of factors
+ * @param b second set of factors
+ * @return new sum instance
+ * @see #addProducts(double[], double[])
+ */
+ public static Sum ofProducts(final double[] a, final double[] b) {
+ return new Sum().addProducts(a, b);
+ }
+}
diff --git a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Summation.java b/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Summation.java
deleted file mode 100644
index 3c26ab9..0000000
--- a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Summation.java
+++ /dev/null
@@ -1,179 +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;
-
-/** Class providing accurate floating-point summations. The methods provided
- * use a compensated summation technique to reduce numerical errors.
- * The approach is based on the <em>Sum2S</em> algorithm described in the
- * 2005 paper <a href="https://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>.
- *
- * <p>Method results follow the standard rules for IEEE 754 addition. For example,
- * if any input value is NaN, the result is NaN.
- */
-public final class Summation {
-
- /** Utility class; no instantiation. */
- private Summation() {}
-
- /** Compute the sum of the input values.
- * @param a first value
- * @param b second value
- * @param c third value
- * @return sum of the input values
- */
- public static double value(final double a, final double b, final double c) {
- double sum = a;
- double comp = 0d;
-
- final double sb = sum + b;
- comp += ExtendedPrecision.twoSumLow(sum, b, sb);
- sum = sb;
-
- final double sc = sum + c;
- comp += ExtendedPrecision.twoSumLow(sum, c, sc);
- sum = sc;
-
- return summationResult(sum, comp);
- }
-
- /** Compute the sum of the input values.
- * @param a first value
- * @param b second value
- * @param c third value
- * @param d fourth value
- * @return sum of the input values
- */
- public static double value(final double a, final double b, final double c, final double d) {
- double sum = a;
- double comp = 0d;
-
- final double sb = sum + b;
- comp += ExtendedPrecision.twoSumLow(sum, b, sb);
- sum = sb;
-
- final double sc = sum + c;
- comp += ExtendedPrecision.twoSumLow(sum, c, sc);
- sum = sc;
-
- final double sd = sum + d;
- comp += ExtendedPrecision.twoSumLow(sum, d, sd);
- sum = sd;
-
- return summationResult(sum, comp);
- }
-
- /** Compute the sum of the input values.
- * @param a first value
- * @param b second value
- * @param c third value
- * @param d fourth value
- * @param e fifth value
- * @return sum of the input values
- */
- public static double value(final double a, final double b, final double c, final double d,
- final double e) {
- double sum = a;
- double comp = 0d;
-
- final double sb = sum + b;
- comp += ExtendedPrecision.twoSumLow(sum, b, sb);
- sum = sb;
-
- final double sc = sum + c;
- comp += ExtendedPrecision.twoSumLow(sum, c, sc);
- sum = sc;
-
- final double sd = sum + d;
- comp += ExtendedPrecision.twoSumLow(sum, d, sd);
- sum = sd;
-
- final double se = sum + e;
- comp += ExtendedPrecision.twoSumLow(sum, e, se);
- sum = se;
-
- return summationResult(sum, comp);
- }
-
- /** Compute the sum of the input values.
- * @param a first value
- * @param b second value
- * @param c third value
- * @param d fourth value
- * @param e fifth value
- * @param f sixth value
- * @return sum of the input values
- */
- public static double value(final double a, final double b, final double c, final double d,
- final double e, final double f) {
- double sum = a;
- double comp = 0d;
-
- final double sb = sum + b;
- comp += ExtendedPrecision.twoSumLow(sum, b, sb);
- sum = sb;
-
- final double sc = sum + c;
- comp += ExtendedPrecision.twoSumLow(sum, c, sc);
- sum = sc;
-
- final double sd = sum + d;
- comp += ExtendedPrecision.twoSumLow(sum, d, sd);
- sum = sd;
-
- final double se = sum + e;
- comp += ExtendedPrecision.twoSumLow(sum, e, se);
- sum = se;
-
- final double sf = sum + f;
- comp += ExtendedPrecision.twoSumLow(sum, f, sf);
- sum = sf;
-
- return summationResult(sum, comp);
- }
-
- /** Compute the sum of the input values.
- * @param a array containing values to sum
- * @return sum of the input values
- */
- public static double value(final double[] a) {
- double sum = 0d;
- double comp = 0d;
-
- for (final double x : a) {
- final double s = sum + x;
- comp += ExtendedPrecision.twoSumLow(sum, x, s);
- sum = s;
- }
-
- return summationResult(sum, comp);
- }
-
- /** Return the final result from a summation operation.
- * @param sum standard sum value
- * @param comp compensation value
- * @return final summation result
- */
- static double summationResult(final double sum, final double comp) {
- // only add comp if finite; otherwise, return the raw sum
- // to comply with standard double addition rules
- return Double.isFinite(comp) ?
- sum + comp :
- sum;
- }
-}
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 99c9876..0000000
--- a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/LinearCombinationTest.java
+++ /dev/null
@@ -1,416 +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 java.math.BigDecimal;
-
-import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.rng.simple.RandomSource;
-import org.junit.jupiter.api.Assertions;
-import org.junit.jupiter.api.Test;
-
-/**
- * Test cases for the {@link LinearCombination} class.
- */
-class LinearCombinationTest {
- @Test
- void testDimensionMismatch() {
- Assertions.assertThrows(IllegalArgumentException.class,
- () -> LinearCombination.value(new double[1], new double[2]));
- }
-
- // MATH-1005
- @Test
- void testSingleElementArray() {
- final double[] a = {1.23456789};
- final double[] b = {98765432.1};
-
- Assertions.assertEquals(a[0] * b[0], LinearCombination.value(a, b));
- }
-
- @Test
- void testTwoSums() {
- final BigDecimal[] aFN = new BigDecimal[] {
- BigDecimal.valueOf(-1321008684645961L),
- BigDecimal.valueOf(-5774608829631843L),
- BigDecimal.valueOf(-7645843051051357L),
- };
- final BigDecimal[] aFD = new BigDecimal[] {
- BigDecimal.valueOf(268435456L),
- BigDecimal.valueOf(268435456L),
- BigDecimal.valueOf(8589934592L)
- };
- final BigDecimal[] bFN = new BigDecimal[] {
- BigDecimal.valueOf(-5712344449280879L),
- BigDecimal.valueOf(-4550117129121957L),
- BigDecimal.valueOf(8846951984510141L)
- };
- final BigDecimal[] bFD = new BigDecimal[] {
- BigDecimal.valueOf(2097152L),
- BigDecimal.valueOf(2097152L),
- BigDecimal.valueOf(131072L)
- };
-
- final int len = aFN.length;
- final double[] a = new double[len];
- final double[] b = new double[len];
- for (int i = 0; i < len; i++) {
- a[i] = aFN[i].doubleValue() / aFD[i].doubleValue();
- b[i] = bFN[i].doubleValue() / bFD[i].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);
- Assertions.assertEquals(abSumInline, abSumArray);
-
- // Compare with arbitrary precision computation.
- BigDecimal result = BigDecimal.ZERO;
- for (int i = 0; i < a.length; i++) {
- result = result.add(aFN[i].divide(aFD[i]).multiply(bFN[i].divide(bFD[i])));
- }
- final double expected = result.doubleValue();
- Assertions.assertEquals(expected, abSumInline, 1e-15);
-
- final double naive = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
- Assertions.assertTrue(Math.abs(naive - abSumInline) > 1.5);
- }
-
- @Test
- void testArrayVsInline() {
- final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_SHI_RO_256_PP);
-
- 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});
- Assertions.assertEquals(sInline, sArray);
-
- // Two sums.
- sInline = LinearCombination.value(u1, v1, u2, v2, u3, v3);
- sArray = LinearCombination.value(new double[] {u1, u2, u3},
- new double[] {v1, v2, v3});
- Assertions.assertEquals(sInline, sArray);
-
- // 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});
- Assertions.assertEquals(sInline, sArray);
- }
- }
-
- @Test
- 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);
-
- Assertions.assertEquals(abSumInline, abSumArray);
- Assertions.assertEquals(-1.8551294182586248737720779899, abSumInline, 1e-15);
-
- final double naive = scaledA[0] * scaledB[0] + scaledA[1] * scaledB[1] + scaledA[2] * scaledB[2];
- Assertions.assertTrue(Math.abs(naive - abSumInline) > 1.5);
- }
-
- @Test
- void testNonFinite() {
- 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},
- {1, Double.MAX_VALUE, 3, 4},
- {1, 2, Double.MAX_VALUE, 4},
- {1, Double.MAX_VALUE / 2, 3, -Double.MAX_VALUE / 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},
- {1, -2, 3, 4},
- {1, -2, 3, 4},
- {1, -2, 3, 4},
- };
-
- Assertions.assertEquals(-3,
- LinearCombination.value(a[0][0], b[0][0],
- a[0][1], b[0][1]));
- Assertions.assertEquals(6,
- LinearCombination.value(a[0][0], b[0][0],
- a[0][1], b[0][1],
- a[0][2], b[0][2]));
- Assertions.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]));
- Assertions.assertEquals(22, LinearCombination.value(a[0], b[0]));
-
- Assertions.assertEquals(Double.NEGATIVE_INFINITY,
- LinearCombination.value(a[1][0], b[1][0],
- a[1][1], b[1][1]));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY,
- LinearCombination.value(a[1][0], b[1][0],
- a[1][1], b[1][1],
- a[1][2], b[1][2]));
- Assertions.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]));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY, LinearCombination.value(a[1], b[1]));
-
- Assertions.assertEquals(-3,
- LinearCombination.value(a[2][0], b[2][0],
- a[2][1], b[2][1]));
- Assertions.assertEquals(Double.POSITIVE_INFINITY,
- LinearCombination.value(a[2][0], b[2][0],
- a[2][1], b[2][1],
- a[2][2], b[2][2]));
- Assertions.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]));
- Assertions.assertEquals(Double.POSITIVE_INFINITY, LinearCombination.value(a[2], b[2]));
-
- Assertions.assertEquals(Double.NEGATIVE_INFINITY,
- LinearCombination.value(a[3][0], b[3][0],
- a[3][1], b[3][1]));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY,
- LinearCombination.value(a[3][0], b[3][0],
- a[3][1], b[3][1],
- a[3][2], b[3][2]));
- Assertions.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]));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY, LinearCombination.value(a[3], b[3]));
-
- Assertions.assertEquals(Double.POSITIVE_INFINITY,
- LinearCombination.value(a[4][0], b[4][0],
- a[4][1], b[4][1]));
- Assertions.assertEquals(Double.POSITIVE_INFINITY,
- LinearCombination.value(a[4][0], b[4][0],
- a[4][1], b[4][1],
- a[4][2], b[4][2]));
- Assertions.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]));
- Assertions.assertEquals(Double.POSITIVE_INFINITY, LinearCombination.value(a[4], b[4]));
-
- Assertions.assertEquals(-3,
- LinearCombination.value(a[5][0], b[5][0],
- a[5][1], b[5][1]));
- Assertions.assertEquals(Double.POSITIVE_INFINITY,
- LinearCombination.value(a[5][0], b[5][0],
- a[5][1], b[5][1],
- a[5][2], b[5][2]));
- Assertions.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]));
- Assertions.assertEquals(Double.POSITIVE_INFINITY, LinearCombination.value(a[5], b[5]));
-
- Assertions.assertEquals(Double.POSITIVE_INFINITY,
- LinearCombination.value(a[6][0], b[6][0],
- a[6][1], b[6][1]));
- Assertions.assertEquals(Double.POSITIVE_INFINITY,
- LinearCombination.value(a[6][0], b[6][0],
- a[6][1], b[6][1],
- a[6][2], b[6][2]));
- Assertions.assertEquals(Double.NaN,
- 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]));
- Assertions.assertEquals(Double.NaN, LinearCombination.value(a[6], b[6]));
-
- Assertions.assertEquals(Double.NaN,
- LinearCombination.value(a[7][0], b[7][0],
- a[7][1], b[7][1]));
- Assertions.assertEquals(Double.NaN,
- LinearCombination.value(a[7][0], b[7][0],
- a[7][1], b[7][1],
- a[7][2], b[7][2]));
- Assertions.assertEquals(Double.NaN,
- 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]));
- Assertions.assertEquals(Double.NaN, LinearCombination.value(a[7], b[7]));
-
- Assertions.assertEquals(Double.NEGATIVE_INFINITY,
- LinearCombination.value(a[8][0], b[8][0],
- a[8][1], b[8][1]));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY,
- LinearCombination.value(a[8][0], b[8][0],
- a[8][1], b[8][1],
- a[8][2], b[8][2]));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY,
- LinearCombination.value(a[8][0], b[8][0],
- a[8][1], b[8][1],
- a[8][2], b[8][2],
- a[8][3], b[8][3]));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY, LinearCombination.value(a[8], b[8]));
-
- Assertions.assertEquals(-3,
- LinearCombination.value(a[9][0], b[9][0],
- a[9][1], b[9][1]));
- Assertions.assertEquals(Double.POSITIVE_INFINITY,
- LinearCombination.value(a[9][0], b[9][0],
- a[9][1], b[9][1],
- a[9][2], b[9][2]));
- Assertions.assertEquals(Double.POSITIVE_INFINITY,
- LinearCombination.value(a[9][0], b[9][0],
- a[9][1], b[9][1],
- a[9][2], b[9][2],
- a[9][3], b[9][3]));
- Assertions.assertEquals(Double.POSITIVE_INFINITY, LinearCombination.value(a[9], b[9]));
-
- Assertions.assertEquals(-Double.MAX_VALUE,
- LinearCombination.value(a[10][0], b[10][0],
- a[10][1], b[10][1]));
- Assertions.assertEquals(-Double.MAX_VALUE,
- LinearCombination.value(a[10][0], b[10][0],
- a[10][1], b[10][1],
- a[10][2], b[10][2]));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY,
- LinearCombination.value(a[10][0], b[10][0],
- a[10][1], b[10][1],
- a[10][2], b[10][2],
- a[10][3], b[10][3]));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY, LinearCombination.value(a[10], b[10]));
- }
-
- /**
- * This creates a scenario where the split product will overflow but the standard
- * precision computation will not. The result is expected to be in extended precision,
- * i.e. the method correctly detects and handles intermediate overflow.
- *
- * <p>Note: This test assumes that LinearCombination computes a split number
- * using Dekker's method. This can result in the high part of the number being
- * greater in magnitude than the the original number due to round-off in the split.
- */
- @Test
- void testOverflow() {
- // Create a simple dot product that is different in high precision and has
- // values that create a high part above the original number. This can be done using
- // a mantissa with almost all bits set to 1.
- final double x = Math.nextDown(2.0);
- final double y = -Math.nextDown(x);
- final double xxMxy = x * x + x * y;
- final double xxMxyHighPrecision = LinearCombination.value(x, x, x, y);
- Assertions.assertNotEquals(xxMxy, xxMxyHighPrecision, "High precision result should be different");
-
- // Scale it close to max value.
- // The current exponent is 0 so the combined scale must be 1023-1 as the
- // partial product x*x and x*y have an exponent 1 higher
- Assertions.assertEquals(0, Math.getExponent(x));
- Assertions.assertEquals(0, Math.getExponent(y));
-
- final double a1 = Math.scalb(x, 1022 - 30);
- final double b1 = Math.scalb(x, 30);
- final double a2 = a1;
- final double b2 = Math.scalb(y, 30);
- // Verify low precision result is scaled and finite
- final double sxxMxy = Math.scalb(xxMxy, 1022);
- Assertions.assertEquals(sxxMxy, a1 * b1 + a2 * b2);
- Assertions.assertTrue(Double.isFinite(sxxMxy));
-
- // High precision result using Dekker's multiplier.
- final double m = (1 << 27) + 1;
- // First demonstrate that Dekker's split will create overflow in the high part.
- double c;
- c = a1 * m;
- final double ha1 = c - (c - a1);
- c = b1 * m;
- final double hb1 = c - (c - b1);
- c = a2 * m;
- final double ha2 = c - (c - a2);
- c = b2 * m;
- final double hb2 = c - (c - b2);
- Assertions.assertTrue(Double.isFinite(ha1));
- Assertions.assertTrue(Double.isFinite(hb1));
- Assertions.assertTrue(Double.isFinite(ha2));
- Assertions.assertTrue(Double.isFinite(hb2));
- // High part should be bigger in magnitude
- Assertions.assertTrue(Math.abs(ha1) > Math.abs(a1));
- Assertions.assertTrue(Math.abs(hb1) > Math.abs(b1));
- Assertions.assertTrue(Math.abs(ha2) > Math.abs(a2));
- Assertions.assertTrue(Math.abs(hb2) > Math.abs(b2));
- Assertions.assertEquals(Double.POSITIVE_INFINITY, ha1 * hb1, "Expected split high part to overflow");
- Assertions.assertEquals(Double.NEGATIVE_INFINITY, ha2 * hb2, "Expected split high part to overflow");
-
- // LinearCombination should detect and handle intermediate overflow and return the
- // high precision result.
- final double expected = Math.scalb(xxMxyHighPrecision, 1022);
- Assertions.assertEquals(expected, LinearCombination.value(a1, b1, a2, b2));
- Assertions.assertEquals(expected, LinearCombination.value(a1, b1, a2, b2, 0, 0));
- Assertions.assertEquals(expected, LinearCombination.value(a1, b1, a2, b2, 0, 0, 0, 0));
- Assertions.assertEquals(expected, LinearCombination.value(new double[] {a1, a2}, new double[] {b1, b2}));
- }
-}
diff --git a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SumTest.java b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SumTest.java
new file mode 100644
index 0000000..70afcd5
--- /dev/null
+++ b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SumTest.java
@@ -0,0 +1,568 @@
+/*
+ * 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 java.math.BigDecimal;
+import java.math.MathContext;
+import java.util.Arrays;
+
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
+
+class SumTest {
+
+ @Test
+ void testSum_simple() {
+ // act/assert
+ Assertions.assertEquals(0d, Sum.create().getAsDouble());
+
+ assertSum(Math.PI, Math.PI);
+ assertSum(Math.PI + Math.E, Math.PI, Math.E);
+
+ assertSum(0, 0, 0, 0);
+ assertSum(6, 1, 2, 3);
+ assertSum(2, 1, -2, 3);
+
+ assertSum(Double.NaN, Double.NaN, 0, 0);
+ assertSum(Double.NaN, 0, Double.NaN, 0);
+ assertSum(Double.NaN, 0, 0, Double.NaN);
+
+ assertSum(Double.NaN, Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, 0);
+
+ assertSum(Double.POSITIVE_INFINITY, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE);
+
+ assertSum(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, 1, 1);
+ assertSum(Double.NEGATIVE_INFINITY, 1, Double.NEGATIVE_INFINITY, 1);
+ }
+
+ @Test
+ void testSumAccuracy() {
+ // arrange
+ final double a = 9.999999999;
+ final double b = Math.scalb(a, -53);
+ final double c = Math.scalb(a, -53);
+ final double d = Math.scalb(a, -27);
+ final double e = Math.scalb(a, -27);
+ final double f = Math.scalb(a, -50);
+
+ // act/assert
+ assertSumExact(a);
+
+ assertSumExact(a, b);
+ assertSumExact(b, a);
+
+ assertSumExact(a, b, c);
+ assertSumExact(c, b, a);
+
+ assertSumExact(a, b, c, d);
+ assertSumExact(d, c, b, a);
+
+ assertSumExact(a, -b, c, -d);
+ assertSumExact(d, -c, b, -a);
+
+ assertSumExact(a, b, c, d, e, f);
+ assertSumExact(f, e, d, c, b, a);
+
+ assertSumExact(a, -b, c, -d, e, f);
+ assertSumExact(f, -e, d, -c, b, -a);
+ }
+
+ @Test
+ void testAdd_sumInstance() {
+ // arrange
+ final double a = Math.PI;
+ final double b = Math.scalb(a, -53);
+ final double c = Math.scalb(a, -53);
+ final double d = Math.scalb(a, -27);
+ final double e = Math.scalb(a, -27);
+ final double f = Math.scalb(a, -50);
+
+ // act/assert
+ Assertions.assertEquals(exactSum(a, b, c, d), Sum.of(a, b, c, d).add(Sum.create()).getAsDouble());
+ Assertions.assertEquals(exactSum(a, a, b, c, d, e, f),
+ Sum.of(a, b)
+ .add(Sum.of(a, c))
+ .add(Sum.of(d, e, f)).getAsDouble());
+
+ final Sum s = Sum.of(a, b);
+ Assertions.assertEquals(exactSum(a, b, a, b), s.add(s).getAsDouble());
+ }
+
+ @Test
+ void testSumOfProducts_dimensionMismatch() {
+ // act/assert
+ final Sum sum = Sum.create();
+ Assertions.assertThrows(IllegalArgumentException.class,
+ () -> sum.addProducts(new double[1], new double[2]));
+
+ Assertions.assertThrows(IllegalArgumentException.class,
+ () -> Sum.ofProducts(new double[1], new double[2]));
+ }
+
+ @Test
+ void testSumOfProducts_singleElement() {
+ final double[] a = {1.23456789};
+ final double[] b = {98765432.1};
+
+ Assertions.assertEquals(a[0] * b[0], Sum.ofProducts(a, b).getAsDouble());
+ }
+
+ @Test
+ void testSumOfProducts() {
+ // arrange
+ final BigDecimal[] aFN = new BigDecimal[] {
+ BigDecimal.valueOf(-1321008684645961L),
+ BigDecimal.valueOf(-5774608829631843L),
+ BigDecimal.valueOf(-7645843051051357L),
+ };
+ final BigDecimal[] aFD = new BigDecimal[] {
+ BigDecimal.valueOf(268435456L),
+ BigDecimal.valueOf(268435456L),
+ BigDecimal.valueOf(8589934592L)
+ };
+ final BigDecimal[] bFN = new BigDecimal[] {
+ BigDecimal.valueOf(-5712344449280879L),
+ BigDecimal.valueOf(-4550117129121957L),
+ BigDecimal.valueOf(8846951984510141L)
+ };
+ final BigDecimal[] bFD = new BigDecimal[] {
+ BigDecimal.valueOf(2097152L),
+ BigDecimal.valueOf(2097152L),
+ BigDecimal.valueOf(131072L)
+ };
+
+ final int len = aFN.length;
+ final double[] a = new double[len];
+ final double[] b = new double[len];
+ for (int i = 0; i < len; i++) {
+ a[i] = aFN[i].doubleValue() / aFD[i].doubleValue();
+ b[i] = bFN[i].doubleValue() / bFD[i].doubleValue();
+ }
+
+ // act
+ final double sum = Sum.ofProducts(a, b).getAsDouble();
+
+ // assert
+ // Compare with arbitrary precision computation.
+ BigDecimal result = BigDecimal.ZERO;
+ for (int i = 0; i < a.length; i++) {
+ result = result.add(aFN[i].divide(aFD[i]).multiply(bFN[i].divide(bFD[i])));
+ }
+ final double expected = result.doubleValue();
+ Assertions.assertEquals(expected, sum, 1e-15);
+
+ final double naive = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+ Assertions.assertTrue(Math.abs(naive - sum) > 1.5);
+ }
+
+ @Test
+ void testSumOfProducts_huge() {
+ // arrange
+ 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);
+ }
+
+ // act
+ final double sum = Sum.ofProducts(scaledA[0], scaledB[0],
+ scaledA[1], scaledB[1],
+ scaledA[2], scaledB[2]).getAsDouble();
+
+ // assert
+ Assertions.assertEquals(-1.8551294182586248737720779899, sum, 1e-15);
+
+ final double naive = scaledA[0] * scaledB[0] + scaledA[1] * scaledB[1] + scaledA[2] * scaledB[2];
+ Assertions.assertTrue(Math.abs(naive - sum) > 1.5);
+ }
+
+ @Test
+ void testSumOfProducts_nonFinite() {
+ // arrange
+ 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},
+ {1, Double.MAX_VALUE, 3, 4},
+ {1, 2, Double.MAX_VALUE, 4},
+ {1, Double.MAX_VALUE / 2, 3, -Double.MAX_VALUE / 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},
+ {1, -2, 3, 4},
+ {1, -2, 3, 4},
+ {1, -2, 3, 4},
+ };
+
+ // act/assert
+ assertSumOfProducts(-3,
+ a[0][0], b[0][0],
+ a[0][1], b[0][1]);
+ assertSumOfProducts(6,
+ a[0][0], b[0][0],
+ a[0][1], b[0][1],
+ a[0][2], b[0][2]);
+ assertSumOfProducts(22, a[0], b[0]);
+
+ assertSumOfProducts(Double.NEGATIVE_INFINITY,
+ a[1][0], b[1][0],
+ a[1][1], b[1][1]);
+ assertSumOfProducts(Double.NEGATIVE_INFINITY,
+ a[1][0], b[1][0],
+ a[1][1], b[1][1],
+ a[1][2], b[1][2]);
+ assertSumOfProducts(Double.NEGATIVE_INFINITY, a[1], b[1]);
+
+ assertSumOfProducts(-3,
+ a[2][0], b[2][0],
+ a[2][1], b[2][1]);
+ assertSumOfProducts(Double.POSITIVE_INFINITY,
+ a[2][0], b[2][0],
+ a[2][1], b[2][1],
+ a[2][2], b[2][2]);
+ assertSumOfProducts(Double.POSITIVE_INFINITY, a[2], b[2]);
+
+ assertSumOfProducts(Double.NEGATIVE_INFINITY,
+ a[3][0], b[3][0],
+ a[3][1], b[3][1]);
+ assertSumOfProducts(Double.NEGATIVE_INFINITY,
+ a[3][0], b[3][0],
+ a[3][1], b[3][1],
+ a[3][2], b[3][2]);
+ assertSumOfProducts(Double.NEGATIVE_INFINITY, a[3], b[3]);
+
+ assertSumOfProducts(Double.POSITIVE_INFINITY,
+ a[4][0], b[4][0],
+ a[4][1], b[4][1]);
+ assertSumOfProducts(Double.POSITIVE_INFINITY,
+ a[4][0], b[4][0],
+ a[4][1], b[4][1],
+ a[4][2], b[4][2]);
+ assertSumOfProducts(Double.POSITIVE_INFINITY, a[4], b[4]);
+
+ assertSumOfProducts(-3,
+ a[5][0], b[5][0],
+ a[5][1], b[5][1]);
+ assertSumOfProducts(Double.POSITIVE_INFINITY,
+ a[5][0], b[5][0],
+ a[5][1], b[5][1],
+ a[5][2], b[5][2]);
+ assertSumOfProducts(Double.POSITIVE_INFINITY, a[5], b[5]);
+
+ assertSumOfProducts(Double.POSITIVE_INFINITY,
+ a[6][0], b[6][0],
+ a[6][1], b[6][1]);
+ assertSumOfProducts(Double.POSITIVE_INFINITY,
+ a[6][0], b[6][0],
+ a[6][1], b[6][1],
+ a[6][2], b[6][2]);
+ assertSumOfProducts(Double.NaN, a[6], b[6]);
+
+ assertSumOfProducts(Double.NaN,
+ a[7][0], b[7][0],
+ a[7][1], b[7][1]);
+ assertSumOfProducts(Double.NaN,
+ a[7][0], b[7][0],
+ a[7][1], b[7][1],
+ a[7][2], b[7][2]);
+ assertSumOfProducts(Double.NaN, a[7], b[7]);
+
+ assertSumOfProducts(Double.NEGATIVE_INFINITY,
+ a[8][0], b[8][0],
+ a[8][1], b[8][1]);
+ assertSumOfProducts(Double.NEGATIVE_INFINITY,
+ a[8][0], b[8][0],
+ a[8][1], b[8][1],
+ a[8][2], b[8][2]);
+ assertSumOfProducts(Double.NEGATIVE_INFINITY, a[8], b[8]);
+
+ assertSumOfProducts(-3,
+ a[9][0], b[9][0],
+ a[9][1], b[9][1]);
+ assertSumOfProducts(Double.POSITIVE_INFINITY,
+ a[9][0], b[9][0],
+ a[9][1], b[9][1],
+ a[9][2], b[9][2]);
+ assertSumOfProducts(Double.POSITIVE_INFINITY, a[9], b[9]);
+
+ assertSumOfProducts(-Double.MAX_VALUE,
+ a[10][0], b[10][0],
+ a[10][1], b[10][1]);
+ assertSumOfProducts(-Double.MAX_VALUE,
+ a[10][0], b[10][0],
+ a[10][1], b[10][1],
+ a[10][2], b[10][2]);
+ assertSumOfProducts(Double.NEGATIVE_INFINITY, a[10], b[10]);
+ }
+
+ /**
+ * This creates a scenario where the split product will overflow but the standard
+ * precision computation will not. The result is expected to be in extended precision,
+ * i.e. the method correctly detects and handles intermediate overflow.
+ *
+ * <p>Note: This test assumes that LinearCombination computes a split number
+ * using Dekker's method. This can result in the high part of the number being
+ * greater in magnitude than the the original number due to round-off in the split.
+ */
+ @Test
+ void testSumOfProducts_overflow() {
+ // Create a simple dot product that is different in high precision and has
+ // values that create a high part above the original number. This can be done using
+ // a mantissa with almost all bits set to 1.
+ final double x = Math.nextDown(2.0);
+ final double y = -Math.nextDown(x);
+ final double xxMxy = x * x + x * y;
+ final double xxMxyHighPrecision = Sum.ofProducts(x, x, x, y).getAsDouble();
+ Assertions.assertNotEquals(xxMxy, xxMxyHighPrecision, "High precision result should be different");
+
+ // Scale it close to max value.
+ // The current exponent is 0 so the combined scale must be 1023-1 as the
+ // partial product x*x and x*y have an exponent 1 higher
+ Assertions.assertEquals(0, Math.getExponent(x));
+ Assertions.assertEquals(0, Math.getExponent(y));
+
+ final double a1 = Math.scalb(x, 1022 - 30);
+ final double b1 = Math.scalb(x, 30);
+ final double a2 = a1;
+ final double b2 = Math.scalb(y, 30);
+ // Verify low precision result is scaled and finite
+ final double sxxMxy = Math.scalb(xxMxy, 1022);
+ Assertions.assertEquals(sxxMxy, a1 * b1 + a2 * b2);
+ Assertions.assertTrue(Double.isFinite(sxxMxy));
+
+ // High precision result using Dekker's multiplier.
+ final double m = (1 << 27) + 1;
+ // First demonstrate that Dekker's split will create overflow in the high part.
+ double c;
+ c = a1 * m;
+ final double ha1 = c - (c - a1);
+ c = b1 * m;
+ final double hb1 = c - (c - b1);
+ c = a2 * m;
+ final double ha2 = c - (c - a2);
+ c = b2 * m;
+ final double hb2 = c - (c - b2);
+ Assertions.assertTrue(Double.isFinite(ha1));
+ Assertions.assertTrue(Double.isFinite(hb1));
+ Assertions.assertTrue(Double.isFinite(ha2));
+ Assertions.assertTrue(Double.isFinite(hb2));
+ // High part should be bigger in magnitude
+ Assertions.assertTrue(Math.abs(ha1) > Math.abs(a1));
+ Assertions.assertTrue(Math.abs(hb1) > Math.abs(b1));
+ Assertions.assertTrue(Math.abs(ha2) > Math.abs(a2));
+ Assertions.assertTrue(Math.abs(hb2) > Math.abs(b2));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, ha1 * hb1, "Expected split high part to overflow");
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, ha2 * hb2, "Expected split high part to overflow");
+
+ // LinearCombination should detect and handle intermediate overflow and return the
+ // high precision result.
+ final double expected = Math.scalb(xxMxyHighPrecision, 1022);
+ assertSumOfProducts(expected, a1, b1, a2, b2);
+ assertSumOfProducts(expected, a1, b1, a2, b2, 0, 0);
+ assertSumOfProducts(expected, a1, b1, a2, b2, 0, 0, 0, 0);
+ }
+
+ @Test
+ void testMixedSingleTermAndProduct() {
+ // arrange
+ final double a = 9.999999999;
+ final double b = Math.scalb(a, -53);
+ final double c = Math.scalb(a, -53);
+ final double d = Math.scalb(a, -27);
+
+ // act/assert
+ Assertions.assertEquals(exactLinearCombination(1, a, -1, b, 2, c, 4, d),
+ Sum.create()
+ .add(a)
+ .add(-b)
+ .addProduct(2, c)
+ .addProduct(d, 4).getAsDouble());
+
+ Assertions.assertEquals(exactLinearCombination(1, a, -1, b, 2, c, 4, d),
+ Sum.create()
+ .addProduct(d, 4)
+ .add(a)
+ .addProduct(2, c)
+ .add(-b).getAsDouble());
+ }
+
+ @Test
+ public void testUnityValuesInProduct() {
+ // arrange
+ final double a = 9.999999999;
+ final double b = Math.scalb(a, -53);
+ final double c = Math.scalb(a, -53);
+ final double d = Math.scalb(a, -27);
+
+ // act/assert
+ Assertions.assertEquals(exactLinearCombination(1, a, -1, b, 2, c, 4, d),
+ Sum.create()
+ .addProduct(1, a)
+ .addProduct(-1, b)
+ .addProduct(2, c)
+ .addProduct(d, 4).getAsDouble());
+
+ Assertions.assertEquals(exactLinearCombination(1, a, -1, b, 2, c, 4, d),
+ Sum.create()
+ .addProduct(a, 1)
+ .addProduct(b, -1)
+ .addProduct(2, c)
+ .addProduct(d, 4).getAsDouble());
+ }
+
+ private static void assertSumExact(final double... values) {
+ final double exact = exactSum(values);
+ assertSum(exact, values);
+ }
+
+ private static void assertSum(final double expected, final double... values) {
+ // check non-array method variants
+ final int len = values.length;
+ if (len == 1) {
+ Assertions.assertEquals(expected, Sum.of(values[0]).getAsDouble());
+ } else if (len == 2) {
+ Assertions.assertEquals(expected, Sum.of(values[0], values[1]).getAsDouble());
+ } else if (len == 3) {
+ Assertions.assertEquals(expected, Sum.of(values[0], values[1], values[2]).getAsDouble());
+ } else if (len == 4) {
+ Assertions.assertEquals(expected, Sum.of(values[0], values[1], values[2], values[3]).getAsDouble());
+ }
+
+ // check use with add()
+ final Sum addAccumulator = Sum.create();
+ for (int i = 0; i < len; ++i) {
+ addAccumulator.add(values[i]);
+ }
+ Assertions.assertEquals(expected, addAccumulator.getAsDouble());
+
+ // check with accept()
+ final Sum acceptAccumulator = Sum.create();
+ for (int i = 0; i < len; ++i) {
+ acceptAccumulator.accept(values[i]);
+ }
+ Assertions.assertEquals(expected, acceptAccumulator.getAsDouble());
+
+ // check using stream
+ final Sum streamAccumulator = Sum.create();
+ Arrays.stream(values).forEach(streamAccumulator);
+ Assertions.assertEquals(expected, streamAccumulator.getAsDouble());
+
+ // check array instance method
+ Assertions.assertEquals(expected, Sum.create().add(values).getAsDouble());
+
+ // check array factory method
+ Assertions.assertEquals(expected, Sum.of(values).getAsDouble());
+ }
+
+ private static void assertSumOfProducts(final double expected, final double... args) {
+ final int halfLen = args.length / 2;
+
+ final double[] a = new double[halfLen];
+ final double[] b = new double[halfLen];
+ for (int i = 0; i < halfLen; ++i) {
+ a[i] = args[2 * i];
+ b[i] = args[(2 * i) + 1];
+ }
+
+ assertSumOfProducts(expected, a, b);
+ }
+
+ private static void assertSumOfProducts(final double expected, final double[] a, final double[] b) {
+ final int len = a.length;
+
+ // check non-array method variants
+ if (len == 2) {
+ Assertions.assertEquals(expected, Sum.ofProducts(a[0], b[0],
+ a[1], b[1]).getAsDouble());
+ } else if (len == 3) {
+ Assertions.assertEquals(expected, Sum.ofProducts(a[0], b[0],
+ a[1], b[1],
+ a[2], b[2]).getAsDouble());
+ }
+
+ // check use of addProduct()
+ final Sum accumulator = Sum.create();
+ for (int i = 0; i < len; ++i) {
+ accumulator.addProduct(a[i], b[i]);
+ }
+ Assertions.assertEquals(expected, accumulator.getAsDouble());
+
+ // check use of array instance method
+ Assertions.assertEquals(expected, Sum.create().addProducts(a, b).getAsDouble());
+
+ // check use of array factory method
+ Assertions.assertEquals(expected, Sum.ofProducts(a, b).getAsDouble());
+ }
+
+ /** Return the double estimation of the exact summation result computed with unlimited precision.
+ * @param values values to add
+ * @return double value closest to the exact result
+ */
+ private static double exactSum(final double... values) {
+ BigDecimal sum = BigDecimal.ZERO;
+ for (double value : values) {
+ sum = sum.add(new BigDecimal(value), MathContext.UNLIMITED);
+ }
+
+ return sum.doubleValue();
+ }
+
+ /** Return the double estimation of the exact linear combination result. Factors are
+ * listed sequentially in the argument array, e.g., {@code a1, b1, a2, b2, ...}.
+ * @param values linear combination input
+ * @return double value closest to the exact result
+ */
+ private static double exactLinearCombination(final double... values) {
+ final int halfLen = values.length / 2;
+
+ BigDecimal sum = BigDecimal.ZERO;
+ for (int i = 0; i < halfLen; ++i) {
+ final BigDecimal a = new BigDecimal(values[2 * i]);
+ final BigDecimal b = new BigDecimal(values[(2 * i) + 1]);
+
+ sum = sum.add(a.multiply(b, MathContext.UNLIMITED));
+ }
+
+ return sum.doubleValue();
+ }
+}
diff --git a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SummationTest.java b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SummationTest.java
deleted file mode 100644
index fd67160..0000000
--- a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SummationTest.java
+++ /dev/null
@@ -1,283 +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 java.math.BigDecimal;
-import java.math.MathContext;
-
-import org.junit.jupiter.api.Assertions;
-import org.junit.jupiter.api.Test;
-
-class SummationTest {
-
- @Test
- void test3n_simple() {
- // act/assert
- Assertions.assertEquals(0, Summation.value(0, 0, 0));
- Assertions.assertEquals(6, Summation.value(1, 2, 3));
- Assertions.assertEquals(2, Summation.value(1, -2, 3));
-
- Assertions.assertEquals(Double.NaN, Summation.value(Double.NaN, 0, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, Double.NaN, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, 0, Double.NaN));
-
- Assertions.assertEquals(Double.NaN, Summation.value(Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, 0));
-
- Assertions.assertEquals(Double.POSITIVE_INFINITY,
- Summation.value(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE));
-
- Assertions.assertEquals(Double.POSITIVE_INFINITY, Summation.value(Double.POSITIVE_INFINITY, 1, 1));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY, Summation.value(1, Double.NEGATIVE_INFINITY, 1));
- }
-
- @Test
- void test3n_accuracy() {
- // arrange
- final double a = 9.999999999;
- final double b = Math.scalb(a, -53);
- final double c = Math.scalb(a, -53);
-
- // act/assert
- assertValue(a, b, c);
- assertValue(c, b, a);
-
- assertValue(a, -b, c);
- assertValue(-c, b, -a);
- }
-
- @Test
- void test4n_simple() {
- // act/assert
- Assertions.assertEquals(0, Summation.value(0, 0, 0, 0));
- Assertions.assertEquals(10, Summation.value(1, 2, 3, 4));
- Assertions.assertEquals(-2, Summation.value(1, -2, 3, -4));
-
- Assertions.assertEquals(Double.NaN, Summation.value(Double.NaN, 0, 0, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, Double.NaN, 0, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, 0, Double.NaN, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, 0, 0, Double.NaN));
-
- Assertions.assertEquals(Double.NaN, Summation.value(Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, 0, 0));
-
- Assertions.assertEquals(Double.POSITIVE_INFINITY,
- Summation.value(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE));
-
- Assertions.assertEquals(Double.POSITIVE_INFINITY, Summation.value(Double.POSITIVE_INFINITY, 1, 1, 1));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY, Summation.value(1, Double.NEGATIVE_INFINITY, 1, 1));
- }
-
- @Test
- void test4n_accuracy() {
- // arrange
- final double a = 9.999999999;
- final double b = Math.scalb(a, -53);
- final double c = Math.scalb(a, -53);
- final double d = Math.scalb(a, -27);
-
- // act/assert
- assertValue(a, b, c, d);
- assertValue(d, c, b, a);
-
- assertValue(a, -b, c, -d);
- assertValue(d, -c, b, -a);
- }
-
- @Test
- void test5n_simple() {
- // act/assert
- Assertions.assertEquals(0, Summation.value(0, 0, 0, 0, 0));
- Assertions.assertEquals(15, Summation.value(1, 2, 3, 4, 5));
- Assertions.assertEquals(3, Summation.value(1, -2, 3, -4, 5));
-
- Assertions.assertEquals(Double.NaN, Summation.value(Double.NaN, 0, 0, 0, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, Double.NaN, 0, 0, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, 0, Double.NaN, 0, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, 0, 0, Double.NaN, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, 0, 0, 0, Double.NaN));
-
- Assertions.assertEquals(Double.NaN, Summation.value(
- Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, 0, 0, 0));
-
- Assertions.assertEquals(Double.POSITIVE_INFINITY, Summation.value(
- Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE));
-
- Assertions.assertEquals(Double.POSITIVE_INFINITY, Summation.value(Double.POSITIVE_INFINITY, 1, 1, 1, 1));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY, Summation.value(1, Double.NEGATIVE_INFINITY, 1, 1, 1));
- }
-
- @Test
- void test5n_accuracy() {
- // arrange
- final double a = 9.999999999;
- final double b = Math.scalb(a, -53);
- final double c = Math.scalb(a, -53);
- final double d = Math.scalb(a, -27);
- final double e = Math.scalb(a, -27);
-
- // act/assert
- assertValue(a, b, c, d, e);
- assertValue(e, d, c, b, a);
-
- assertValue(a, -b, c, -d, e);
- assertValue(-e, d, -c, b, -a);
- }
-
- @Test
- void test6n_simple() {
- // act/assert
- Assertions.assertEquals(0, Summation.value(0, 0, 0, 0, 0, 0));
- Assertions.assertEquals(21, Summation.value(1, 2, 3, 4, 5, 6));
- Assertions.assertEquals(-3, Summation.value(1, -2, 3, -4, 5, -6));
-
- Assertions.assertEquals(Double.NaN, Summation.value(Double.NaN, 0, 0, 0, 0, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, Double.NaN, 0, 0, 0, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, 0, Double.NaN, 0, 0, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, 0, 0, Double.NaN, 0, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, 0, 0, 0, Double.NaN, 0));
- Assertions.assertEquals(Double.NaN, Summation.value(0, 0, 0, 0, 0, Double.NaN));
-
- Assertions.assertEquals(Double.NaN, Summation.value(
- Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, 0, 0, 0, 0));
-
- Assertions.assertEquals(Double.POSITIVE_INFINITY, Summation.value(
- Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE,
- Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE));
-
- Assertions.assertEquals(Double.POSITIVE_INFINITY, Summation.value(Double.POSITIVE_INFINITY, 1, 1, 1, 1, 1));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY, Summation.value(1, Double.NEGATIVE_INFINITY, 1, 1, 1, 1));
- }
-
- @Test
- void test6n_accuracy() {
- // arrange
- final double a = 9.999999999;
- final double b = Math.scalb(a, -53);
- final double c = Math.scalb(a, -53);
- final double d = Math.scalb(a, -27);
- final double e = Math.scalb(a, -27);
- final double f = Math.scalb(a, -50);
-
- // act/assert
- assertValue(a, b, c, d, e, f);
- assertValue(f, e, d, c, b, a);
-
- assertValue(a, -b, c, -d, e, f);
- assertValue(f, -e, d, -c, b, -a);
- }
-
- @Test
- void testArray_simple() {
- // act/assert
- Assertions.assertEquals(0, Summation.value(new double[0]));
- Assertions.assertEquals(-1, Summation.value(new double[] {-1}));
- Assertions.assertEquals(0, Summation.value(new double[] {0, 0, 0}));
- Assertions.assertEquals(6, Summation.value(new double[] {1, 2, 3}));
- Assertions.assertEquals(2, Summation.value(new double[] {1, -2, 3}));
-
- Assertions.assertEquals(Double.MAX_VALUE, Summation.value(new double[] {Double.MAX_VALUE}));
- Assertions.assertEquals(Double.MIN_VALUE, Summation.value(new double[] {Double.MIN_VALUE}));
-
- Assertions.assertEquals(Double.NaN, Summation.value(new double[] {0d, Double.NaN}));
- Assertions.assertEquals(Double.NaN, Summation.value(new double[] {Double.POSITIVE_INFINITY, Double.NaN}));
- Assertions.assertEquals(Double.NaN,
- Summation.value(new double[] {Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY}));
-
- Assertions.assertEquals(Double.POSITIVE_INFINITY,
- Summation.value(new double[] {Double.MAX_VALUE, Double.MAX_VALUE}));
- Assertions.assertEquals(Double.POSITIVE_INFINITY,
- Summation.value(new double[] {Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY}));
- Assertions.assertEquals(Double.NEGATIVE_INFINITY,
- Summation.value(new double[] {Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY}));
- }
-
- @Test
- void testArray_accuracy() {
- // arrange
- final double a = 9.999999999;
- final double b = Math.scalb(a, -53);
- final double c = Math.scalb(a, -53);
- final double d = Math.scalb(a, -27);
- final double e = Math.scalb(a, -27);
- final double f = Math.scalb(a, -50);
-
- // act/assert
- assertArrayValue(a);
-
- assertArrayValue(a, b);
- assertArrayValue(b, a);
-
- assertArrayValue(a, b, c);
- assertArrayValue(c, b, a);
-
- assertArrayValue(a, b, c, d);
- assertArrayValue(d, c, b, a);
-
- assertArrayValue(a, -b, c, -d);
- assertArrayValue(d, -c, b, -a);
-
- assertArrayValue(a, b, c, d, e, f);
- assertArrayValue(f, e, d, c, b, a);
-
- assertArrayValue(a, -b, c, -d, e, f);
- assertArrayValue(f, -e, d, -c, b, -a);
- }
-
- private static void assertValue(final double a, final double b, final double c) {
- final double computed = Summation.value(a, b, c);
- assertComputedValue(computed, a, b, c);
- }
-
- private static void assertValue(final double a, final double b, final double c, final double d) {
- final double computed = Summation.value(a, b, c, d);
- assertComputedValue(computed, a, b, c, d);
- }
-
- private static void assertValue(final double a, final double b, final double c, final double d,
- final double e) {
- final double computed = Summation.value(a, b, c, d, e);
- assertComputedValue(computed, a, b, c, d, e);
- }
-
- private static void assertValue(final double a, final double b, final double c, final double d,
- final double e, final double f) {
- final double computed = Summation.value(a, b, c, d, e, f);
- assertComputedValue(computed, a, b, c, d, e, f);
- }
-
- private static void assertArrayValue(final double... values) {
- final double computed = Summation.value(values);
- assertComputedValue(computed, values);
- }
-
- private static void assertComputedValue(final double computed, final double... values) {
- final double exact = computeExact(values);
- Assertions.assertEquals(exact, computed);
- }
-
- /** Return the double estimation of the exact summation result computed with unlimited precision.
- * @param values values to add
- * @return double value closest to the exact result
- */
- private static double computeExact(final double... values) {
- BigDecimal sum = BigDecimal.ZERO;
- for (double value : values) {
- sum = sum.add(new BigDecimal(value), MathContext.UNLIMITED);
- }
-
- return sum.doubleValue();
- }
-}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/EuclideanNormAlgorithms.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/EuclideanNormAlgorithms.java
index ec8f689..5f65f84 100644
--- a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/EuclideanNormAlgorithms.java
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/EuclideanNormAlgorithms.java
@@ -20,6 +20,8 @@ import java.math.BigDecimal;
import java.math.MathContext;
import java.util.function.ToDoubleFunction;
+import org.apache.commons.numbers.core.Sum;
+
/** Class containing various Euclidean norm computation methods for comparison.
*/
public final class EuclideanNormAlgorithms {
@@ -370,7 +372,7 @@ public final class EuclideanNormAlgorithms {
}
rescale = 0x1.0p-600;
}
- return Math.sqrt(org.apache.commons.numbers.core.LinearCombination.value(x, x)) * rescale;
+ return Math.sqrt(Sum.ofProducts(x, x).getAsDouble()) * rescale;
}
}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/LinearCombinationPerformance.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/LinearCombinationPerformance.java
index 0590842..9873e22 100644
--- a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/LinearCombinationPerformance.java
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/LinearCombinationPerformance.java
@@ -22,7 +22,7 @@ import java.util.Arrays;
import java.util.concurrent.TimeUnit;
import java.util.function.IntFunction;
-import org.apache.commons.numbers.core.LinearCombination;
+import org.apache.commons.numbers.core.Sum;
import org.apache.commons.numbers.examples.jmh.core.LinearCombination.FourD;
import org.apache.commons.numbers.examples.jmh.core.LinearCombination.ND;
import org.apache.commons.numbers.examples.jmh.core.LinearCombination.ThreeD;
@@ -230,10 +230,15 @@ public class LinearCombinationPerformance {
@Setup
public void setup() {
if ("current".endsWith(name)) {
- twod = LinearCombination::value;
- threed = LinearCombination::value;
- fourd = LinearCombination::value;
- nd = LinearCombination::value;
+ twod = (a1, b1, a2, b2) -> Sum.ofProducts(a1, b1, a2, b2).getAsDouble();
+ threed = (a1, b1, a2, b2, a3, b3) -> Sum.ofProducts(a1, b1, a2, b2, a3, b3).getAsDouble();
+ fourd = (a1, b1, a2, b2, a3, b3, a4, b4) ->
+ Sum.create()
+ .addProduct(a1, b1)
+ .addProduct(a2, b2)
+ .addProduct(a3, b3)
+ .addProduct(a4, b4).getAsDouble();
+ nd = (a, b) -> Sum.ofProducts(a, b).getAsDouble();
return;
}
// All implementations below are expected to implement all the interfaces.
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/NormPerformance.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/NormPerformance.java
index 1792bb3..c6475ca 100644
--- a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/NormPerformance.java
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/NormPerformance.java
@@ -184,7 +184,7 @@ public class NormPerformance {
eval(v -> Math.hypot(v[0], v[1]), input, bh);
}
- /** Compute the performance of the {@link Norm#of(double, double)} method.
+ /** Compute the performance of the {@link Norm#L2} 2D method.
* @param input benchmark input
* @param bh blackhole
*/
@@ -193,7 +193,7 @@ public class NormPerformance {
eval(v -> Norm.L2.of(v[0], v[1]), input, bh);
}
- /** Compute the performance of the {@link Norm#of(double, double, double)} method.
+ /** Compute the performance of the {@link Norm#L2} 3D norm computation.
* @param input benchmark input
* @param bh blackhole
*/
@@ -202,7 +202,7 @@ public class NormPerformance {
eval(v -> Norm.L2.of(v[0], v[1], v[2]), input, bh);
}
- /** Compute the performance of the {@link Norm#of(double[])} method.
+ /** Compute the performance of the {@link Norm#L2} array norm method.
* @param input benchmark input
* @param bh blackhole
*/
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/SumPerformance.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/SumPerformance.java
new file mode 100644
index 0000000..0dd6732
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/SumPerformance.java
@@ -0,0 +1,183 @@
+/*
+ * 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.examples.jmh.core;
+
+import java.util.concurrent.TimeUnit;
+import java.util.function.ToDoubleBiFunction;
+import java.util.function.ToDoubleFunction;
+
+import org.apache.commons.numbers.core.Sum;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.simple.RandomSource;
+import org.openjdk.jmh.annotations.Benchmark;
+import org.openjdk.jmh.annotations.BenchmarkMode;
+import org.openjdk.jmh.annotations.Fork;
+import org.openjdk.jmh.annotations.Measurement;
+import org.openjdk.jmh.annotations.Mode;
+import org.openjdk.jmh.annotations.OutputTimeUnit;
+import org.openjdk.jmh.annotations.Param;
+import org.openjdk.jmh.annotations.Scope;
+import org.openjdk.jmh.annotations.Setup;
+import org.openjdk.jmh.annotations.State;
+import org.openjdk.jmh.annotations.Warmup;
+import org.openjdk.jmh.infra.Blackhole;
+
+/**
+ * Executes a benchmark to measure the speed of operations in the {@link Sum} class.
+ */
+@BenchmarkMode(Mode.AverageTime)
+@OutputTimeUnit(TimeUnit.NANOSECONDS)
+@Warmup(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
+@Measurement(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
+@State(Scope.Benchmark)
+@Fork(value = 1, jvmArgs = {"-server", "-Xms512M", "-Xmx512M"})
+public class SumPerformance {
+ /**
+ * The seed to use to create the random benchmark input.
+ * Using a fixed seed ensures the same values are created across benchmarks.
+ */
+ private static final long SEED = System.currentTimeMillis();
+
+ /** Class providing double arrays for benchmarks.
+ */
+ @State(Scope.Benchmark)
+ public static class ArrayInput {
+
+ /** Number of array samples. */
+ @Param("100000")
+ private int samples;
+
+ /** Number of values in each input array. */
+ @Param("50")
+ private int len;
+
+ /** Minimum possible double exponent. */
+ @Param("-550")
+ private int minExp;
+
+ /** Maximum possible double exponent. */
+ @Param("+550")
+ private int maxExp;
+
+ /** Range of exponents within a single array. */
+ @Param("26")
+ private int expRange;
+
+ /** First set of input arrays. */
+ private double[][] a;
+
+ /** Second set of input arrays. */
+ private double[][] b;
+
+ /** Get the first set of input arrays.
+ * @return first set of input arrays
+ */
+ public double[][] getA() {
+ return a;
+ }
+
+ /** Get the second set of input arrays.
+ * @return second set of input arrays
+ */
+ public double[][] getB() {
+ return b;
+ }
+
+ /** Create the input arrays for the instance. */
+ @Setup
+ public void createArrays() {
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP, SEED);
+
+ a = new double[samples][];
+ b = new double[samples][];
+ for (int i = 0; i < samples; ++i) {
+ // pick a general range for the array element exponents and then
+ // create values within that range
+ final int vMidExp = rng.nextInt(maxExp - minExp + 1) + minExp;
+ final int vExpRadius = expRange / 2;
+ final int vMinExp = vMidExp - vExpRadius;
+ final int vMaxExp = vMidExp + vExpRadius;
+
+ a[i] = DoubleUtils.randomArray(len, vMinExp, vMaxExp, rng);
+ b[i] = DoubleUtils.randomArray(len, vMinExp, vMaxExp, rng);
+ }
+ }
+ }
+
+ /** Run a benchmark for a function that accepts a single array and produces a double result.
+ * @param input benchmark input
+ * @param bh data sink
+ * @param fn function to benchmark
+ */
+ private static void runSingle(final ArrayInput input, final Blackhole bh,
+ final ToDoubleFunction<double[]> fn) {
+ final double[][] a = input.getA();
+ for (int i = 0; i < a.length; ++i) {
+ bh.consume(fn.applyAsDouble(a[i]));
+ }
+ }
+
+ /** Run a benchmark for a function that accepts a two arrays and produces a double result.
+ * @param input benchmark input
+ * @param bh data sink
+ * @param fn function to benchmark
+ */
+ private static void runDouble(final ArrayInput input, final Blackhole bh,
+ final ToDoubleBiFunction<double[], double[]> fn) {
+ final double[][] a = input.getA();
+ final double[][] b = input.getB();
+ for (int i = 0; i < a.length; ++i) {
+ bh.consume(fn.applyAsDouble(a[i], b[i]));
+ }
+ }
+
+ /** Benchmark baseline for functions that use a single input array.
+ * @param input benchmark input
+ * @param bh data sink
+ */
+ @Benchmark
+ public void baselineSingle(final ArrayInput input, final Blackhole bh) {
+ runSingle(input, bh, a -> 0d);
+ }
+
+ /** Benchmark baseline for functions that use two input arrays.
+ * @param input benchmark input
+ * @param bh data sink
+ */
+ @Benchmark
+ public void baselineDouble(final ArrayInput input, final Blackhole bh) {
+ runDouble(input, bh, (a, b) -> 0d);
+ }
+
+ /** Benchmark testing {@link Sum} addition performance.
+ * @param input benchmark input
+ * @param bh data sink
+ */
+ @Benchmark
+ public void sum(final ArrayInput input, final Blackhole bh) {
+ runSingle(input, bh, a -> Sum.of(a).getAsDouble());
+ }
+
+ /** Benchmark testing {@link Sum} linear combination performance.
+ * @param input benchmark input
+ * @param bh data sink
+ */
+ @Benchmark
+ public void sumOfProducts(final ArrayInput input, final Blackhole bh) {
+ runDouble(input, bh, (a, b) -> Sum.ofProducts(a, b).getAsDouble());
+ }
+}