This is an automated email from the ASF dual-hosted git repository.
erans pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit a77a8eb22b9527c50e50eac106a6821676a63fba
Author: Gilles Sadowski <gillese...@gmail.com>
AuthorDate: Sun Jun 13 00:16:58 2021 +0200
NUMBERS-156: Norm API.
---
.../org/apache/commons/numbers/angle/CosAngle.java | 4 +-
.../org/apache/commons/numbers/core/Norms.java | 448 ---------------
.../org/apache/commons/numbers/core/NormsTest.java | 625 ---------------------
3 files changed, 2 insertions(+), 1075 deletions(-)
diff --git
a/commons-numbers-angle/src/main/java/org/apache/commons/numbers/angle/CosAngle.java
b/commons-numbers-angle/src/main/java/org/apache/commons/numbers/angle/CosAngle.java
index cfa569d..d9d7dd4 100644
---
a/commons-numbers-angle/src/main/java/org/apache/commons/numbers/angle/CosAngle.java
+++
b/commons-numbers-angle/src/main/java/org/apache/commons/numbers/angle/CosAngle.java
@@ -17,7 +17,7 @@
package org.apache.commons.numbers.angle;
import org.apache.commons.numbers.core.LinearCombination;
-import org.apache.commons.numbers.core.Norms;
+import org.apache.commons.numbers.core.Norm;
/**
* Computes the cosine of the angle between two vectors.
@@ -39,7 +39,7 @@ public final class CosAngle {
*/
public static double value(double[] v1,
double[] v2) {
- return LinearCombination.value(v1, v2) / Norms.euclidean(v1) /
Norms.euclidean(v2);
+ return LinearCombination.value(v1, v2) / Norm.L2.of(v1) /
Norm.L2.of(v2);
}
}
diff --git
a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Norms.java
b/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Norms.java
deleted file mode 100644
index 2a357ca..0000000
---
a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/Norms.java
+++ /dev/null
@@ -1,448 +0,0 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-package org.apache.commons.numbers.core;
-
-/** Class providing 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-core/src/test/java/org/apache/commons/numbers/core/NormsTest.java
b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/NormsTest.java
deleted file mode 100644
index aa1496d..0000000
---
a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/NormsTest.java
+++ /dev/null
@@ -1,625 +0,0 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-package org.apache.commons.numbers.core;
-
-import java.math.BigDecimal;
-import java.math.MathContext;
-import 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;
-
- private static final BigDecimal BD_MAX_VALUE = new
BigDecimal(Double.MAX_VALUE);
- private static final BigDecimal BD_MIN_NORMAL = new
BigDecimal(Double.MIN_NORMAL);
-
- /** The scale, used to scale the sqrt of the sum of squares. */
- private static final double SCALE = 0x1.0p200;
-
- /** The scale squared, used to scale the sum of squares. */
- private static final BigDecimal SCALE2 = new BigDecimal(SCALE * SCALE);
-
- @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, 3L);
-
- // 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(new BigDecimal(d).pow(2), ctx);
- }
-
- // Java 9+:
- // sum.sqrt(ctx).doubleValue()
-
- // Require the sum to be in the range of a double for conversion
before sqrt().
- // We scale by a power of 2. Rescaling uses the square root of this
which is also
- // a power of 2 and can be accumulated for exact rescaling.
- double rescale = 1.0;
- if (sum.compareTo(BD_MIN_NORMAL) < 0) {
- while (sum.compareTo(BD_MIN_NORMAL) < 0) {
- sum = sum.multiply(SCALE2);
- rescale /= SCALE;
- }
- } else if (sum.compareTo(BD_MAX_VALUE) > 0) {
- while (sum.compareTo(BD_MAX_VALUE) > 0) {
- sum = sum.divide(SCALE2);
- rescale *= SCALE;
- }
- }
-
- return Math.sqrt(sum.doubleValue()) * rescale;
- }
-
- /** 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;
- }
- }
-}
