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:24 UTC
[commons-numbers] branch master updated (d4c00eb -> 93c45e9)
This is an automated email from the ASF dual-hosted git repository.
erans pushed a change to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git.
from d4c00eb Javadoc.
new 06a613d NUMBERS-163: replacing LinearCombination and Summation with Sum class
new 93c45e9 NUMBERS-163: Nits.
The 2 revisions listed above as "new" are entirely new to this
repository and will be described in separate emails. The revisions
listed as "add" were already present in the repository and have only
been added to this reference.
Summary of changes:
.../org/apache/commons/numbers/angle/CosAngle.java | 5 +-
.../commons/numbers/core/LinearCombination.java | 222 --------
.../java/org/apache/commons/numbers/core/Norm.java | 17 +-
.../java/org/apache/commons/numbers/core/Sum.java | 205 ++++++++
.../org/apache/commons/numbers/core/Summation.java | 179 -------
.../numbers/core/LinearCombinationTest.java | 416 ---------------
.../org/apache/commons/numbers/core/SumTest.java | 556 +++++++++++++++++++++
.../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, 969 insertions(+), 1122 deletions(-)
delete mode 100644 commons-numbers-core/src/main/java/org/apache/commons/numbers/core/LinearCombination.java
create mode 100644 commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Sum.java
delete mode 100644 commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Summation.java
delete mode 100644 commons-numbers-core/src/test/java/org/apache/commons/numbers/core/LinearCombinationTest.java
create mode 100644 commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SumTest.java
delete mode 100644 commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SummationTest.java
create mode 100644 commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/core/SumPerformance.java
[commons-numbers] 02/02: NUMBERS-163: Nits.
Posted by er...@apache.org.
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 93c45e9dcde7f283c620508daacaf36383925596
Author: Gilles Sadowski <gi...@gmail.com>
AuthorDate: Sun Jun 20 04:10:59 2021 +0200
NUMBERS-163: Nits.
Javadoc style to match style in (most) other files.
Removed factory methods that were just one-liner replacements.
Added "varargs" (any number of terms) methods.
Explicit initialization of "comp" instance variable.
Closes #97.
---
.../java/org/apache/commons/numbers/core/Norm.java | 7 +-
.../java/org/apache/commons/numbers/core/Sum.java | 239 ++++++++-------------
.../org/apache/commons/numbers/core/SumTest.java | 16 +-
3 files changed, 99 insertions(+), 163 deletions(-)
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 255215a..92587fc 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,10 @@ public enum Norm {
private static double manhattan(final double x,
final double y,
final double z) {
- return Sum.of(Math.abs(x),
- Math.abs(y),
- Math.abs(z)).getAsDouble();
+ return Sum.of(Math.abs(x))
+ .add(Math.abs(y))
+ .add(Math.abs(z))
+ .getAsDouble();
}
/** Computes the Manhattan 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
index 9939e62..390649a 100644
--- 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
@@ -19,40 +19,42 @@ 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">
+/**
+ * Class providing accurate floating-point sums and linear combinations.
+ *
+ * In order to reduce errors, compensated summations and multiplications
+ * are performed according to the <em>Sum2S</em> and <em>Dot2S</em>
+ * algorithms described in
+ * <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>.
+ * and Shin'ichi Oishi (<em>SIAM J. Sci. Comput</em>, 2005).
*
- * <p>Method results follow the standard rules for IEEE 754 addition. For example,
- * if any input value is NaN, the result is NaN.
+ * <p>Results follow the IEEE 754 rules for addition: For example, if any
+ * input value is {@link Double#NaN}, the result is {@link Double#NaN}.
*/
-public final class Sum implements DoubleSupplier, DoubleConsumer {
-
+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
+ /**
+ * Constructs a new instance with the given initial value.
+ *
+ * @param initialValue Initial value.
*/
private Sum(final double initialValue) {
- this.sum = initialValue;
+ sum = initialValue;
+ comp = 0d;
}
- /** Add a single term to this sum.
- * @param t value to add
- * @return this instance
+ /**
+ * Adds 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;
@@ -62,11 +64,13 @@ public final class Sum implements DoubleSupplier, DoubleConsumer {
return this;
}
- /** Add an array of value to the sum.
- * @param terms terms to add
- * @return this instance
+ /**
+ * Adds values from the given array to the sum.
+ *
+ * @param terms Terms to add.
+ * @return this instance.
*/
- public Sum add(final double[] terms) {
+ public Sum add(final double... terms) {
for (double t : terms) {
add(t);
}
@@ -74,12 +78,15 @@ public final class Sum implements DoubleSupplier, DoubleConsumer {
return this;
}
- /** Add the high-accuracy product \(a b\) to this sum.
- * @param a first factor
- * @param b second factor
+ /**
+ * Adds the high-accuracy product \( a b \) to this sum.
+ *
+ * @param a Factor
+ * @param b Factor.
* @return this instance
*/
- public Sum addProduct(final double a, final double b) {
+ public Sum addProduct(final double a,
+ final double b) {
final double ab = a * b;
final double pLow = ExtendedPrecision.productLow(a, b, ab);
@@ -90,18 +97,20 @@ public final class Sum implements DoubleSupplier, DoubleConsumer {
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
+ /**
+ * Adds \( \sum_i a_i b_i \) to this 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) {
+ 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);
+ throw new IllegalArgumentException("Dimension mismatch: " +
+ a.length + " != " + b.length);
}
for (int i = 0; i < len; ++i) {
@@ -111,22 +120,27 @@ public final class Sum implements DoubleSupplier, DoubleConsumer {
return this;
}
- /** Add another sum to this sum.
- * @param other sum to add
- * @return this instance
+ /**
+ * Adds 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
+ // 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);
+ return add(s).add(c);
}
- /** Add a single term to this sum. This is equivalent to {@link #add(double)}.
- * @param value value to add
+ /**
+ * Adds a single term to this sum.
+ * This is equivalent to {@link #add(double)}.
+ *
+ * @param value Value to add.
+ *
* @see #add(double)
*/
@Override
@@ -134,125 +148,58 @@ public final class Sum implements DoubleSupplier, DoubleConsumer {
add(value);
}
- /** Get the sum value.
- * @return sum value as a double
+ /**
+ * Gets the sum value.
+ *
+ * @return the sum value.
*/
@Override
public double getAsDouble() {
- // compute and return the high precision sum if it is finite; otherwise,
- // return the standard IEEE754 result
+ // High-precision value if it is finite, standard IEEE754 result otherwise.
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
+ /**
+ * Creates a new instance with an initial value of zero.
+ *
+ * @return a new instance.
*/
public static Sum create() {
- return new Sum();
+ return new Sum(0d);
}
- /** Return a new sum instance containing a single value.
- * @param a value
- * @return new sum instance
- * @see #add(double)
+ /**
+ * Creates an instance initialized to the given value.
+ *
+ * @param a Initial value.
+ * @return a new instance.
*/
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)
+ /**
+ * Creates an instance containing the sum of the given values.
+ *
+ * @param values Values to add.
+ * @return a new instance.
*/
- 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);
+ public static Sum of(final double... values) {
+ return create().add(values);
}
- /** 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[])
+ /**
+ * Creates a new instance containing \( \sum_i a_i b_i \).
+ *
+ * @param a Factors.
+ * @param b Factors.
+ * @return a new instance.
*/
- public static Sum ofProducts(final double[] a, final double[] b) {
- return new Sum().addProducts(a, b);
+ public static Sum ofProducts(final double[] a,
+ final double[] b) {
+ return create().addProducts(a, b);
}
}
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
index 70afcd5..5d38260 100644
--- 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
@@ -193,9 +193,7 @@ class SumTest {
}
// act
- final double sum = Sum.ofProducts(scaledA[0], scaledB[0],
- scaledA[1], scaledB[1],
- scaledA[2], scaledB[2]).getAsDouble();
+ final double sum = Sum.ofProducts(scaledA, scaledB).getAsDouble();
// assert
Assertions.assertEquals(-1.8551294182586248737720779899, sum, 1e-15);
@@ -352,7 +350,7 @@ class SumTest {
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();
+ final double xxMxyHighPrecision = Sum.create().addProduct(x, x).addProduct(x, y).getAsDouble();
Assertions.assertNotEquals(xxMxy, xxMxyHighPrecision, "High precision result should be different");
// Scale it close to max value.
@@ -510,16 +508,6 @@ class SumTest {
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) {
[commons-numbers] 01/02: NUMBERS-163: replacing LinearCombination
and Summation with Sum class
Posted by er...@apache.org.
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());
+ }
+}