You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ma...@apache.org on 2021/06/08 01:02:16 UTC
[commons-numbers] branch master updated (64262fd -> 63ad23b)
This is an automated email from the ASF dual-hosted git repository.
mattjuntunen pushed a change to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git.
from 64262fd Track changes.
new 67c9c58 adding legitimate floating point equality check to spotbugs exclude
new 63ad23b NUMBERS-156: replacing SafeNorm with Norms and associated Summation class
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:
.travis.yml | 2 +-
commons-numbers-arrays/pom.xml | 4 +
.../apache/commons/numbers/arrays/CosAngle.java | 2 +-
.../commons/numbers/arrays/ExtendedPrecision.java | 21 +
.../org/apache/commons/numbers/arrays/Norms.java | 448 ++++++++++++
.../apache/commons/numbers/arrays/SafeNorm.java | 91 ---
.../apache/commons/numbers/arrays/Summation.java | 179 +++++
.../commons/numbers/arrays/DoubleTestUtils.java | 69 ++
.../numbers/arrays/ExtendedPrecisionTest.java | 17 +
.../apache/commons/numbers/arrays/NormsTest.java | 597 ++++++++++++++++
.../commons/numbers/arrays/SafeNormTest.java | 107 ---
.../commons/numbers/arrays/SummationTest.java | 283 ++++++++
commons-numbers-examples/examples-jmh/pom.xml | 39 ++
.../numbers/examples/jmh/arrays/DoubleUtils.java | 72 ++
.../arrays/EuclideanNormAlgorithmPerformance.java | 225 +++++++
.../jmh/arrays/EuclideanNormAlgorithms.java | 749 +++++++++++++++++++++
.../jmh/arrays/EuclideanNormEvaluator.java | 249 +++++++
.../examples/jmh/arrays/NormsPerformance.java | 198 ++++++
.../jmh/arrays/EuclideanNormAccuracyTest.java | 141 ++++
pom.xml | 9 +
.../resources/spotbugs/spotbugs-exclude-filter.xml | 5 +
21 files changed, 3307 insertions(+), 200 deletions(-)
create mode 100644 commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Norms.java
delete mode 100644 commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/SafeNorm.java
create mode 100644 commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Summation.java
create mode 100644 commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/DoubleTestUtils.java
create mode 100644 commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/NormsTest.java
delete mode 100644 commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/SafeNormTest.java
create mode 100644 commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/SummationTest.java
create mode 100644 commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/DoubleUtils.java
create mode 100644 commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAlgorithmPerformance.java
create mode 100644 commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAlgorithms.java
create mode 100644 commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormEvaluator.java
create mode 100644 commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/NormsPerformance.java
create mode 100644 commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAccuracyTest.java
[commons-numbers] 01/02: adding legitimate floating point equality
check to spotbugs exclude
Posted by ma...@apache.org.
This is an automated email from the ASF dual-hosted git repository.
mattjuntunen pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit 67c9c5807231ace9b6a26c293ab462e8dd88624e
Author: Matt Juntunen <ma...@apache.org>
AuthorDate: Mon Jun 7 06:36:37 2021 -0400
adding legitimate floating point equality check to spotbugs exclude
---
src/main/resources/spotbugs/spotbugs-exclude-filter.xml | 5 +++++
1 file changed, 5 insertions(+)
diff --git a/src/main/resources/spotbugs/spotbugs-exclude-filter.xml b/src/main/resources/spotbugs/spotbugs-exclude-filter.xml
index c3012e8..071037f 100644
--- a/src/main/resources/spotbugs/spotbugs-exclude-filter.xml
+++ b/src/main/resources/spotbugs/spotbugs-exclude-filter.xml
@@ -41,6 +41,11 @@
<Method name="sqrt"/>
<BugPattern name="FE_FLOATING_POINT_EQUALITY"/>
</Match>
+ <Match>
+ <Class name="org.apache.commons.numbers.angle.Angle$Normalizer"/>
+ <Method name="applyAsDouble"/>
+ <BugPattern name="FE_FLOATING_POINT_EQUALITY"/>
+ </Match>
<Match>
<!-- Benchmark state classes can expose internal representations -->
[commons-numbers] 02/02: NUMBERS-156: replacing SafeNorm with Norms
and associated Summation class
Posted by ma...@apache.org.
This is an automated email from the ASF dual-hosted git repository.
mattjuntunen pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit 63ad23b3bdcc8a98b8fdfd24093ec8ac8e3923b8
Author: Matt Juntunen <ma...@apache.org>
AuthorDate: Mon Jun 7 06:38:27 2021 -0400
NUMBERS-156: replacing SafeNorm with Norms and associated Summation class
---
.travis.yml | 2 +-
commons-numbers-arrays/pom.xml | 4 +
.../apache/commons/numbers/arrays/CosAngle.java | 2 +-
.../commons/numbers/arrays/ExtendedPrecision.java | 21 +
.../org/apache/commons/numbers/arrays/Norms.java | 448 ++++++++++++
.../apache/commons/numbers/arrays/SafeNorm.java | 91 ---
.../apache/commons/numbers/arrays/Summation.java | 179 +++++
.../commons/numbers/arrays/DoubleTestUtils.java | 69 ++
.../numbers/arrays/ExtendedPrecisionTest.java | 17 +
.../apache/commons/numbers/arrays/NormsTest.java | 597 ++++++++++++++++
.../commons/numbers/arrays/SafeNormTest.java | 107 ---
.../commons/numbers/arrays/SummationTest.java | 283 ++++++++
commons-numbers-examples/examples-jmh/pom.xml | 39 ++
.../numbers/examples/jmh/arrays/DoubleUtils.java | 72 ++
.../arrays/EuclideanNormAlgorithmPerformance.java | 225 +++++++
.../jmh/arrays/EuclideanNormAlgorithms.java | 749 +++++++++++++++++++++
.../jmh/arrays/EuclideanNormEvaluator.java | 249 +++++++
.../examples/jmh/arrays/NormsPerformance.java | 198 ++++++
.../jmh/arrays/EuclideanNormAccuracyTest.java | 141 ++++
pom.xml | 9 +
.../resources/spotbugs/spotbugs-exclude-filter.xml | 2 +-
21 files changed, 3303 insertions(+), 201 deletions(-)
diff --git a/.travis.yml b/.travis.yml
index c51bd7d..783355f 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -18,7 +18,7 @@ cache:
directories:
- $HOME/.m2
jdk:
- - openjdk8
+ - openjdk9
# The default build is so verbose the travis log file can exceed the allowed limit.
# Run a custom install that hides all the Maven ‘Downloading’ messages.
diff --git a/commons-numbers-arrays/pom.xml b/commons-numbers-arrays/pom.xml
index d27922a..7dcc5d5 100644
--- a/commons-numbers-arrays/pom.xml
+++ b/commons-numbers-arrays/pom.xml
@@ -32,6 +32,10 @@
<description>Array utilities.</description>
<properties>
+ <!-- JDK 9+ required for BigDecimal.sqrt() method in unit tests. -->
+ <maven.compiler.testSource>1.9</maven.compiler.testSource>
+ <maven.compiler.testTarget>1.9</maven.compiler.testTarget>
+
<!-- The Java Module System Name -->
<commons.module.name>org.apache.commons.numbers.arrays</commons.module.name>
<!-- This value must reflect the current name of the base package. -->
diff --git a/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/CosAngle.java b/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/CosAngle.java
index e728d21..39d0a4f 100644
--- a/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/CosAngle.java
+++ b/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/CosAngle.java
@@ -36,7 +36,7 @@ public final class CosAngle {
*/
public static double value(double[] v1,
double[] v2) {
- return LinearCombination.value(v1, v2) / SafeNorm.value(v1) / SafeNorm.value(v2);
+ return LinearCombination.value(v1, v2) / Norms.euclidean(v1) / Norms.euclidean(v2);
}
}
diff --git a/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/ExtendedPrecision.java b/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/ExtendedPrecision.java
index da7a72d..18854e5 100644
--- a/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/ExtendedPrecision.java
+++ b/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/ExtendedPrecision.java
@@ -206,6 +206,27 @@ final class ExtendedPrecision {
/**
* Compute the low part of the double length number {@code (z,zz)} for the exact
+ * square of {@code x} using Dekker's mult12 algorithm. The standard precision product
+ * {@code x*x} must be provided. The number {@code x} is split into high and low parts
+ * using Dekker's algorithm.
+ *
+ * <p>Warning: This method does not perform scaling in Dekker's split and large
+ * finite numbers can create NaN results.
+ *
+ * @param x Number to square
+ * @param xx Standard precision product {@code x*x}
+ * @return the low part of the square double length number
+ */
+ static double squareLowUnscaled(double x, double xx) {
+ // Split the numbers using Dekker's algorithm without scaling
+ final double hx = highPartUnscaled(x);
+ final double lx = x - hx;
+
+ return productLow(hx, lx, hx, lx, xx);
+ }
+
+ /**
+ * Compute the low part of the double length number {@code (z,zz)} for the exact
* product of {@code x} and {@code y} using Dekker's mult12 algorithm. The standard
* precision product {@code x*y} must be provided. The numbers {@code x} and {@code y}
* should already be split into low and high parts.
diff --git a/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Norms.java b/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Norms.java
new file mode 100644
index 0000000..1e43ec9
--- /dev/null
+++ b/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Norms.java
@@ -0,0 +1,448 @@
+/*
+ * 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.arrays;
+
+/** Class providing methods to compute various norm values.
+ *
+ * <p>This class uses a variety of techniques to increase numerical accuracy
+ * and reduce errors. A primary source for the included algorithms is 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>.
+ * @see <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)">Norm</a>
+ */
+public final class Norms {
+
+ /** Threshold for scaling small numbers. This value is chosen such that doubles
+ * set to this value can be squared without underflow. Values less than this must
+ * be scaled up.
+ */
+ private static final double SMALL_THRESH = 0x1.0p-511;
+
+ /** Threshold for scaling large numbers. This value is chosen such that 2^31 doubles
+ * set to this value can be squared and added without overflow. Values greater than
+ * this must be scaled down.
+ */
+ private static final double LARGE_THRESH = 0x1.0p+496;
+
+ /** Threshold for scaling up a single value by {@link #SCALE_UP} without risking overflow
+ * when the value is squared.
+ */
+ private static final double SAFE_SCALE_UP_THRESH = 0x1.0p-100;
+
+ /** Value used to scale down large numbers. */
+ private static final double SCALE_DOWN = 0x1.0p-600;
+
+ /** Value used to scale up small numbers. */
+ private static final double SCALE_UP = 0x1.0p+600;
+
+ /** Threshold for the difference between the exponents of two Euclidean 2D input values
+ * where the larger value dominates the calculation.
+ */
+ private static final int EXP_DIFF_THRESHOLD_2D = 54;
+
+ /** Utility class; no instantiation. */
+ private Norms() {}
+
+ /** Compute the Manhattan norm (also known as the Taxicab norm or L1 norm) of the arguments.
+ * The result is equal to \(|x| + |y|\), i.e., the sum of the absolute values of the arguments.
+ *
+ * <p>Special cases:
+ * <ul>
+ * <li>If either value is NaN, then the result is NaN.</li>
+ * <li>If either value is infinite and the other value is not NaN, then the result is positive infinity.</li>
+ * </ul>
+ * @param x first input value
+ * @param y second input value
+ * @return Manhattan norm
+ * @see <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Taxicab_norm_or_Manhattan_norm">Manhattan norm</a>
+ */
+ public static double manhattan(final double x, final double y) {
+ return Math.abs(x) + Math.abs(y);
+ }
+
+ /** Compute the Manhattan norm (also known as the Taxicab norm or L1 norm) of the arguments.
+ * The result is equal to \(|x| + |y| + |z|\), i.e., the sum of the absolute values of the arguments.
+ *
+ * <p>Special cases:
+ * <ul>
+ * <li>If any value is NaN, then the result is NaN.</li>
+ * <li>If any value is infinite and no value is NaN, then the result is positive infinity.</li>
+ * </ul>
+ * @param x first input value
+ * @param y second input value
+ * @param z third input value
+ * @return Manhattan norm
+ * @see <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Taxicab_norm_or_Manhattan_norm">Manhattan norm</a>
+ */
+ public static double manhattan(final double x, final double y, final double z) {
+ return Summation.value(
+ Math.abs(x),
+ Math.abs(y),
+ Math.abs(z));
+ }
+
+ /** Compute the Manhattan norm (also known as the Taxicab norm or L1 norm) of the given values.
+ * The result is equal to \(|v_0| + ... + |v_i|\), i.e., the sum of the absolute values of the input elements.
+ *
+ * <p>Special cases:
+ * <ul>
+ * <li>If any value is NaN, then the result is NaN.</li>
+ * <li>If any value is infinite and no value is NaN, then the result is positive infinity.</li>
+ * <li>If the array is empty, then the result is 0.</li>
+ * </ul>
+ * @param v input values
+ * @return Manhattan norm
+ * @see <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Taxicab_norm_or_Manhattan_norm">Manhattan norm</a>
+ */
+ public static double manhattan(final double[] v) {
+ double sum = 0d;
+ double comp = 0d;
+
+ 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;
+ }
+
+ return Summation.summationResult(sum, comp);
+ }
+
+ /** Compute the Euclidean norm (also known as the L2 norm) of the arguments. The result is equal to
+ * \(\sqrt{x^2 + y^2}\). This method correctly handles the possibility of overflow or underflow
+ * during the computation.
+ *
+ * <p>Special cases:
+ * <ul>
+ * <li>If either value is NaN, then the result is NaN.</li>
+ * <li>If either value is infinite and the other value is not NaN, then the result is positive infinity.</li>
+ * </ul>
+ *
+ * <p><strong>Comparison with Math.hypot()</strong>
+ * <p>While not guaranteed to return the same result, this method does provide similar error bounds to
+ * the JDK's Math.hypot() method and may run faster on some JVMs.
+ * @param x first input
+ * @param y second input
+ * @return Euclidean norm
+ * @see <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Euclidean_norm">Euclidean norm</a>
+ */
+ public static double euclidean(final double x, final double y) {
+ final double xabs = Math.abs(x);
+ final double yabs = Math.abs(y);
+
+ final double max;
+ final double min;
+ // the compare method considers NaN greater than other values, meaning that our
+ // check for if the max is finite later on will detect NaNs correctly
+ if (Double.compare(xabs, yabs) > 0) {
+ max = xabs;
+ min = yabs;
+ } else {
+ max = yabs;
+ min = xabs;
+ }
+
+ // if the max is not finite, then one of the inputs must not have
+ // been finite
+ if (!Double.isFinite(max)) {
+ // let the standard multiply operation determine whether to return NaN or infinite
+ return xabs * yabs;
+ } else if (Math.getExponent(max) - Math.getExponent(min) > EXP_DIFF_THRESHOLD_2D) {
+ // value is completely dominated by max; just return max
+ return max;
+ }
+
+ // compute the scale and rescale values
+ final double scale;
+ final double rescale;
+ if (max > LARGE_THRESH) {
+ scale = SCALE_DOWN;
+ rescale = SCALE_UP;
+ } else if (max < SAFE_SCALE_UP_THRESH) {
+ scale = SCALE_UP;
+ rescale = SCALE_DOWN;
+ } else {
+ scale = 1d;
+ rescale = 1d;
+ }
+
+ double sum = 0d;
+ double comp = 0d;
+
+ // add scaled x
+ double sx = xabs * scale;
+ final double px = sx * sx;
+ comp += ExtendedPrecision.squareLowUnscaled(sx, px);
+ final double sumPx = sum + px;
+ comp += ExtendedPrecision.twoSumLow(sum, px, sumPx);
+ sum = sumPx;
+
+ // add scaled y
+ double sy = yabs * scale;
+ final double py = sy * sy;
+ comp += ExtendedPrecision.squareLowUnscaled(sy, py);
+ final double sumPy = sum + py;
+ comp += ExtendedPrecision.twoSumLow(sum, py, sumPy);
+ sum = sumPy;
+
+ return Math.sqrt(sum + comp) * rescale;
+ }
+
+ /** Compute the Euclidean norm (also known as the L2 norm) of the arguments. The result is equal to
+ * \(\sqrt{x^2 + y^2 + z^2}\). This method correctly handles the possibility of overflow or underflow
+ * during the computation.
+ *
+ * <p>Special cases:
+ * <ul>
+ * <li>If any value is NaN, then the result is NaN.</li>
+ * <li>If any value is infinite and no value is NaN, then the result is positive infinity.</li>
+ * </ul>
+ * @param x first input
+ * @param y second input
+ * @param z third input
+ * @return Euclidean norm
+ * @see <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Euclidean_norm">Euclidean norm</a>
+ */
+ public static double euclidean(final double x, final double y, final double z) {
+ final double xabs = Math.abs(x);
+ final double yabs = Math.abs(y);
+ final double zabs = Math.abs(z);
+
+ final double max = Math.max(Math.max(xabs, yabs), zabs);
+
+ // if the max is not finite, then one of the inputs must not have
+ // been finite
+ if (!Double.isFinite(max)) {
+ // let the standard multiply operation determine whether to return NaN or infinite
+ return xabs * yabs * zabs;
+ }
+
+ // compute the scale and rescale values
+ final double scale;
+ final double rescale;
+ if (max > LARGE_THRESH) {
+ scale = SCALE_DOWN;
+ rescale = SCALE_UP;
+ } else if (max < SAFE_SCALE_UP_THRESH) {
+ scale = SCALE_UP;
+ rescale = SCALE_DOWN;
+ } else {
+ scale = 1d;
+ rescale = 1d;
+ }
+
+ double sum = 0d;
+ double comp = 0d;
+
+ // add scaled x
+ double sx = xabs * scale;
+ final double px = sx * sx;
+ comp += ExtendedPrecision.squareLowUnscaled(sx, px);
+ final double sumPx = sum + px;
+ comp += ExtendedPrecision.twoSumLow(sum, px, sumPx);
+ sum = sumPx;
+
+ // add scaled y
+ double sy = yabs * scale;
+ final double py = sy * sy;
+ comp += ExtendedPrecision.squareLowUnscaled(sy, py);
+ final double sumPy = sum + py;
+ comp += ExtendedPrecision.twoSumLow(sum, py, sumPy);
+ sum = sumPy;
+
+ // add scaled z
+ final double sz = zabs * scale;
+ final double pz = sz * sz;
+ comp += ExtendedPrecision.squareLowUnscaled(sz, pz);
+ final double sumPz = sum + pz;
+ comp += ExtendedPrecision.twoSumLow(sum, pz, sumPz);
+ sum = sumPz;
+
+ return Math.sqrt(sum + comp) * rescale;
+ }
+
+ /** Compute the Euclidean norm (also known as the L2 norm) of the given values. The result is equal to
+ * \(\sqrt{v_0^2 + ... + v_{n-1}^2}\). This method correctly handles the possibility of overflow or underflow
+ * during the computation.
+ *
+ * <p>Special cases:
+ * <ul>
+ * <li>If any value is NaN, then the result is NaN.</li>
+ * <li>If any value is infinite and no value is NaN, then the result is positive infinity.</li>
+ * <li>If the array is empty, then the result is 0.</li>
+ * </ul>
+ * @param v input values
+ * @return Euclidean norm
+ * @see <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Euclidean_norm">Euclidean norm</a>
+ */
+ public static double euclidean(final double[] v) {
+ // sum of big, normal and small numbers
+ double s1 = 0;
+ double s2 = 0;
+ double s3 = 0;
+
+ // sum compensation values
+ double c1 = 0;
+ double c2 = 0;
+ double c3 = 0;
+
+ for (int i = 0; i < v.length; ++i) {
+ final double x = Math.abs(v[i]);
+ if (!Double.isFinite(x)) {
+ // not finite; determine whether to return NaN or positive infinity
+ return euclideanNormSpecial(v, i);
+ } else if (x > LARGE_THRESH) {
+ // scale down
+ final double sx = x * SCALE_DOWN;
+
+ // compute the product and product compensation
+ final double p = sx * sx;
+ final double cp = ExtendedPrecision.squareLowUnscaled(sx, p);
+
+ // compute the running sum and sum compensation
+ final double s = s1 + p;
+ final double cs = ExtendedPrecision.twoSumLow(s1, p, s);
+
+ // update running totals
+ c1 += cp + cs;
+ s1 = s;
+ } else if (x < SMALL_THRESH) {
+ // scale up
+ final double sx = x * SCALE_UP;
+
+ // compute the product and product compensation
+ final double p = sx * sx;
+ final double cp = ExtendedPrecision.squareLowUnscaled(sx, p);
+
+ // compute the running sum and sum compensation
+ final double s = s3 + p;
+ final double cs = ExtendedPrecision.twoSumLow(s3, p, s);
+
+ // update running totals
+ c3 += cp + cs;
+ s3 = s;
+ } else {
+ // no scaling
+ // compute the product and product compensation
+ final double p = x * x;
+ final double cp = ExtendedPrecision.squareLowUnscaled(x, p);
+
+ // compute the running sum and sum compensation
+ final double s = s2 + p;
+ final double cs = ExtendedPrecision.twoSumLow(s2, p, s);
+
+ // update running totals
+ c2 += cp + cs;
+ s2 = s;
+ }
+ }
+
+ // The highest sum is the significant component. Add the next significant.
+ // Note that the "x * SCALE_DOWN * SCALE_DOWN" expressions must be executed
+ // in the order given. If the two scale factors are multiplied together first,
+ // they will underflow to zero.
+ if (s1 != 0) {
+ // add s1, s2, c1, c2
+ final double s2Adj = s2 * SCALE_DOWN * SCALE_DOWN;
+ final double sum = s1 + s2Adj;
+ final double comp = ExtendedPrecision.twoSumLow(s1, s2Adj, sum) + c1 + (c2 * SCALE_DOWN * SCALE_DOWN);
+ return Math.sqrt(sum + comp) * SCALE_UP;
+ } else if (s2 != 0) {
+ // add s2, s3, c2, c3
+ final double s3Adj = s3 * SCALE_DOWN * SCALE_DOWN;
+ final double sum = s2 + s3Adj;
+ final double comp = ExtendedPrecision.twoSumLow(s2, s3Adj, sum) + c2 + (c3 * SCALE_DOWN * SCALE_DOWN);
+ return Math.sqrt(sum + comp);
+ }
+ // add s3, c3
+ return Math.sqrt(s3 + c3) * SCALE_DOWN;
+ }
+
+ /** Return a Euclidean norm value for special cases of non-finite input.
+ * @param v input vector
+ * @param start index to start examining the input vector from
+ * @return Euclidean norm special value
+ */
+ private static double euclideanNormSpecial(final double[] v, final int start) {
+ for (int i = start; i < v.length; ++i) {
+ if (Double.isNaN(v[i])) {
+ return Double.NaN;
+ }
+ }
+ return Double.POSITIVE_INFINITY;
+ }
+
+ /** Compute the maximum norm (also known as the infinity norm or L<sub>inf</sub> norm) of the arguments.
+ * The result is equal to \(\max{(|x|, |y|)}\), i.e., the maximum of the absolute values of the arguments.
+ *
+ * <p>Special cases:
+ * <ul>
+ * <li>If either value is NaN, then the result is NaN.</li>
+ * <li>If either value is infinite and the other value is not NaN, then the result is positive infinity.</li>
+ * </ul>
+ * @param x first input
+ * @param y second input
+ * @return maximum norm
+ * @see <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Maximum_norm_(special_case_of:_infinity_norm,_uniform_norm,_or_supremum_norm)">Maximum norm</a>
+ */
+ public static double maximum(final double x, final double y) {
+ return Math.max(Math.abs(x), Math.abs(y));
+ }
+
+ /** Compute the maximum norm (also known as the infinity norm or L<sub>inf</sub> norm) of the arguments.
+ * The result is equal to \(\max{(|x|, |y|, |z|)}\), i.e., the maximum of the absolute values of the arguments.
+ *
+ * <p>Special cases:
+ * <ul>
+ * <li>If any value is NaN, then the result is NaN.</li>
+ * <li>If any value is infinite and no value is NaN, then the result is positive infinity.</li>
+ * </ul>
+ * @param x first input
+ * @param y second input
+ * @param z third input
+ * @return maximum norm
+ * @see <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Maximum_norm_(special_case_of:_infinity_norm,_uniform_norm,_or_supremum_norm)">Maximum norm</a>
+ */
+ public static double maximum(final double x, final double y, final double z) {
+ return Math.max(
+ Math.abs(x),
+ Math.max(Math.abs(y), Math.abs(z)));
+ }
+
+ /** Compute the maximum norm (also known as the infinity norm or L<sub>inf</sub> norm) of the given values.
+ * The result is equal to \(\max{(|v_0|, \ldots, |v_{n-1}|)}\), i.e., the maximum of the absolute values of the
+ * input elements.
+ *
+ * <p>Special cases:
+ * <ul>
+ * <li>If any value is NaN, then the result is NaN.</li>
+ * <li>If any value is infinite and no value is NaN, then the result is positive infinity.</li>
+ * <li>If the array is empty, then the result is 0.</li>
+ * </ul>
+ * @param v input values
+ * @return maximum norm
+ * @see <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Maximum_norm_(special_case_of:_infinity_norm,_uniform_norm,_or_supremum_norm)">Maximum norm</a>
+ */
+ public static double maximum(final double[] v) {
+ double max = 0d;
+ for (int i = 0; i < v.length; ++i) {
+ max = Math.max(max, Math.abs(v[i]));
+ }
+ return max;
+ }
+}
diff --git a/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/SafeNorm.java b/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/SafeNorm.java
deleted file mode 100644
index 0656bd6..0000000
--- a/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/SafeNorm.java
+++ /dev/null
@@ -1,91 +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.arrays;
-
-/**
- * Computes the Cartesian norm (2-norm), handling both overflow and underflow.
- * Translation of the <a href="http://www.netlib.org/minpack">minpack</a>
- * "enorm" subroutine.
- */
-public final class SafeNorm {
- /** Constant. */
- private static final double R_DWARF = 3.834e-20;
- /** Constant. */
- private static final double R_GIANT = 1.304e+19;
-
- /** Private constructor. */
- private SafeNorm() {
- // intentionally empty.
- }
-
- /**
- * @param v Cartesian coordinates.
- * @return the 2-norm of the vector.
- */
- public static double value(double[] v) {
- double s1 = 0;
- double s2 = 0;
- double s3 = 0;
- double x1max = 0;
- double x3max = 0;
- final double floatn = v.length;
- final double agiant = R_GIANT / floatn;
- for (int i = 0; i < v.length; i++) {
- final double xabs = Math.abs(v[i]);
- if (xabs < R_DWARF || xabs > agiant) {
- if (xabs > R_DWARF) {
- if (xabs > x1max) {
- final double r = x1max / xabs;
- s1 = 1 + s1 * r * r;
- x1max = xabs;
- } else {
- final double r = xabs / x1max;
- s1 += r * r;
- }
- } else {
- if (xabs > x3max) {
- final double r = x3max / xabs;
- s3 = 1 + s3 * r * r;
- x3max = xabs;
- } else {
- if (xabs != 0) {
- final double r = xabs / x3max;
- s3 += r * r;
- }
- }
- }
- } else {
- s2 += xabs * xabs;
- }
- }
- double norm;
- if (s1 != 0) {
- norm = x1max * Math.sqrt(s1 + (s2 / x1max) / x1max);
- } else {
- if (s2 == 0) {
- norm = x3max * Math.sqrt(s3);
- } else {
- if (s2 >= x3max) {
- norm = Math.sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));
- } else {
- norm = Math.sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
- }
- }
- }
- return norm;
- }
-}
diff --git a/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Summation.java b/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Summation.java
new file mode 100644
index 0000000..da1bec0
--- /dev/null
+++ b/commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Summation.java
@@ -0,0 +1,179 @@
+/*
+ * 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.arrays;
+
+/** 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-arrays/src/test/java/org/apache/commons/numbers/arrays/DoubleTestUtils.java b/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/DoubleTestUtils.java
new file mode 100644
index 0000000..3735434
--- /dev/null
+++ b/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/DoubleTestUtils.java
@@ -0,0 +1,69 @@
+/*
+ * 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.arrays;
+
+import org.apache.commons.rng.UniformRandomProvider;
+
+/** Class providing test utilities related to doubles.
+ */
+final class DoubleTestUtils {
+
+ /** Utility class; no instantiation. */
+ private DoubleTestUtils() {}
+
+ /** Compute the difference in ULP between two arguments of the same sign.
+ * @param a first argument
+ * @param b second argument
+ * @return ULP difference between the arguments
+ */
+ public static int computeUlpDifference(final double a, final double b) {
+ return (int) (Double.doubleToLongBits(a) - Double.doubleToLongBits(b));
+ }
+
+ /** Construct an array of length {@code len} containing random double values with exponents between
+ * {@code minExp} and {@code maxExp}.
+ * @param len vector length
+ * @param minExp minimum element exponent value
+ * @param maxExp maximum element exponent value
+ * @param rng random number generator
+ * @return random vector array
+ */
+ public static double[] randomArray(final int len, final int minExp, final int maxExp,
+ final UniformRandomProvider rng) {
+ final double[] v = new double[len];
+ for (int i = 0; i < v.length; ++i) {
+ v[i] = randomDouble(minExp, maxExp, rng);
+ }
+ return v;
+ }
+
+ /** Construct a random double with an exponent in the range {@code [minExp, maxExp]}.
+ * @param minExp minimum exponent; must be less than {@code maxExp}
+ * @param maxExp maximum exponent; must be greater than {@code minExp}
+ * @param rng random number generator
+ * @return random double value with an exponent in the specified range
+ */
+ public static double randomDouble(final int minExp, final int maxExp, final UniformRandomProvider rng) {
+ // Create random doubles using random bits in the sign bit and the mantissa.
+ final long mask = ((1L << 52) - 1) | 1L << 63;
+ final long bits = rng.nextLong() & mask;
+ // The exponent must be unsigned so + 1023 to the signed exponent
+ final int expRange = maxExp - minExp + 1;
+ final long exp = rng.nextInt(expRange) + minExp + 1023;
+ return Double.longBitsToDouble(bits | (exp << 52));
+ }
+}
diff --git a/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/ExtendedPrecisionTest.java b/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/ExtendedPrecisionTest.java
index 262d6c6..c9647fb 100644
--- a/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/ExtendedPrecisionTest.java
+++ b/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/ExtendedPrecisionTest.java
@@ -102,6 +102,23 @@ class ExtendedPrecisionTest {
final double lo2 = a - hi2;
Assertions.assertEquals(a, hi2);
Assertions.assertEquals(0, lo2);
+
Assertions.assertTrue(Math.abs(hi2) > Math.abs(lo2));
}
+
+ @Test
+ void testSquareLowUnscaled() {
+ assertSquareLowUnscaled(0.0, 1.0);
+ assertSquareLowUnscaled(0.0, -1.0);
+ assertSquareLowUnscaled(Math.fma(Math.PI, Math.PI, -Math.PI * Math.PI), Math.PI);
+
+ assertSquareLowUnscaled(Double.NaN, Double.POSITIVE_INFINITY);
+ assertSquareLowUnscaled(Double.NaN, Double.NEGATIVE_INFINITY);
+ assertSquareLowUnscaled(Double.NaN, Double.NaN);
+ assertSquareLowUnscaled(Double.NaN, Double.MAX_VALUE);
+ }
+
+ private static void assertSquareLowUnscaled(final double expected, final double x) {
+ Assertions.assertEquals(expected, ExtendedPrecision.squareLowUnscaled(x, x * x));
+ }
}
diff --git a/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/NormsTest.java b/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/NormsTest.java
new file mode 100644
index 0000000..acb5610
--- /dev/null
+++ b/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/NormsTest.java
@@ -0,0 +1,597 @@
+/*
+ * 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.arrays;
+
+import java.math.BigDecimal;
+import java.math.MathContext;
+import java.util.Arrays;
+import java.util.function.ToDoubleFunction;
+
+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;
+
+class NormsTest {
+
+ private static final int SMALL_THRESH_EXP = -511;
+
+ private static final int LARGE_THRESH_EXP = +496;
+
+ private static final int RAND_VECTOR_CNT = 1_000;
+
+ private static final int MAX_ULP_ERR = 1;
+
+ private static final double HYPOT_COMPARE_EPS = 1e-2;
+
+ @Test
+ void testManhattan_2d() {
+ // act/assert
+ Assertions.assertEquals(0d, Norms.manhattan(0d, -0d));
+ Assertions.assertEquals(3d, Norms.manhattan(-1d, 2d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, Norms.manhattan(Double.MAX_VALUE, Double.MAX_VALUE));
+
+ Assertions.assertEquals(Double.NaN, Norms.manhattan(Double.NaN, 1d));
+ Assertions.assertEquals(Double.NaN, Norms.manhattan(1d, Double.NaN));
+ Assertions.assertEquals(Double.NaN, Norms.manhattan(Double.POSITIVE_INFINITY, Double.NaN));
+
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.manhattan(Double.POSITIVE_INFINITY, 0d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.manhattan(0d, Double.POSITIVE_INFINITY));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.manhattan(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY));
+ }
+
+ @Test
+ void testManhattan_3d() {
+ // act/assert
+ Assertions.assertEquals(0d, Norms.manhattan(0d, -0d, 0d));
+ Assertions.assertEquals(6d, Norms.manhattan(-1d, 2d, -3d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, Norms.manhattan(Double.MAX_VALUE, Double.MAX_VALUE, 0d));
+
+ Assertions.assertEquals(Double.NaN, Norms.manhattan(Double.NaN, -2d, 1d));
+ Assertions.assertEquals(Double.NaN, Norms.manhattan(-2d, Double.NaN, 1d));
+ Assertions.assertEquals(Double.NaN, Norms.manhattan(-2d, 1d, Double.NaN));
+ Assertions.assertEquals(Double.NaN, Norms.manhattan(-2d, Double.POSITIVE_INFINITY, Double.NaN));
+
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.manhattan(Double.POSITIVE_INFINITY, 2d, -4d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.manhattan(Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, -4d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.manhattan(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY));
+ }
+
+ @Test
+ void testManhattan_array() {
+ // act/assert
+ Assertions.assertEquals(0d, Norms.manhattan(new double[0]));
+ Assertions.assertEquals(0d, Norms.manhattan(new double[] {0d, -0d}));
+ Assertions.assertEquals(6d, Norms.manhattan(new double[] {-1d, 2d, -3d}));
+ Assertions.assertEquals(10d, Norms.manhattan(new double[] {-1d, 2d, -3d, 4d}));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.manhattan(new double[] {Double.MAX_VALUE, Double.MAX_VALUE}));
+
+ Assertions.assertEquals(Double.NaN, Norms.manhattan(new double[] {-2d, Double.NaN, 1d}));
+ Assertions.assertEquals(Double.NaN, Norms.manhattan(new double[] {Double.POSITIVE_INFINITY, Double.NaN, 1d}));
+
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.manhattan(new double[] {Double.POSITIVE_INFINITY, 0d}));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.manhattan(new double[] {Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY}));
+ }
+
+ @Test
+ void testEuclidean_2d_simple() {
+ // act/assert
+ Assertions.assertEquals(0d, Norms.euclidean(0d, 0d));
+ Assertions.assertEquals(1d, Norms.euclidean(1d, 0d));
+ Assertions.assertEquals(1d, Norms.euclidean(0d, 1d));
+ Assertions.assertEquals(5d, Norms.euclidean(-3d, 4d));
+ Assertions.assertEquals(Double.MIN_VALUE, Norms.euclidean(0d, Double.MIN_VALUE));
+ Assertions.assertEquals(Double.MAX_VALUE, Norms.euclidean(Double.MAX_VALUE, 0d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, Norms.euclidean(Double.MAX_VALUE, Double.MAX_VALUE));
+
+ Assertions.assertEquals(Math.sqrt(2), Norms.euclidean(1d, -1d));
+
+ Assertions.assertEquals(Double.NaN, Norms.euclidean(Double.NaN, -2d));
+ Assertions.assertEquals(Double.NaN, Norms.euclidean(Double.NaN, Double.POSITIVE_INFINITY));
+ Assertions.assertEquals(Double.NaN, Norms.euclidean(-2d, Double.NaN));
+ Assertions.assertEquals(Double.NaN,
+ Norms.euclidean(Double.NaN, Double.NEGATIVE_INFINITY));
+
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.euclidean(1d, Double.NEGATIVE_INFINITY));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.euclidean(Double.POSITIVE_INFINITY, -1d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.euclidean(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY));
+ }
+
+ @Test
+ void testEuclidean_2d_scaled() {
+ // arrange
+ final double[] ones = new double[] {1, 1};
+ final double[] multiplesOfTen = new double[] {1, 10};
+ final ToDoubleFunction<double[]> fn = v -> Norms.euclidean(v[0], v[1]);
+
+ // act/assert
+ checkScaledEuclideanNorm(ones, 0, fn);
+ checkScaledEuclideanNorm(ones, LARGE_THRESH_EXP, fn);
+ checkScaledEuclideanNorm(ones, LARGE_THRESH_EXP + 1, fn);
+ checkScaledEuclideanNorm(ones, -100, fn);
+ checkScaledEuclideanNorm(ones, -101, fn);
+ checkScaledEuclideanNorm(ones, SMALL_THRESH_EXP, fn);
+ checkScaledEuclideanNorm(ones, SMALL_THRESH_EXP - 1, fn);
+
+
+ checkScaledEuclideanNorm(multiplesOfTen, 0, fn);
+ checkScaledEuclideanNorm(multiplesOfTen, -100, fn);
+ checkScaledEuclideanNorm(multiplesOfTen, -101, fn);
+ checkScaledEuclideanNorm(multiplesOfTen, LARGE_THRESH_EXP - 1, fn);
+ checkScaledEuclideanNorm(multiplesOfTen, SMALL_THRESH_EXP, fn);
+ }
+
+ @Test
+ void testEuclidean_2d_dominantValue() {
+ // act/assert
+ Assertions.assertEquals(Math.PI, Norms.euclidean(-Math.PI, 0x1.0p-55));
+ Assertions.assertEquals(Math.PI, Norms.euclidean(0x1.0p-55, -Math.PI));
+ }
+
+ @Test
+ void testEuclidean_2d_random() {
+ // arrange
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP, 1L);
+
+ // act/assert
+ checkEuclideanRandom(2, rng, v -> Norms.euclidean(v[0], v[1]));
+ }
+
+ @Test
+ void testEuclidean_2d_vsArray() {
+ // arrange
+ final double[][] inputs = {
+ {-4.074598908124454E-9, 9.897869969944898E-28},
+ {1.3472131556526359E-27, -9.064577177323565E9},
+ {-3.9219339341360245E149, -7.132522817112096E148},
+ {-1.4888098520466735E153, -2.9099184907796666E150},
+ {-8.659395144898396E-152, -1.123275532302136E-150},
+ {-3.660198254902351E-152, -6.656524053354807E-153}
+ };
+
+ // act/assert
+ for (final double[] input : inputs) {
+ Assertions.assertEquals(Norms.euclidean(input), Norms.euclidean(input[0], input[1]),
+ () -> "Expected inline method result to equal array result for input " + Arrays.toString(input));
+ }
+ }
+
+ @Test
+ void testEuclidean_2d_vsHypot() {
+ // arrange
+ final int samples = 1000;
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP, 2L);
+
+ // act/assert
+ assertEuclidean2dVersusHypot(-10, +10, samples, rng);
+ assertEuclidean2dVersusHypot(0, +20, samples, rng);
+ assertEuclidean2dVersusHypot(-20, 0, samples, rng);
+ assertEuclidean2dVersusHypot(-20, +20, samples, rng);
+ assertEuclidean2dVersusHypot(-100, +100, samples, rng);
+ assertEuclidean2dVersusHypot(LARGE_THRESH_EXP - 10, LARGE_THRESH_EXP + 10, samples, rng);
+ assertEuclidean2dVersusHypot(SMALL_THRESH_EXP - 10, SMALL_THRESH_EXP + 10, samples, rng);
+ assertEuclidean2dVersusHypot(-600, +600, samples, rng);
+ }
+
+ /** Assert that the Norms euclidean 2D computation produces similar error behavior to Math.hypot().
+ * @param minExp minimum exponent for random inputs
+ * @param maxExp maximum exponent for random inputs
+ * @param samples sample count
+ * @param rng random number generator
+ */
+ private static void assertEuclidean2dVersusHypot(final int minExp, final int maxExp, final int samples,
+ final UniformRandomProvider rng) {
+ // generate random inputs
+ final double[][] inputs = new double[samples][];
+ for (int i = 0; i < samples; ++i) {
+ inputs[i] = DoubleTestUtils.randomArray(2, minExp, maxExp, rng);
+ }
+
+ // compute exact results
+ final double[] exactResults = new double[samples];
+ for (int i = 0; i < samples; ++i) {
+ exactResults[i] = exactEuclideanNorm(inputs[i]);
+ }
+
+ // compute the std devs
+ final UlpErrorStats hypotStats = computeUlpErrorStats(inputs, exactResults, v -> Math.hypot(v[0], v[1]));
+ final UlpErrorStats normStats = computeUlpErrorStats(inputs, exactResults, v -> Norms.euclidean(v[0], v[1]));
+
+ // ensure that we are within the ballpark of Math.hypot
+ Assertions.assertTrue(normStats.getMean() <= (hypotStats.getMean() + HYPOT_COMPARE_EPS),
+ () -> "Expected 2D norm result to have similar error mean to Math.hypot(): hypot error mean= " +
+ hypotStats.getMean() + ", norm error mean= " + normStats.getMean());
+
+ Assertions.assertTrue(normStats.getStdDev() <= (hypotStats.getStdDev() + HYPOT_COMPARE_EPS),
+ () -> "Expected 2D norm result to have similar std deviation to Math.hypot(): hypot std dev= " +
+ hypotStats.getStdDev() + ", norm std dev= " + normStats.getStdDev());
+ }
+
+ @Test
+ void testEuclidean_3d_simple() {
+ // act/assert
+ Assertions.assertEquals(0d, Norms.euclidean(0d, 0d, 0d));
+ Assertions.assertEquals(1d, Norms.euclidean(1d, 0d, 0d));
+ Assertions.assertEquals(1d, Norms.euclidean(0d, 1d, 0d));
+ Assertions.assertEquals(1d, Norms.euclidean(0d, 0d, 1d));
+ Assertions.assertEquals(5 * Math.sqrt(2), Norms.euclidean(-3d, -4d, 5d));
+ Assertions.assertEquals(Double.MIN_VALUE, Norms.euclidean(0d, 0d, Double.MIN_VALUE));
+ Assertions.assertEquals(Double.MAX_VALUE, Norms.euclidean(Double.MAX_VALUE, 0d, 0d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.euclidean(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE));
+
+ Assertions.assertEquals(Math.sqrt(3), Norms.euclidean(1d, -1d, 1d));
+
+ Assertions.assertEquals(Double.NaN, Norms.euclidean(Double.NaN, -2d, 0d));
+ Assertions.assertEquals(Double.NaN, Norms.euclidean(-2d, Double.NaN, 0d));
+ Assertions.assertEquals(Double.NaN, Norms.euclidean(-2d, 0d, Double.NaN));
+ Assertions.assertEquals(Double.NaN,
+ Norms.euclidean(Double.POSITIVE_INFINITY, Double.NaN, 1d));
+
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.euclidean(Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, 1d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.euclidean(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY));
+ }
+
+ @Test
+ void testEuclidean_3d_scaled() {
+ // arrange
+ final double[] ones = new double[] {1, 1, 1};
+ final double[] multiplesOfTen = new double[] {1, 10, 100};
+ final ToDoubleFunction<double[]> fn = v -> Norms.euclidean(v[0], v[1], v[2]);
+
+ // act/assert
+ checkScaledEuclideanNorm(ones, 0, fn);
+ checkScaledEuclideanNorm(ones, LARGE_THRESH_EXP, fn);
+ checkScaledEuclideanNorm(ones, LARGE_THRESH_EXP + 1, fn);
+ checkScaledEuclideanNorm(ones, -100, fn);
+ checkScaledEuclideanNorm(ones, -101, fn);
+ checkScaledEuclideanNorm(ones, SMALL_THRESH_EXP, fn);
+ checkScaledEuclideanNorm(ones, SMALL_THRESH_EXP - 1, fn);
+
+ checkScaledEuclideanNorm(multiplesOfTen, 0, fn);
+ checkScaledEuclideanNorm(multiplesOfTen, -100, fn);
+ checkScaledEuclideanNorm(multiplesOfTen, -101, fn);
+ checkScaledEuclideanNorm(multiplesOfTen, LARGE_THRESH_EXP - 1, fn);
+ checkScaledEuclideanNorm(multiplesOfTen, SMALL_THRESH_EXP - 1, fn);
+ }
+
+ @Test
+ void testEuclidean_3d_random() {
+ // arrange
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP, 1L);
+
+ // act/assert
+ checkEuclideanRandom(3, rng, v -> Norms.euclidean(v[0], v[1], v[2]));
+ }
+
+ @Test
+ void testEuclidean_3d_vsArray() {
+ // arrange
+ final double[][] inputs = {
+ {-4.074598908124454E-9, 9.897869969944898E-28, 7.849935157082846E-14},
+ {1.3472131556526359E-27, -9.064577177323565E9, 323771.526282239},
+ {-3.9219339341360245E149, -7.132522817112096E148, -3.427334456813165E147},
+ {-1.4888098520466735E153, -2.9099184907796666E150, 1.0144962310234785E152},
+ {-8.659395144898396E-152, -1.123275532302136E-150, -2.151505326692001E-152},
+ {-3.660198254902351E-152, -6.656524053354807E-153, -3.198606556986218E-154}
+ };
+
+ // act/assert
+ for (final double[] input : inputs) {
+ Assertions.assertEquals(Norms.euclidean(input), Norms.euclidean(input[0], input[1], input[2]),
+ () -> "Expected inline method result to equal array result for input " + Arrays.toString(input));
+ }
+ }
+
+ @Test
+ void testEuclidean_array_simple() {
+ // act/assert
+ Assertions.assertEquals(0d, Norms.euclidean(new double[0]));
+ Assertions.assertEquals(5d, Norms.euclidean(new double[] {-3d, 4d}));
+
+ Assertions.assertEquals(Math.sqrt(2), Norms.euclidean(new double[] {1d, -1d}));
+ Assertions.assertEquals(Math.sqrt(3), Norms.euclidean(new double[] {1d, -1d, 1d}));
+ Assertions.assertEquals(2, Norms.euclidean(new double[] {1d, -1d, 1d, -1d}));
+
+ final double[] longVec = new double[] {-0.9, 8.7, -6.5, -4.3, -2.1, 0, 1.2, 3.4, -5.6, 7.8, 9.0};
+ Assertions.assertEquals(directEuclideanNorm(longVec), Norms.euclidean(longVec));
+
+ Assertions.assertEquals(Double.MIN_VALUE, Norms.euclidean(new double[] {0d, Double.MIN_VALUE}));
+ Assertions.assertEquals(Double.MAX_VALUE, Norms.euclidean(new double[] {Double.MAX_VALUE, 0d}));
+
+ final double[] maxVec = new double[1000];
+ Arrays.fill(maxVec, Double.MAX_VALUE);
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, Norms.euclidean(maxVec));
+
+ final double[] largeThreshVec = new double[1000];
+ Arrays.fill(largeThreshVec, 0x1.0p496);
+ Assertions.assertEquals(Math.sqrt(largeThreshVec.length) * largeThreshVec[0], Norms.euclidean(largeThreshVec));
+
+ Assertions.assertEquals(Double.NaN, Norms.euclidean(new double[] {-2d, Double.NaN, 1d}));
+ Assertions.assertEquals(Double.NaN,
+ Norms.euclidean(new double[] {Double.POSITIVE_INFINITY, Double.NaN}));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.euclidean(new double[] {Double.POSITIVE_INFINITY, 1, 0}));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.euclidean(new double[] {Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY}));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.euclidean(new double[] {Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY}));
+ }
+
+ @Test
+ void testEuclidean_array_scaled() {
+ // arrange
+ final double[] ones = new double[] {1, 1, 1, 1};
+ final double[] multiplesOfTen = new double[] {1, 10, 100, 1000};
+
+ // act/assert
+ checkScaledEuclideanNorm(ones, 0, Norms::euclidean);
+ checkScaledEuclideanNorm(ones, LARGE_THRESH_EXP, Norms::euclidean);
+ checkScaledEuclideanNorm(ones, LARGE_THRESH_EXP + 1, Norms::euclidean);
+ checkScaledEuclideanNorm(ones, SMALL_THRESH_EXP, Norms::euclidean);
+ checkScaledEuclideanNorm(ones, SMALL_THRESH_EXP - 1, Norms::euclidean);
+
+ checkScaledEuclideanNorm(multiplesOfTen, 1, Norms::euclidean);
+ checkScaledEuclideanNorm(multiplesOfTen, LARGE_THRESH_EXP - 1, Norms::euclidean);
+ checkScaledEuclideanNorm(multiplesOfTen, SMALL_THRESH_EXP - 1, Norms::euclidean);
+ }
+
+ @Test
+ void testEuclidean_array_random() {
+ // arrange
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP, 1L);
+
+ // act/assert
+ checkEuclideanRandom(2, rng, Norms::euclidean);
+ checkEuclideanRandom(3, rng, Norms::euclidean);
+ checkEuclideanRandom(4, rng, Norms::euclidean);
+ checkEuclideanRandom(10, rng, Norms::euclidean);
+ checkEuclideanRandom(100, rng, Norms::euclidean);
+ }
+
+ @Test
+ void testMaximum_2d() {
+ // act/assert
+ Assertions.assertEquals(0d, Norms.maximum(0d, -0d));
+ Assertions.assertEquals(2d, Norms.maximum(1d, -2d));
+ Assertions.assertEquals(3d, Norms.maximum(3d, 1d));
+ Assertions.assertEquals(Double.MAX_VALUE, Norms.maximum(Double.MAX_VALUE, Double.MAX_VALUE));
+
+ Assertions.assertEquals(Double.NaN, Norms.maximum(Double.NaN, 0d));
+ Assertions.assertEquals(Double.NaN, Norms.maximum(0d, Double.NaN));
+ Assertions.assertEquals(Double.NaN, Norms.maximum(Double.POSITIVE_INFINITY, Double.NaN));
+
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, Norms.maximum(Double.POSITIVE_INFINITY, 0d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, Norms.maximum(Double.NEGATIVE_INFINITY, 0d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, Norms.maximum(0d, Double.NEGATIVE_INFINITY));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.maximum(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY));
+ }
+
+ @Test
+ void testMaximum_3d() {
+ // act/assert
+ Assertions.assertEquals(0d, Norms.maximum(0d, -0d, 0d));
+ Assertions.assertEquals(3d, Norms.maximum(1d, -2d, 3d));
+ Assertions.assertEquals(4d, Norms.maximum(-4d, -2d, 3d));
+ Assertions.assertEquals(Double.MAX_VALUE, Norms.maximum(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE));
+
+ Assertions.assertEquals(Double.NaN, Norms.maximum(Double.NaN, 3d, 0d));
+ Assertions.assertEquals(Double.NaN, Norms.maximum(3d, Double.NaN, 0d));
+ Assertions.assertEquals(Double.NaN, Norms.maximum(3d, 0d, Double.NaN));
+ Assertions.assertEquals(Double.NaN, Norms.maximum(Double.POSITIVE_INFINITY, 0d, Double.NaN));
+
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, Norms.maximum(Double.POSITIVE_INFINITY, 0d, 1d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, Norms.maximum(0d, Double.POSITIVE_INFINITY, 1d));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, Norms.maximum(0d, 1d, Double.NEGATIVE_INFINITY));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.maximum(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY));
+ }
+
+ @Test
+ void testMaximum_array() {
+ // act/assert
+ Assertions.assertEquals(0d, Norms.maximum(new double[0]));
+ Assertions.assertEquals(0d, Norms.maximum(new double[] {0d, -0d}));
+ Assertions.assertEquals(3d, Norms.maximum(new double[] {-1d, 2d, -3d}));
+ Assertions.assertEquals(4d, Norms.maximum(new double[] {-1d, 2d, -3d, 4d}));
+ Assertions.assertEquals(Double.MAX_VALUE, Norms.maximum(new double[] {Double.MAX_VALUE, Double.MAX_VALUE}));
+
+ Assertions.assertEquals(Double.NaN, Norms.maximum(new double[] {-2d, Double.NaN, 1d}));
+ Assertions.assertEquals(Double.NaN, Norms.maximum(new double[] {Double.POSITIVE_INFINITY, Double.NaN}));
+
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.maximum(new double[] {0d, Double.POSITIVE_INFINITY}));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY,
+ Norms.maximum(new double[] {Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY}));
+ }
+
+ /** Check a number of random vectors of length {@code len} with various exponent
+ * ranges.
+ * @param len vector array length
+ * @param rng random number generator
+ * @param fn euclidean norm test function
+ */
+ private static void checkEuclideanRandom(final int len, final UniformRandomProvider rng,
+ final ToDoubleFunction<double[]> fn) {
+ checkEuclideanRandom(len, +600, +620, rng, fn);
+ checkEuclideanRandom(len, LARGE_THRESH_EXP - 10, LARGE_THRESH_EXP + 10, rng, fn);
+ checkEuclideanRandom(len, +400, +420, rng, fn);
+ checkEuclideanRandom(len, +100, +120, rng, fn);
+ checkEuclideanRandom(len, -10, +10, rng, fn);
+ checkEuclideanRandom(len, -120, -100, rng, fn);
+ checkEuclideanRandom(len, -420, -400, rng, fn);
+ checkEuclideanRandom(len, SMALL_THRESH_EXP - 10, SMALL_THRESH_EXP + 10, rng, fn);
+ checkEuclideanRandom(len, -620, -600, rng, fn);
+
+ checkEuclideanRandom(len, -600, +600, rng, fn);
+ }
+
+ /** Check a number of random vectors of length {@code len} with elements containing
+ * exponents in the range {@code [minExp, maxExp]}.
+ * @param len vector array length
+ * @param minExp min exponent
+ * @param maxExp max exponent
+ * @param rng random number generator
+ * @param fn euclidean norm test function
+ */
+ private static void checkEuclideanRandom(final int len, final int minExp, final int maxExp,
+ final UniformRandomProvider rng, final ToDoubleFunction<double[]> fn) {
+ for (int i = 0; i < RAND_VECTOR_CNT; ++i) {
+ // arrange
+ final double[] v = DoubleTestUtils.randomArray(len, minExp, maxExp, rng);
+
+ final double exact = exactEuclideanNorm(v);
+ final double direct = directEuclideanNorm(v);
+
+ // act
+ final double actual = fn.applyAsDouble(v);
+
+ // assert
+ Assertions.assertTrue(Double.isFinite(actual), () ->
+ "Computed norm was not finite; vector= " + Arrays.toString(v) + ", exact= " + exact +
+ ", direct= " + direct + ", actual= " + actual);
+
+ final int ulpError = Math.abs(DoubleTestUtils.computeUlpDifference(exact, actual));
+
+ Assertions.assertTrue(ulpError <= MAX_ULP_ERR, () ->
+ "Computed norm ulp error exceeds bounds; vector= " + Arrays.toString(v) +
+ ", exact= " + exact + ", actual= " + actual + ", ulpError= " + ulpError);
+ }
+ }
+
+ /** Assert that {@code directNorm(v) * 2^scaleExp = fn(v * 2^scaleExp)}.
+ * @param v unscaled vector
+ * @param scaleExp scale factor exponent
+ * @param fn euclidean norm function
+ */
+ private static void checkScaledEuclideanNorm(final double[] v, final int scaleExp,
+ final ToDoubleFunction<double[]> fn) {
+
+ final double scale = Math.scalb(1d, scaleExp);
+ final double[] scaledV = new double[v.length];
+ for (int i = 0; i < v.length; ++i) {
+ scaledV[i] = v[i] * scale;
+ }
+
+ final double norm = directEuclideanNorm(v);
+ final double scaledNorm = fn.applyAsDouble(scaledV);
+
+ Assertions.assertEquals(norm * scale, scaledNorm);
+ }
+
+ /** Direct euclidean norm computation.
+ * @param v array
+ * @return euclidean norm using direct summation.
+ */
+ private static double directEuclideanNorm(final double[] v) {
+ double n = 0;
+ for (int i = 0; i < v.length; i++) {
+ n += v[i] * v[i];
+ }
+ return Math.sqrt(n);
+ }
+
+ /** Compute the exact double value of the vector norm using BigDecimals
+ * with a math context of {@link MathContext#DECIMAL128}.
+ * @param v array
+ * @return euclidean norm using BigDecimal with MathContext.DECIMAL128
+ */
+ private static double exactEuclideanNorm(final double[] v) {
+ final MathContext ctx = MathContext.DECIMAL128;
+
+ BigDecimal sum = BigDecimal.ZERO;
+ for (final double d : v) {
+ sum = sum.add(BigDecimal.valueOf(d).pow(2), ctx);
+ }
+
+ return sum.sqrt(ctx).doubleValue();
+ }
+
+ /** Compute statistics for the ulp error of {@code fn} for the given inputs and
+ * array of exact results.
+ * @param inputs sample inputs
+ * @param exactResults array containing the exact expected results
+ * @param fn function to perform the computation
+ * @return ulp error statistics
+ */
+ private static UlpErrorStats computeUlpErrorStats(final double[][] inputs, final double[] exactResults,
+ final ToDoubleFunction<double[]> fn) {
+
+ // compute the ulp errors for each input
+ final int[] ulpErrors = new int[inputs.length];
+ int sum = 0;
+ for (int i = 0; i < inputs.length; ++i) {
+ final double exact = exactResults[i];
+ final double actual = fn.applyAsDouble(inputs[i]);
+
+ final int error = DoubleTestUtils.computeUlpDifference(exact, actual);
+ ulpErrors[i] = error;
+ sum += error;
+ }
+
+ // compute the mean
+ final double mean = sum / (double) ulpErrors.length;
+
+ // compute the std dev
+ double diffSumSq = 0d;
+ double diff;
+ for (int ulpError : ulpErrors) {
+ diff = ulpError - mean;
+ diffSumSq += diff * diff;
+ }
+
+ final double stdDev = Math.sqrt(diffSumSq / (inputs.length - 1));
+
+ return new UlpErrorStats(mean, stdDev);
+ }
+
+ /** Class containing ULP error statistics. */
+ private static final class UlpErrorStats {
+
+ private final double mean;
+
+ private final double stdDev;
+
+ UlpErrorStats(final double mean, final double stdDev) {
+ this.mean = mean;
+ this.stdDev = stdDev;
+ }
+
+ public double getMean() {
+ return mean;
+ }
+
+ public double getStdDev() {
+ return stdDev;
+ }
+ }
+}
diff --git a/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/SafeNormTest.java b/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/SafeNormTest.java
deleted file mode 100644
index de81135..0000000
--- a/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/SafeNormTest.java
+++ /dev/null
@@ -1,107 +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.arrays;
-
-import org.junit.jupiter.api.Assertions;
-import org.junit.jupiter.api.Test;
-
-/**
- * Test cases for the {@link SafeNorm} class.
- */
-class SafeNormTest {
-
- @Test
- void testTiny() {
- final double s = 1e-320;
- final double[] v = new double[] {s, s};
- Assertions.assertEquals(Math.sqrt(2) * s, SafeNorm.value(v));
- }
-
- @Test
- void testBig() {
- final double s = 1e300;
- final double[] v = new double[] {s, s};
- Assertions.assertEquals(Math.sqrt(2) * s, SafeNorm.value(v));
- }
-
- @Test
- void testOne3D() {
- final double s = 1;
- final double[] v = new double[] {s, s, s};
- Assertions.assertEquals(Math.sqrt(3), SafeNorm.value(v));
- }
-
- @Test
- void testUnit3D() {
- Assertions.assertEquals(1, SafeNorm.value(new double[] {1, 0, 0}));
- Assertions.assertEquals(1, SafeNorm.value(new double[] {0, 1, 0}));
- Assertions.assertEquals(1, SafeNorm.value(new double[] {0, 0, 1}));
- }
-
- @Test
- void testSimple() {
- final double[] v = new double[] {-0.9, 8.7, -6.5, -4.3, -2.1, 0, 1.2, 3.4, -5.6, 7.8, 9.0};
- final double expected = direct(v);
- Assertions.assertEquals(expected, SafeNorm.value(v));
- }
-
- @Test
- void testZero() {
- final double[] v = new double[] {0, 0, 0, 0, 0};
- Assertions.assertEquals(0d, SafeNorm.value(v));
- }
-
- @Test
- void testNaN() {
- final double[] v = new double[] {0, Double.NaN, 0};
- Assertions.assertEquals(Double.NaN, SafeNorm.value(v));
- }
-
- @Test
- void testInf() {
- final double[] v = new double[] {0, Double.NEGATIVE_INFINITY, 0};
- Assertions.assertEquals(Double.POSITIVE_INFINITY, SafeNorm.value(v));
- }
-
- @Test
- void testOverflow() {
- final double[] v = new double[] {0, Double.MAX_VALUE, Double.MAX_VALUE, 0};
- Assertions.assertEquals(Double.POSITIVE_INFINITY, SafeNorm.value(v));
- }
-
- @Test
- void testTinyAndSmallNormal() {
- // Ensure the sum of the squared values for 'normal' values (1e-19*1e-19)
- // is less than largest tiny value (1e-20)
- final double[] v = new double[] {1e-20, 1e-19};
- Assertions.assertEquals(Math.sqrt(101) * 1e-20, SafeNorm.value(v));
- }
-
- /**
- * Direct computation.
- *
- * @param v Array.
- * @return the norm using direct summation.
- */
- private double direct(double[] v) {
- double n = 0;
- for (int i = 0; i < v.length; i++) {
- n += v[i] * v[i];
- }
- return Math.sqrt(n);
- }
-}
diff --git a/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/SummationTest.java b/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/SummationTest.java
new file mode 100644
index 0000000..07d3128
--- /dev/null
+++ b/commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/SummationTest.java
@@ -0,0 +1,283 @@
+/*
+ * 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.arrays;
+
+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(BigDecimal.valueOf(value), MathContext.UNLIMITED);
+ }
+
+ return sum.doubleValue();
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/pom.xml b/commons-numbers-examples/examples-jmh/pom.xml
index cb14b49..40b10a8 100644
--- a/commons-numbers-examples/examples-jmh/pom.xml
+++ b/commons-numbers-examples/examples-jmh/pom.xml
@@ -80,6 +80,10 @@
</dependencies>
<properties>
+ <!-- JDK 9+ required for BigDecimal.sqrt() method. -->
+ <maven.compiler.source>1.9</maven.compiler.source>
+ <maven.compiler.target>1.9</maven.compiler.target>
+
<!-- OSGi -->
<commons.osgi.symbolicName>org.apache.commons.numbers.examples.jmh</commons.osgi.symbolicName>
<commons.osgi.export>org.apache.commons.numbers.examples.jmh</commons.osgi.export>
@@ -96,8 +100,43 @@
<project.mainClass>org.openjdk.jmh.Main</project.mainClass>
<!-- Disable analysis for benchmarking code. -->
<pmd.skip>true</pmd.skip>
+ <!-- Disable JDK compatibility check for benchmarking code. Also required since no signature
+ projects exist for JDK 9+. -->
+ <animal.sniffer.skip>true</animal.sniffer.skip>
</properties>
+ <build>
+ <plugins>
+ <plugin>
+ <!-- NOTE: javadoc config must also be set under <reporting> -->
+ <groupId>org.apache.maven.plugins</groupId>
+ <artifactId>maven-javadoc-plugin</artifactId>
+ <configuration>
+ <!-- Set source to JDK 8 to prevent issues with lack of module information. -->
+ <source>1.8</source>
+ <!-- Enable MathJax -->
+ <additionalOptions>${doclint.javadoc.qualifier} ${allowscript.javadoc.qualifier} -header '<script type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/${numbers.mathjax.version}/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>'</additionalOptions>
+ </configuration>
+ </plugin>
+ </plugins>
+ </build>
+
+ <reporting>
+ <plugins>
+ <plugin>
+ <!-- NOTE: javadoc config must also be set under <build> -->
+ <groupId>org.apache.maven.plugins</groupId>
+ <artifactId>maven-javadoc-plugin</artifactId>
+ <configuration>
+ <!-- Set source to JDK 8 to prevent issues with lack of module information. -->
+ <source>1.8</source>
+ <!-- Enable MathJax -->
+ <additionalOptions>${doclint.javadoc.qualifier} ${allowscript.javadoc.qualifier} -header '<script type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/${numbers.mathjax.version}/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>'</additionalOptions>
+ </configuration>
+ </plugin>
+ </plugins>
+ </reporting>
+
<profiles>
<profile>
<!-- Run a named benchmark from maven. The class to run can be specified as a property
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/DoubleUtils.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/DoubleUtils.java
new file mode 100644
index 0000000..9224898
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/DoubleUtils.java
@@ -0,0 +1,72 @@
+/*
+ * 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.arrays;
+
+import org.apache.commons.rng.UniformRandomProvider;
+
+/** Class containing utility methods for working with doubles.
+ */
+final class DoubleUtils {
+
+ /** No instantiation. */
+ private DoubleUtils() {}
+
+ /** Create a random double value with exponent in the range {@code [minExp, maxExp]}.
+ * @param minExp minimum exponent; must be less than {@code maxExp}
+ * @param maxExp maximum exponent; must be greater than {@code minExp}
+ * @param rng random number generator
+ * @return random double
+ */
+ static double random(final int minExp, final int maxExp, final UniformRandomProvider rng) {
+ // Create random doubles using random bits in the sign bit and the mantissa.
+ final long mask = ((1L << 52) - 1) | 1L << 63;
+ final long bits = rng.nextLong() & mask;
+ // The exponent must be unsigned so + 1023 to the signed exponent
+ final long exp = rng.nextInt(maxExp - minExp + 1) + minExp + 1023;
+ return Double.longBitsToDouble(bits | (exp << 52));
+ }
+
+ /** Create an array with the given length containing random doubles with exponents in the range
+ * {@code [minExp, maxExp]}.
+ * @param len array length
+ * @param minExp minimum exponent; must be less than {@code maxExp}
+ * @param maxExp maximum exponent; must be greater than {@code minExp}
+ * @param rng random number generator
+ * @return array of random doubles
+ */
+ static double[] randomArray(final int len, final int minExp, final int maxExp,
+ final UniformRandomProvider rng) {
+ final double[] arr = new double[len];
+ for (int i = 0; i < arr.length; ++i) {
+ arr[i] = random(minExp, maxExp, rng);
+ }
+ return arr;
+ }
+
+ /** Return a new array containing each element of {@code v} multiplied by {@code s}.
+ * @param v input array
+ * @param s scale factor
+ * @return new array containing each element of {@code v} multiplied by {@code s}
+ */
+ static double[] scalarMultiply(final double[] v, final double s) {
+ final double[] sv = new double[v.length];
+ for (int i = 0; i < v.length; ++i) {
+ sv[i] = s * v[i];
+ }
+ return sv;
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAlgorithmPerformance.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAlgorithmPerformance.java
new file mode 100644
index 0000000..43fb5e5
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAlgorithmPerformance.java
@@ -0,0 +1,225 @@
+/*
+ * 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.arrays;
+
+import java.util.concurrent.TimeUnit;
+import java.util.function.ToDoubleFunction;
+
+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;
+
+/**
+ * Execute benchmarks for the algorithms in the {@link EuclideanNormAlgorithms} 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 EuclideanNormAlgorithmPerformance {
+
+ /** String indicating double exponents with very low negative values, likely to underflow. */
+ private static final String LOW = "low";
+
+ /** String indicating double exponents with mid-level values which will not overflow or underflow. */
+ private static final String MID = "mid";
+
+ /** String indicating double exponents with very high positive values, likely to overflow. */
+ private static final String HIGH = "high";
+
+ /** String indicating double exponents over a very wide range of values. */
+ private static final String FULL = "full";
+
+ /** Class providing input vectors for benchmarks.
+ */
+ @State(Scope.Benchmark)
+ public static class VectorArrayInput {
+
+ /** The number of samples. */
+ @Param("100000")
+ private int samples;
+
+ /** The length of each vector. */
+ @Param("100")
+ private int vectorLength;
+
+ /** The type of double values placed in the vector arrays. */
+ @Param({LOW, MID, HIGH, FULL})
+ private String type;
+
+ /** Array of input vectors. */
+ private double[][] vectors;
+
+ /** Get the input vectors.
+ * @return input vectors
+ */
+ public double[][] getVectors() {
+ return vectors;
+ }
+
+ /** Create the input vectors for the instance.
+ */
+ @Setup
+ public void createVectors() {
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP);
+
+ int minExp;
+ int maxExp;
+
+ switch (type) {
+ case LOW:
+ minExp = -530;
+ maxExp = -510;
+ break;
+ case MID:
+ minExp = -10;
+ maxExp = +10;
+ break;
+ case HIGH:
+ minExp = +510;
+ maxExp = +530;
+ break;
+ default:
+ throw new IllegalArgumentException("Invalid vector type: " + type);
+ }
+
+ vectors = new double[samples][];
+ for (int i = 0; i < vectors.length; ++i) {
+ vectors[i] = DoubleUtils.randomArray(vectorLength, minExp, maxExp, rng);
+ }
+ }
+ }
+
+ /** Evaluate a norm computation method with the given input.
+ * @param fn function to evaluate
+ * @param input computation input
+ * @param bh blackhole
+ */
+ private static void eval(final ToDoubleFunction<double[]> fn, final VectorArrayInput input,
+ final Blackhole bh) {
+ final double[][] vectors = input.getVectors();
+ for (int i = 0; i < vectors.length; ++i) {
+ bh.consume(fn.applyAsDouble(vectors[i]));
+ }
+ }
+
+ /** Compute the performance of the {@link EuclideanNormAlgorithms.Exact} class.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void exact(final VectorArrayInput input, final Blackhole bh) {
+ eval(new EuclideanNormAlgorithms.Exact(), input, bh);
+ }
+
+ /** Compute the performance of the {@link EuclideanNormAlgorithms.Direct} class.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void direct(final VectorArrayInput input, final Blackhole bh) {
+ eval(new EuclideanNormAlgorithms.Direct(), input, bh);
+ }
+
+ /** Compute the performance of the {@link EuclideanNormAlgorithms.Enorm} class.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void enorm(final VectorArrayInput input, final Blackhole bh) {
+ eval(new EuclideanNormAlgorithms.Enorm(), input, bh);
+ }
+
+ /** Compute the performance of the {@link EuclideanNormAlgorithms.EnormMod} class.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void enormMod(final VectorArrayInput input, final Blackhole bh) {
+ eval(new EuclideanNormAlgorithms.EnormMod(), input, bh);
+ }
+
+ /** Compute the performance of the {@link EuclideanNormAlgorithms.EnormModKahan} class.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void enormModKahan(final VectorArrayInput input, final Blackhole bh) {
+ eval(new EuclideanNormAlgorithms.EnormModKahan(), input, bh);
+ }
+
+ /** Compute the performance of the {@link EuclideanNormAlgorithms.EnormModExt} class.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void enormModExt(final VectorArrayInput input, final Blackhole bh) {
+ eval(new EuclideanNormAlgorithms.EnormModExt(), input, bh);
+ }
+
+ /** Compute the performance of the {@link EuclideanNormAlgorithms.ExtendedPrecisionLinearCombination} class.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void extLinear(final VectorArrayInput input, final Blackhole bh) {
+ eval(new EuclideanNormAlgorithms.ExtendedPrecisionLinearCombination(), input, bh);
+ }
+
+
+ /** Compute the performance of the {@link EuclideanNormAlgorithms.ExtendedPrecisionLinearCombinationMod} class.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void extLinearMod(final VectorArrayInput input, final Blackhole bh) {
+ eval(new EuclideanNormAlgorithms.ExtendedPrecisionLinearCombinationMod(), input, bh);
+ }
+
+ /** Compute the performance of the {@link EuclideanNormAlgorithms.ExtendedPrecisionLinearCombinationSinglePass}
+ * class.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void extLinearSinglePass(final VectorArrayInput input, final Blackhole bh) {
+ eval(new EuclideanNormAlgorithms.ExtendedPrecisionLinearCombinationSinglePass(), input, bh);
+ }
+
+ /** Compute the performance of the {@link EuclideanNormAlgorithms.ExtendedPrecisionLinearCombinationSqrt2}
+ * class.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void extLinearSqrt2(final VectorArrayInput input, final Blackhole bh) {
+ eval(new EuclideanNormAlgorithms.ExtendedPrecisionLinearCombinationSqrt2(), input, bh);
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAlgorithms.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAlgorithms.java
new file mode 100644
index 0000000..69008b7
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAlgorithms.java
@@ -0,0 +1,749 @@
+/*
+ * 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.arrays;
+
+import java.math.BigDecimal;
+import java.math.MathContext;
+import java.util.function.ToDoubleFunction;
+
+/** Class containing various Euclidean norm computation methods for comparison.
+ */
+public final class EuclideanNormAlgorithms {
+
+ /** No instantiation. */
+ private EuclideanNormAlgorithms() {}
+
+ /** Exact computation method using {@link BigDecimal} and {@link MathContext#DECIMAL128}.
+ */
+ static final class Exact implements ToDoubleFunction<double[]> {
+
+ /** {@inheritDoc} */
+ @Override
+ public double applyAsDouble(final double[] v) {
+ // compute the sum of squares
+ final MathContext ctx = MathContext.DECIMAL128;
+
+ BigDecimal sum = BigDecimal.ZERO;
+ BigDecimal n;
+ for (int i = 0; i < v.length; ++i) {
+ n = BigDecimal.valueOf(Double.isFinite(v[i]) ? v[i] : 0d);
+ sum = sum.add(n.multiply(n, ctx), ctx);
+ }
+
+ return sum.sqrt(ctx).doubleValue();
+ }
+ }
+
+ /** Direct computation method that simply computes the sums of squares and takes
+ * the square root with no special handling of values.
+ */
+ static final class Direct implements ToDoubleFunction<double[]> {
+
+ /** {@inheritDoc} */
+ @Override
+ public double applyAsDouble(final double[] v) {
+ double n = 0;
+ for (int i = 0; i < v.length; i++) {
+ n += v[i] * v[i];
+ }
+ return Math.sqrt(n);
+ }
+ }
+
+ /** Translation of the <a href="http://www.netlib.org/minpack">minpack</a>
+ * "enorm" subroutine. This method handles overflow and underflow.
+ */
+ static final class Enorm implements ToDoubleFunction<double[]> {
+
+ /** Constant. */
+ private static final double R_DWARF = 3.834e-20;
+ /** Constant. */
+ private static final double R_GIANT = 1.304e+19;
+
+ /** {@inheritDoc} */
+ @Override
+ public double applyAsDouble(final double[] v) {
+ double s1 = 0;
+ double s2 = 0;
+ double s3 = 0;
+ double x1max = 0;
+ double x3max = 0;
+ final double floatn = v.length;
+ final double agiant = R_GIANT / floatn;
+ for (int i = 0; i < v.length; i++) {
+ final double xabs = Math.abs(v[i]);
+ if (xabs < R_DWARF || xabs > agiant) {
+ if (xabs > R_DWARF) {
+ if (xabs > x1max) {
+ final double r = x1max / xabs;
+ s1 = 1 + s1 * r * r;
+ x1max = xabs;
+ } else {
+ final double r = xabs / x1max;
+ s1 += r * r;
+ }
+ } else {
+ if (xabs > x3max) {
+ final double r = x3max / xabs;
+ s3 = 1 + s3 * r * r;
+ x3max = xabs;
+ } else {
+ if (xabs != 0) {
+ final double r = xabs / x3max;
+ s3 += r * r;
+ }
+ }
+ }
+ } else {
+ s2 += xabs * xabs;
+ }
+ }
+ double norm;
+ if (s1 != 0) {
+ norm = x1max * Math.sqrt(s1 + (s2 / x1max) / x1max);
+ } else {
+ if (s2 == 0) {
+ norm = x3max * Math.sqrt(s3);
+ } else {
+ if (s2 >= x3max) {
+ norm = Math.sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));
+ } else {
+ norm = Math.sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
+ }
+ }
+ }
+ return norm;
+ }
+ }
+
+ /** Modified version of {@link Enorm} created by Alex Herbert.
+ */
+ static final class EnormMod implements ToDoubleFunction<double[]> {
+
+ /** {@inheritDoc} */
+ @Override
+ public double applyAsDouble(final double[] v) {
+ // Sum of big, normal and small numbers with 2-fold extended precision summation
+ double s1 = 0;
+ double s2 = 0;
+ double s3 = 0;
+ for (int i = 0; i < v.length; i++) {
+ final double x = Math.abs(v[i]);
+ if (!(x <= Double.MAX_VALUE)) {
+ return x;
+ } else if (x > 0x1.0p500) {
+ // Scale down big numbers
+ s1 += square(x * 0x1.0p-600);
+ } else if (x < 0x1.0p-500) {
+ // Scale up small numbers
+ s3 += square(x * 0x1.0p600);
+ } else {
+ // Unscaled
+ s2 += square(x);
+ }
+ }
+ // The highest sum is the significant component. Add the next significant.
+ if (s1 != 0) {
+ return Math.sqrt(s1 + s2 * 0x1.0p-600 * 0x1.0p-600) * 0x1.0p600;
+ } else if (s2 != 0) {
+ return Math.sqrt(s2 + s3 * 0x1.0p-600 * 0x1.0p-600);
+ }
+ return Math.sqrt(s3) * 0x1.0p-600;
+ }
+
+ /** Compute the square of {@code x}.
+ * @param x input value
+ * @return square of {@code x}
+ */
+ private static double square(final double x) {
+ return x * x;
+ }
+ }
+
+ /** Version of {@link EnormMod} using Kahan summation.
+ */
+ static final class EnormModKahan implements ToDoubleFunction<double[]> {
+
+ /** Threshold for scaling small numbers. */
+ private static final double SMALL_THRESH = 0x1.0p-500;
+
+ /** Threshold for scaling large numbers. */
+ private static final double LARGE_THRESH = 0x1.0p+500;
+
+ /** Value used to scale down large numbers. */
+ private static final double SCALE_DOWN = 0x1.0p-600;
+
+ /** Value used to scale up small numbers. */
+ private static final double SCALE_UP = 0x1.0p+600;
+
+ /** {@inheritDoc} */
+ @Override
+ public double applyAsDouble(final double[] v) {
+ // Sum of big, normal and small numbers
+ double s1 = 0;
+ double s2 = 0;
+ double s3 = 0;
+ double c1 = 0;
+ double c2 = 0;
+ double c3 = 0;
+ for (int i = 0; i < v.length; i++) {
+ final double x = Math.abs(v[i]);
+ if (x > LARGE_THRESH) {
+ // Scale down big numbers
+ final double y = square(x * SCALE_DOWN) - c1;
+ final double t = s1 + y;
+ c1 = (t - s1) - y;
+ s1 = t;
+
+ } else if (x < SMALL_THRESH) {
+ // Scale up small numbers
+ final double y = square(x * SCALE_UP) - c3;
+ final double t = s3 + y;
+ c3 = (t - s3) - y;
+ s3 = t;
+ } else {
+ // Unscaled
+ final double y = square(x) - c2;
+ final double t = s2 + y;
+ c2 = (t - s2) - y;
+ s2 = t;
+ }
+ }
+ // The highest sum is the significant component. Add the next significant.
+ // Add the scaled compensation then the scaled sum.
+ if (s1 != 0) {
+ double y = c2 * SCALE_DOWN * SCALE_DOWN - c1;
+ final double t = s1 + y;
+ c1 = (t - s1) - y;
+ y = s2 * SCALE_DOWN * SCALE_DOWN - c1;
+ return Math.sqrt(t + y) * SCALE_UP;
+ } else if (s2 != 0) {
+ double y = c3 * SCALE_DOWN * SCALE_DOWN - c2;
+ final double t = s2 + y;
+ c2 = (t - s2) - y;
+ y = s3 * SCALE_DOWN * SCALE_DOWN - c2;
+ return Math.sqrt(t + y);
+ }
+ return Math.sqrt(s3) * SCALE_DOWN;
+ }
+
+ /** Compute the square of {@code x}.
+ * @param x input value
+ * @return square of {@code x}
+ */
+ private static double square(final double x) {
+ return x * x;
+ }
+ }
+
+ /** Version of {@link EnormMod} using extended precision summation.
+ */
+ static final class EnormModExt implements ToDoubleFunction<double[]> {
+
+ /** Threshold for scaling small numbers. */
+ private static final double SMALL_THRESH = 0x1.0p-500;
+
+ /** Threshold for scaling large numbers. */
+ private static final double LARGE_THRESH = 0x1.0p+500;
+
+ /** Value used to scale down large numbers. */
+ private static final double SCALE_DOWN = 0x1.0p-600;
+
+ /** Value used to scale up small numbers. */
+ private static final double SCALE_UP = 0x1.0p+600;
+
+ /** {@inheritDoc} */
+ @Override
+ public double applyAsDouble(final double[] v) {
+ // Sum of big, normal and small numbers with 2-fold extended precision summation
+ double s1 = 0;
+ double s2 = 0;
+ double s3 = 0;
+ double c1 = 0;
+ double c2 = 0;
+ double c3 = 0;
+ for (int i = 0; i < v.length; i++) {
+ final double x = Math.abs(v[i]);
+ if (!(x <= Double.MAX_VALUE)) {
+ return x;
+ } else if (x > LARGE_THRESH) {
+ // Scale down big numbers
+ final double y = square(x * SCALE_DOWN);
+ final double t = s1 + y;
+ c1 += DoublePrecision.twoSumLow(s1, y, t);
+ s1 = t;
+ } else if (x < SMALL_THRESH) {
+ // Scale up small numbers
+ final double y = square(x * SCALE_UP);
+ final double t = s3 + y;
+ c3 += DoublePrecision.twoSumLow(s3, y, t);
+ s3 = t;
+ } else {
+ // Unscaled
+ final double y = square(x);
+ final double t = s2 + y;
+ c2 += DoublePrecision.twoSumLow(s2, y, t);
+ s2 = t;
+ }
+ }
+ // The highest sum is the significant component. Add the next significant.
+ // Adapted from LinearCombination dot2s summation.
+ if (s1 != 0) {
+ s2 = s2 * SCALE_DOWN * SCALE_DOWN;
+ c2 = c2 * SCALE_DOWN * SCALE_DOWN;
+ double sum = s1 + s2;
+ c1 += DoublePrecision.twoSumLow(s1, s2, sum) + c2;
+ return Math.sqrt(sum + c1) * SCALE_UP;
+ } else if (s2 != 0) {
+ s3 = s3 * SCALE_DOWN * SCALE_DOWN;
+ c3 = c3 * SCALE_DOWN * SCALE_DOWN;
+ double sum = s2 + s3;
+ c2 += DoublePrecision.twoSumLow(s2, s3, sum) + c3;
+ return Math.sqrt(sum + c2);
+ }
+ return Math.sqrt(s3) * SCALE_DOWN;
+ }
+
+ /** Compute the square of {@code x}.
+ * @param x input value
+ * @return square of {@code x}
+ */
+ private static double square(final double x) {
+ return x * x;
+ }
+ }
+
+ /** Euclidean norm computation algorithm that uses {@link LinearCombinations} to perform
+ * an extended precision summation.
+ */
+ static final class ExtendedPrecisionLinearCombination implements ToDoubleFunction<double[]> {
+
+ /** {@inheritDoc} */
+ @Override
+ public double applyAsDouble(final double[] v) {
+ // Find the magnitude limits ignoring zero
+ double max = 0;
+ double min = Double.POSITIVE_INFINITY;
+ for (int i = 0; i < v.length; i++) {
+ final double x = Math.abs(v[i]);
+ if (Double.isNaN(x)) {
+ return x;
+ } else if (x > max) {
+ max = x;
+ } else if (x < min && x != 0) {
+ min = x;
+ }
+ }
+ // Edge cases
+ if (max == 0 || max == Double.POSITIVE_INFINITY) {
+ return max;
+ }
+ // Use scaling if required
+ double[] x = v;
+ double rescale = 1;
+ if (max > 0x1.0p500) {
+ // Too big so scale down
+ x = x.clone();
+ for (int i = 0; i < x.length; i++) {
+ x[i] *= 0x1.0p-600;
+ }
+ rescale = 0x1.0p600;
+ } else if (min < 0x1.0p-500 && max < 0x1.0p100) {
+ // Too small so scale up
+ x = x.clone();
+ for (int i = 0; i < x.length; i++) {
+ x[i] *= 0x1.0p600;
+ }
+ rescale = 0x1.0p-600;
+ }
+ return Math.sqrt(org.apache.commons.numbers.arrays.LinearCombination.value(x, x)) * rescale;
+ }
+ }
+
+ /** Modification of {@link ExtendedPrecisionLinearCombination} that uses an optimized version of the
+ * linear combination computation.
+ */
+ static final class ExtendedPrecisionLinearCombinationMod implements ToDoubleFunction<double[]> {
+
+ /** Threshold for scaling small numbers. */
+ private static final double SMALL_THRESH = 0x1.0p-500;
+
+ /** Threshold for scaling large numbers. */
+ private static final double LARGE_THRESH = 0x1.0p+500;
+
+ /** Value used to scale down large numbers. */
+ private static final double SCALE_DOWN = 0x1.0p-600;
+
+ /** Value used to scale up small numbers. */
+ private static final double SCALE_UP = 0x1.0p+600;
+
+ /** {@inheritDoc} */
+ @Override
+ public double applyAsDouble(final double[] v) {
+ // Find the magnitude limits ignoring zero
+ double max = 0;
+ double min = Double.POSITIVE_INFINITY;
+ for (int i = 0; i < v.length; i++) {
+ final double x = Math.abs(v[i]);
+ if (Double.isNaN(x)) {
+ return x;
+ } else if (x > max) {
+ max = x;
+ } else if (x < min && x != 0) {
+ min = x;
+ }
+ }
+ // Edge cases
+ if (max == 0 || max == Double.POSITIVE_INFINITY) {
+ return max;
+ }
+
+ // Below here no value is infinite or NaN.
+
+ // Use scaling if required
+ double scale = 1;
+ double rescale = 1;
+ if (max > LARGE_THRESH) {
+ // Too big so scale down
+ scale = SCALE_DOWN;
+ rescale = SCALE_UP;
+ } else if (min < SMALL_THRESH && max < 0x1.0p100) {
+ // Too small so scale up
+ scale = SCALE_UP;
+ rescale = SCALE_DOWN;
+ }
+
+ // Same as LinearCombination but with scaling.
+ // Splitting is safe due to scaling and only one term requires splitting.
+
+ // Implement dot2s (Algorithm 5.4) from Ogita et al (2005).
+ final int len = v.length;
+
+ // p is the standard scalar product sum.
+ // s is the sum of round-off parts.
+ double a = v[0] * scale;
+ double p = a * a;
+ double s = productLowUnscaled(a, p);
+
+ // Remaining split products added to the current sum and round-off sum.
+ for (int i = 1; i < len; i++) {
+ a = v[i] * scale;
+ final double h = a * a;
+ final double r = productLowUnscaled(a, h);
+
+ final double x = p + h;
+ // s_i = s_(i-1) + (q_i + r_i)
+ s += DoublePrecision.twoSumLow(p, h, x) + r;
+ p = x;
+ }
+ p += s;
+
+ return Math.sqrt(p) * rescale;
+ }
+
+ /**
+ * Compute the low part of the double length number {@code (z,zz)} for the exact
+ * square of {@code x} using Dekker's mult12 algorithm. The standard
+ * precision product {@code x*x} must be provided. The number {@code x}
+ * is split into high and low parts using Dekker's algorithm.
+ *
+ * <p>Warning: This method does not perform scaling in Dekker's split and large
+ * finite numbers can create NaN results.
+ *
+ * @param x The factor.
+ * @param xx Square of the factor (x * x).
+ * @return the low part of the product double length number
+ */
+ private static double productLowUnscaled(double x, double xx) {
+ // Split the numbers using Dekker's algorithm without scaling
+ final double hx = DoublePrecision.highPartUnscaled(x);
+ final double lx = x - hx;
+
+ return DoublePrecision.productLow(hx, lx, hx, lx, xx);
+ }
+ }
+
+ /** Modification of {@link ExtendedPrecisionLinearCombination} that only uses a single pass through
+ * the input array.
+ */
+ static final class ExtendedPrecisionLinearCombinationSinglePass implements ToDoubleFunction<double[]> {
+
+ /** Threshold for scaling small numbers. */
+ private static final double SMALL_THRESH = 0x1.0p-500;
+
+ /** Threshold for scaling large numbers. */
+ private static final double LARGE_THRESH = 0x1.0p+500;
+
+ /** Value used to scale down large numbers. */
+ private static final double SCALE_DOWN = 0x1.0p-600;
+
+ /** Value used to scale up small numbers. */
+ private static final double SCALE_UP = 0x1.0p+600;
+
+ /** {@inheritDoc} */
+ @Override
+ public double applyAsDouble(final double[] v) {
+ double s1 = 0;
+ double s2 = 0;
+ double s3 = 0;
+ double c1 = 0;
+ double c2 = 0;
+ double c3 = 0;
+ for (int i = 0; i < v.length; ++i) {
+ final double x = Math.abs(v[i]);
+ if (Double.isNaN(x)) {
+ // found an NaN; no use in continuing
+ return x;
+ } else if (x > LARGE_THRESH) {
+ // scale down
+ final double sx = x * SCALE_DOWN;
+
+ // compute the product and product correction
+ final double p = sx * sx;
+ final double cp = productLowUnscaled(sx, p);
+
+ // compute the running sum and sum correction
+ final double s = s1 + p;
+ final double cs = DoublePrecision.twoSumLow(s1, p, s);
+
+ // update running totals
+ c1 += cp + cs;
+ s1 = s;
+ } else if (x < SMALL_THRESH) {
+ // scale up
+ final double sx = x * SCALE_UP;
+
+ // compute the product and product correction
+ final double p = sx * sx;
+ final double cp = productLowUnscaled(sx, p);
+
+ // compute the running sum and sum correction
+ final double s = s3 + p;
+ final double cs = DoublePrecision.twoSumLow(s3, p, s);
+
+ // update running totals
+ c3 += cp + cs;
+ s3 = s;
+ } else {
+ // no scaling
+ // compute the product and product correction
+ final double p = x * x;
+ final double cp = productLowUnscaled(x, p);
+
+ // compute the running sum and sum correction
+ final double s = s2 + p;
+ final double cs = DoublePrecision.twoSumLow(s2, p, s);
+
+ // update running totals
+ c2 += cp + cs;
+ s2 = s;
+ }
+ }
+
+ if (s1 != 0) {
+ // add s1, s2, c1, c2
+ s2 = s2 * SCALE_DOWN * SCALE_DOWN;
+ c2 = c2 * SCALE_DOWN * SCALE_DOWN;
+ final double sum = s1 + s2;
+ c1 += DoublePrecision.twoSumLow(s1, s2, sum) + c2;
+ return Math.sqrt(sum + c1) * SCALE_UP;
+ } else if (s2 != 0) {
+ // add s2, s3, c2, c3
+ s3 = s3 * SCALE_DOWN * SCALE_DOWN;
+ c3 = c3 * SCALE_DOWN * SCALE_DOWN;
+ final double sum = s2 + s3;
+ c2 += DoublePrecision.twoSumLow(s2, s3, sum) + c3;
+ return Math.sqrt(sum + c2);
+ }
+ // add s3, c3
+ return Math.sqrt(s3 + c3) * SCALE_DOWN;
+ }
+
+ /**
+ * Compute the low part of the double length number {@code (z,zz)} for the exact
+ * square of {@code x} using Dekker's mult12 algorithm. The standard
+ * precision product {@code x*x} must be provided. The number {@code x}
+ * is split into high and low parts using Dekker's algorithm.
+ *
+ * <p>Warning: This method does not perform scaling in Dekker's split and large
+ * finite numbers can create NaN results.
+ *
+ * @param x The factor.
+ * @param xx Square of the factor (x * x).
+ * @return the low part of the product double length number
+ */
+ private static double productLowUnscaled(double x, double xx) {
+ // Split the numbers using Dekker's algorithm without scaling
+ final double hx = DoublePrecision.highPartUnscaled(x);
+ final double lx = x - hx;
+
+ return DoublePrecision.productLow(hx, lx, hx, lx, xx);
+ }
+ }
+
+ /** Modification of {@link ExtendedPrecisionLinearCombination} that only uses a single pass through
+ * the input array as well as an extended precision sqrt computation.
+ */
+ static final class ExtendedPrecisionLinearCombinationSqrt2 implements ToDoubleFunction<double[]> {
+
+ /** Threshold for scaling small numbers. */
+ private static final double SMALL_THRESH = 0x1.0p-500;
+
+ /** Threshold for scaling large numbers. */
+ private static final double LARGE_THRESH = 0x1.0p+500;
+
+ /** Value used to scale down large numbers. */
+ private static final double SCALE_DOWN = 0x1.0p-600;
+
+ /** Value used to scale up small numbers. */
+ private static final double SCALE_UP = 0x1.0p+600;
+
+ /** {@inheritDoc} */
+ @Override
+ public double applyAsDouble(final double[] v) {
+ double s1 = 0;
+ double s2 = 0;
+ double s3 = 0;
+ double c1 = 0;
+ double c2 = 0;
+ double c3 = 0;
+ for (int i = 0; i < v.length; ++i) {
+ final double x = Math.abs(v[i]);
+ if (Double.isNaN(x)) {
+ // found an NaN; no use in continuing
+ return x;
+ } else if (x > LARGE_THRESH) {
+ // scale down
+ final double sx = x * SCALE_DOWN;
+
+ // compute the product and product correction
+ final double p = sx * sx;
+ final double cp = productLowUnscaled(sx, p);
+
+ // compute the running sum and sum correction
+ final double s = s1 + p;
+ final double cs = DoublePrecision.twoSumLow(s1, p, s);
+
+ // update running totals
+ c1 += cp + cs;
+ s1 = s;
+ } else if (x < SMALL_THRESH) {
+ // scale up
+ final double sx = x * SCALE_UP;
+
+ // compute the product and product correction
+ final double p = sx * sx;
+ final double cp = productLowUnscaled(sx, p);
+
+ // compute the running sum and sum correction
+ final double s = s3 + p;
+ final double cs = DoublePrecision.twoSumLow(s3, p, s);
+
+ // update running totals
+ c3 += cp + cs;
+ s3 = s;
+ } else {
+ // no scaling
+ // compute the product and product correction
+ final double p = x * x;
+ final double cp = productLowUnscaled(x, p);
+
+ // compute the running sum and sum correction
+ final double s = s2 + p;
+ final double cs = DoublePrecision.twoSumLow(s2, p, s);
+
+ // update running totals
+ c2 += cp + cs;
+ s2 = s;
+ }
+ }
+
+ if (s1 != 0) {
+ // add s1, s2, c1, c2
+ s2 = s2 * SCALE_DOWN * SCALE_DOWN;
+ c2 = c2 * SCALE_DOWN * SCALE_DOWN;
+ final double sum = s1 + s2;
+ c1 += DoublePrecision.twoSumLow(s1, s2, sum) + c2;
+
+ final double f = sum + c1;
+ final double ff = DoublePrecision.twoSumLow(sum, c1, f);
+ return sqrt2(f, ff) * SCALE_UP;
+ } else if (s2 != 0) {
+ // add s2, s3, c2, c3
+ s3 = s3 * SCALE_DOWN * SCALE_DOWN;
+ c3 = c3 * SCALE_DOWN * SCALE_DOWN;
+ final double sum = s2 + s3;
+ c2 += DoublePrecision.twoSumLow(s2, s3, sum) + c3;
+
+ final double f = sum + c2;
+ final double ff = DoublePrecision.twoSumLow(sum, c2, f);
+ return sqrt2(f, ff);
+ }
+ // add s3, c3
+ final double f = s3 + c3;
+ final double ff = DoublePrecision.twoSumLow(s3, c3, f);
+ return sqrt2(f, ff) * SCALE_DOWN;
+ }
+
+ /**
+ * Compute the low part of the double length number {@code (z,zz)} for the exact
+ * square of {@code x} using Dekker's mult12 algorithm. The standard
+ * precision product {@code x*x} must be provided. The number {@code x}
+ * is split into high and low parts using Dekker's algorithm.
+ *
+ * <p>Warning: This method does not perform scaling in Dekker's split and large
+ * finite numbers can create NaN results.
+ *
+ * @param x The factor.
+ * @param xx Square of the factor (x * x).
+ * @return the low part of the product double length number
+ */
+ private static double productLowUnscaled(double x, double xx) {
+ // Split the numbers using Dekker's algorithm without scaling
+ final double hx = DoublePrecision.highPartUnscaled(x);
+ final double lx = x - hx;
+
+ return DoublePrecision.productLow(hx, lx, hx, lx, xx);
+ }
+
+ /**
+ * Compute the extended precision square root from the split number
+ * {@code x, xx}.
+ * This is a modification of Dekker's sqrt2 algorithm to ignore the
+ * roundoff of the square root.
+ *
+ * @param x the high part
+ * @param xx the low part
+ * @return the double
+ */
+ private static double sqrt2(final double x, final double xx) {
+ if (x > 0) {
+ double c = Math.sqrt(x);
+ double u = c * c;
+ //double uu = ExtendedPrecision.productLow(c, c, u);
+ // Here we use the optimised version:
+ double uu = productLowUnscaled(c, u);
+ double cc = (x - u - uu + xx) * 0.5 / c;
+ // Extended precision sqrt (y, yy)
+ // y = c + cc
+ // yy = c - y + cc (ignored)
+ return c + cc;
+ }
+ return x;
+ }
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormEvaluator.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormEvaluator.java
new file mode 100644
index 0000000..71f5ca6
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormEvaluator.java
@@ -0,0 +1,249 @@
+/*
+ * 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.arrays;
+
+import java.math.BigDecimal;
+import java.math.MathContext;
+import java.util.HashMap;
+import java.util.LinkedHashMap;
+import java.util.Map;
+import java.util.function.ToDoubleFunction;
+
+/**
+ * Class used to evaluate the accuracy of different norm computation
+ * methods.
+ */
+public class EuclideanNormEvaluator {
+
+ /** Map of names to norm computation methods. */
+ private final Map<String, ToDoubleFunction<double[]>> methods = new LinkedHashMap<>();
+
+ /** Add a computation method to be evaluated.
+ * @param name method name
+ * @param method computation method
+ * @return this instance
+ */
+ public EuclideanNormEvaluator addMethod(final String name, final ToDoubleFunction<double[]> method) {
+ methods.put(name, method);
+ return this;
+ }
+
+ /** Evaluate the configured computation methods against the given array of input vectors.
+ * @param inputs array of input vectors
+ * @return map of evaluation results keyed by method name
+ */
+ public Map<String, Stats> evaluate(final double[][] inputs) {
+
+ final Map<String, StatsAccumulator> accumulators = new HashMap<>();
+ for (final String name : methods.keySet()) {
+ accumulators.put(name, new StatsAccumulator(inputs.length * 2));
+ }
+
+ for (int i = 0; i < inputs.length; ++i) {
+ // compute the norm in a forward and reverse directions to include
+ // summation artifacts
+ final double[] vec = inputs[i];
+
+ final double[] reverseVec = new double[vec.length];
+ for (int j = 0; j < vec.length; ++j) {
+ reverseVec[vec.length - 1 - j] = vec[j];
+ }
+
+ final double exact = computeExact(vec);
+
+ for (final Map.Entry<String, ToDoubleFunction<double[]>> entry : methods.entrySet()) {
+ final ToDoubleFunction<double[]> fn = entry.getValue();
+
+ final StatsAccumulator acc = accumulators.get(entry.getKey());
+
+ final double forwardSample = fn.applyAsDouble(vec);
+ acc.report(exact, forwardSample);
+
+ final double reverseSample = fn.applyAsDouble(reverseVec);
+ acc.report(exact, reverseSample);
+ }
+ }
+
+ final Map<String, Stats> stats = new LinkedHashMap<>();
+ for (final String name : methods.keySet()) {
+ stats.put(name, accumulators.get(name).computeStats());
+ }
+
+ return stats;
+ }
+
+ /** Compute the exact double value of the vector norm using BigDecimals
+ * with a math context of {@link MathContext#DECIMAL128}.
+ * @param vec input vector
+ * @return euclidean norm
+ */
+ private static double computeExact(final double[] vec) {
+ final MathContext ctx = MathContext.DECIMAL128;
+
+ BigDecimal sum = BigDecimal.ZERO;
+ for (final double v : vec) {
+ sum = sum.add(BigDecimal.valueOf(v).pow(2), ctx);
+ }
+
+ return sum.sqrt(ctx).doubleValue();
+ }
+
+ /** Compute the ulp difference between two values of the same sign.
+ * @param a first input
+ * @param b second input
+ * @return ulp difference between the arguments
+ */
+ private static int computeUlpDifference(final double a, final double b) {
+ return (int) (Double.doubleToLongBits(a) - Double.doubleToLongBits(b));
+ }
+
+ /** Class containing evaluation statistics for a single computation method.
+ */
+ public static final class Stats {
+
+ /** Mean ulp error. */
+ private final double ulpErrorMean;
+
+ /** Ulp error standard deviation. */
+ private final double ulpErrorStdDev;
+
+ /** Ulp error minimum value. */
+ private final double ulpErrorMin;
+
+ /** Ulp error maximum value. */
+ private final double ulpErrorMax;
+
+ /** Number of failed computations. */
+ private final int failCount;
+
+ /** Construct a new instance.
+ * @param ulpErrorMean ulp error mean
+ * @param ulpErrorStdDev ulp error standard deviation
+ * @param ulpErrorMin ulp error minimum value
+ * @param ulpErrorMax ulp error maximum value
+ * @param failCount number of failed computations
+ */
+ Stats(final double ulpErrorMean, final double ulpErrorStdDev, final double ulpErrorMin,
+ final double ulpErrorMax, final int failCount) {
+ this.ulpErrorMean = ulpErrorMean;
+ this.ulpErrorStdDev = ulpErrorStdDev;
+ this.ulpErrorMin = ulpErrorMin;
+ this.ulpErrorMax = ulpErrorMax;
+ this.failCount = failCount;
+ }
+
+ /** Get the ulp error mean.
+ * @return ulp error mean
+ */
+ public double getUlpErrorMean() {
+ return ulpErrorMean;
+ }
+
+ /** Get the ulp error standard deviation.
+ * @return ulp error standard deviation
+ */
+ public double getUlpErrorStdDev() {
+ return ulpErrorStdDev;
+ }
+
+ /** Get the ulp error minimum value.
+ * @return ulp error minimum value
+ */
+ public double getUlpErrorMin() {
+ return ulpErrorMin;
+ }
+
+ /** Get the ulp error maximum value.
+ * @return ulp error maximum value
+ */
+ public double getUlpErrorMax() {
+ return ulpErrorMax;
+ }
+
+ /** Get the number of failed computations, meaning the number of
+ * computations that overflowed or underflowed.
+ * @return number of failed computations
+ */
+ public int getFailCount() {
+ return failCount;
+ }
+ }
+
+ /** Class used to accumulate statistics during a norm evaluation run.
+ */
+ private static final class StatsAccumulator {
+
+ /** Sample index. */
+ private int sampleIdx;
+
+ /** Array of ulp errors for each sample. */
+ private final double[] ulpErrors;
+
+ /** Construct a new instance.
+ * @param count number of samples to be accumulated
+ */
+ StatsAccumulator(final int count) {
+ ulpErrors = new double[count];
+ }
+
+ /** Report a computation result.
+ * @param expected expected result
+ * @param actual actual result
+ */
+ public void report(final double expected, final double actual) {
+ ulpErrors[sampleIdx++] = Double.isFinite(actual) && actual != 0.0 ?
+ computeUlpDifference(expected, actual) :
+ Double.NaN;
+ }
+
+ /** Compute the final statistics for the run.
+ * @return statistics object
+ */
+ public Stats computeStats() {
+ int successCount = 0;
+ double sum = 0d;
+ double min = Double.POSITIVE_INFINITY;
+ double max = Double.NEGATIVE_INFINITY;
+
+ for (double ulpError : ulpErrors) {
+ if (Double.isFinite(ulpError)) {
+ ++successCount;
+ min = Math.min(ulpError, min);
+ max = Math.max(ulpError, max);
+ sum += ulpError;
+ }
+ }
+
+ final double mean = sum / successCount;
+
+ double diffSumSq = 0d;
+ double diff;
+ for (double ulpError : ulpErrors) {
+ if (Double.isFinite(ulpError)) {
+ diff = ulpError - mean;
+ diffSumSq += diff * diff;
+ }
+ }
+
+ final double stdDev = successCount > 1 ?
+ Math.sqrt(diffSumSq / (successCount - 1)) :
+ 0d;
+
+ return new Stats(mean, stdDev, min, max, ulpErrors.length - successCount);
+ }
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/NormsPerformance.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/NormsPerformance.java
new file mode 100644
index 0000000..2550530
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/NormsPerformance.java
@@ -0,0 +1,198 @@
+/*
+ * 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.arrays;
+
+import java.util.concurrent.TimeUnit;
+import java.util.function.ToDoubleFunction;
+
+import org.apache.commons.numbers.arrays.Norms;
+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;
+
+/**
+ * Execute benchmarks for the methods in the {@link Norms} 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 NormsPerformance {
+
+ /** Class providing input vectors for benchmarks.
+ */
+ @State(Scope.Benchmark)
+ public static class VectorArrayInput {
+
+ /** Number of vector samples. */
+ @Param("100000")
+ private int samples;
+
+ /** Minimum possible double exponent. */
+ @Param("-550")
+ private int minExp;
+
+ /** Maximum possible double exponent. */
+ @Param("+550")
+ private int maxExp;
+
+ /** Range of exponents within a single vector. */
+ @Param("26")
+ private int vectorExpRange;
+
+ /** Array of input vectors. */
+ private double[][] vectors;
+
+ /** Get the input vectors.
+ * @return input vectors
+ */
+ public double[][] getVectors() {
+ return vectors;
+ }
+
+ /** Create the input vectors for the instance.
+ */
+ @Setup
+ public void createVectors() {
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP);
+
+ vectors = new double[samples][];
+ for (int i = 0; i < vectors.length; ++i) {
+ // pick a general range for the vector element exponents and then
+ // create values within that range
+ final int vMidExp = rng.nextInt(maxExp - minExp + 1) + minExp;
+ final int vExpRadius = vectorExpRange / 2;
+ final int vMinExp = vMidExp - vExpRadius;
+ final int vMaxExp = vMidExp + vExpRadius;
+
+ vectors[i] = DoubleUtils.randomArray(getLength(), vMinExp, vMaxExp, rng);
+ }
+ }
+
+ /** Get the length of the input vectors.
+ * @return input vector length
+ */
+ protected int getLength() {
+ return 3;
+ }
+ }
+
+ /** Class providing 2D input vectors for benchmarks.
+ */
+ @State(Scope.Benchmark)
+ public static class VectorArrayInput2D extends VectorArrayInput {
+
+ /** {@inheritDoc} */
+ @Override
+ protected int getLength() {
+ return 2;
+ }
+ }
+
+ /** Evaluate a norm computation method with the given input.
+ * @param fn function to evaluate
+ * @param input computation input
+ * @param bh blackhole
+ */
+ private static void eval(final ToDoubleFunction<double[]> fn, final VectorArrayInput input,
+ final Blackhole bh) {
+ final double[][] vectors = input.getVectors();
+ for (int i = 0; i < vectors.length; ++i) {
+ bh.consume(fn.applyAsDouble(vectors[i]));
+ }
+ }
+
+ /** Compute the Euclidean norm directly with no checks for overflow or underflow.
+ * @param v input vector
+ * @return Euclidean norm
+ */
+ private static double directEuclideanNorm(final double[] v) {
+ double n = 0;
+ for (int i = 0; i < v.length; i++) {
+ n += v[i] * v[i];
+ }
+ return Math.sqrt(n);
+ }
+
+ /** Compute a baseline performance metric with a method that does nothing.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void baseline(final VectorArrayInput input, final Blackhole bh) {
+ eval(v -> 0d, input, bh);
+ }
+
+ /** Compute a baseline performance metric using direct computation of the
+ * Euclidean norm.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void directEuclideanArray(final VectorArrayInput input, final Blackhole bh) {
+ eval(NormsPerformance::directEuclideanNorm, input, bh);
+ }
+
+ /** Compute a baseline performance metric using {@link Math#hypot(double, double)}.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void hypot(final VectorArrayInput2D input, final Blackhole bh) {
+ eval(v -> Math.hypot(v[0], v[1]), input, bh);
+ }
+
+ /** Compute the performance of the {@link Norms#euclidean(double, double)} method.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void euclidean2d(final VectorArrayInput2D input, final Blackhole bh) {
+ eval(v -> Norms.euclidean(v[0], v[1]), input, bh);
+ }
+
+ /** Compute the performance of the {@link Norms#euclidean(double, double, double)} method.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void euclidean3d(final VectorArrayInput input, final Blackhole bh) {
+ eval(v -> Norms.euclidean(v[0], v[1], v[2]), input, bh);
+ }
+
+ /** Compute the performance of the {@link Norms#euclidean(double[])} method.
+ * @param input benchmark input
+ * @param bh blackhole
+ */
+ @Benchmark
+ public void euclideanArray(final VectorArrayInput input, final Blackhole bh) {
+ eval(Norms::euclidean, input, bh);
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAccuracyTest.java b/commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAccuracyTest.java
new file mode 100644
index 0000000..7d30ee2
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAccuracyTest.java
@@ -0,0 +1,141 @@
+/*
+ * 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.arrays;
+
+import java.io.BufferedWriter;
+import java.io.IOException;
+import java.nio.file.Files;
+import java.nio.file.Paths;
+import java.util.Map;
+
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.simple.RandomSource;
+import org.junit.jupiter.api.Disabled;
+import org.junit.jupiter.api.Test;
+
+/**
+ * Test the accuracy of the algorithms in the {@link EuclideanNorms} class.
+ */
+class EuclideanNormAccuracyTest {
+
+ /** Length of vectors to compute norms for. */
+ private static final int VECTOR_LENGTH = 100;
+
+ /** Number of samples per evaluation. */
+ private static final int SAMPLE_COUNT = 100_000;
+
+ /** Report the relative error of various Euclidean norm computation methods and write
+ * the results to a csv file. This is not a test.
+ * @throws IOException if an I/O error occurs
+ */
+ @Test
+ @Disabled("This method is used to output a report of the accuracy of implementations.")
+ void reportUlpErrors() throws IOException {
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP);
+
+ final EuclideanNormEvaluator eval = new EuclideanNormEvaluator();
+ eval.addMethod("direct", new EuclideanNormAlgorithms.Direct())
+ .addMethod("enorm", new EuclideanNormAlgorithms.Enorm())
+ .addMethod("enormMod", new EuclideanNormAlgorithms.EnormMod())
+ .addMethod("enormModKahan", new EuclideanNormAlgorithms.EnormModKahan())
+ .addMethod("enormModExt", new EuclideanNormAlgorithms.EnormModExt())
+ .addMethod("extLinear", new EuclideanNormAlgorithms.ExtendedPrecisionLinearCombination())
+ .addMethod("extLinearMod", new EuclideanNormAlgorithms.ExtendedPrecisionLinearCombinationMod())
+ .addMethod("extLinearSinglePass", new EuclideanNormAlgorithms.ExtendedPrecisionLinearCombinationSinglePass())
+ .addMethod("extLinearSqrt2", new EuclideanNormAlgorithms.ExtendedPrecisionLinearCombinationSqrt2());
+
+ try (BufferedWriter writer = Files.newBufferedWriter(Paths.get("target/norms.csv"))) {
+ writer.write("name, input type, error mean, error std dev, error min, error max, failed");
+ writer.newLine();
+
+ // construct a baseline array of random vectors
+ final double[][] baseInput = createInputVectors(-10, +10, rng);
+
+ evaluate("high", scaleInputVectors(baseInput, 0x1.0p+530), eval, writer);
+ evaluate("high-thresh", scaleInputVectors(baseInput, 0x1.0p+500), eval, writer);
+ evaluate("mid", baseInput, eval, writer);
+ evaluate("low-thresh", scaleInputVectors(baseInput, 0x1.0p-500), eval, writer);
+ evaluate("low", scaleInputVectors(baseInput, 0x1.0p-530), eval, writer);
+
+ final double[][] fullInput = createInputVectors(-550, +550, rng);
+ evaluate("full", fullInput, eval, writer);
+ }
+ }
+
+ /** Perform a single evaluation run and write the results to {@code writer}.
+ * @param inputType type of evaluation input
+ * @param inputs input vectors
+ * @param eval evaluator
+ * @param writer output writer
+ * @throws IOException if an I/O error occurs
+ */
+ private static void evaluate(final String inputType, final double[][] inputs, final EuclideanNormEvaluator eval,
+ final BufferedWriter writer) throws IOException {
+ final Map<String, EuclideanNormEvaluator.Stats> resultMap = eval.evaluate(inputs);
+ writeResults(inputType, resultMap, writer);
+ }
+
+ /** Write evaluation results to the given writer instance.
+ * @param inputType type of evaluation input
+ * @param statsMap evaluation result map
+ * @param writer writer instance
+ * @throws IOException if an I/O error occurs
+ */
+ private static void writeResults(final String inputType, final Map<String, EuclideanNormEvaluator.Stats> resultMap,
+ final BufferedWriter writer) throws IOException {
+ for (Map.Entry<String, EuclideanNormEvaluator.Stats> entry : resultMap.entrySet()) {
+ EuclideanNormEvaluator.Stats stats = entry.getValue();
+
+ writer.write(String.format("%s,%s,%.3g,%.3g,%.3g,%.3g,%d",
+ entry.getKey(),
+ inputType,
+ stats.getUlpErrorMean(),
+ stats.getUlpErrorStdDev(),
+ stats.getUlpErrorMin(),
+ stats.getUlpErrorMax(),
+ stats.getFailCount()));
+ writer.newLine();
+ }
+ }
+
+ /** Create an array of random input vectors with exponents in the range {@code [minExp, maxExp]}.
+ * @param minExp minimum exponent
+ * @param maxExp maximum exponent
+ * @param rng random number generator
+ * @return array of random input vectors
+ */
+ private static double[][] createInputVectors(final int minExp, final int maxExp, final UniformRandomProvider rng) {
+ final double[][] input = new double[SAMPLE_COUNT][];
+ for (int i = 0; i < input.length; ++i) {
+ input[i] = DoubleUtils.randomArray(VECTOR_LENGTH, minExp, maxExp, rng);
+ }
+ return input;
+ }
+
+ /** Return a copy of {@code inputs} with each element scaled by {@code s}.
+ * @param inputs input vectors
+ * @param s scale factor
+ * @return copy of {@code inputs} with each element scaled by {@code s}
+ */
+ private static double[][] scaleInputVectors(final double[][] inputs, final double s) {
+ final double[][] scaled = new double[inputs.length][];
+ for (int i = 0; i < inputs.length; ++i) {
+ scaled[i] = DoubleUtils.scalarMultiply(inputs[i], s);
+ }
+ return scaled;
+ }
+}
diff --git a/pom.xml b/pom.xml
index aa77692..76ededa 100644
--- a/pom.xml
+++ b/pom.xml
@@ -644,6 +644,15 @@
<groupId>org.eluder.coveralls</groupId>
<artifactId>coveralls-maven-plugin</artifactId>
<version>${commons.coveralls.version}</version>
+ <!-- Workaround for coveralls plugin JDK 9 compatibility issue.
+ https://github.com/trautonen/coveralls-maven-plugin/issues/112 -->
+ <dependencies>
+ <dependency>
+ <groupId>javax.xml.bind</groupId>
+ <artifactId>jaxb-api</artifactId>
+ <version>2.2.3</version>
+ </dependency>
+ </dependencies>
<configuration>
<timestampFormat>${commons.coveralls.timestampFormat}</timestampFormat>
</configuration>
diff --git a/src/main/resources/spotbugs/spotbugs-exclude-filter.xml b/src/main/resources/spotbugs/spotbugs-exclude-filter.xml
index 071037f..479aab1 100644
--- a/src/main/resources/spotbugs/spotbugs-exclude-filter.xml
+++ b/src/main/resources/spotbugs/spotbugs-exclude-filter.xml
@@ -41,7 +41,7 @@
<Method name="sqrt"/>
<BugPattern name="FE_FLOATING_POINT_EQUALITY"/>
</Match>
- <Match>
+ <Match>
<Class name="org.apache.commons.numbers.angle.Angle$Normalizer"/>
<Method name="applyAsDouble"/>
<BugPattern name="FE_FLOATING_POINT_EQUALITY"/>