aherbert commented on a change in pull request #92: URL: https://github.com/apache/commons-numbers/pull/92#discussion_r646129861
########## File path: commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Norms.java ########## @@ -0,0 +1,451 @@ +/* + * 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. Primary sources for the included algorithms are 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> and the + * <a href="http://www.netlib.org/minpack";>minpack</a> "enorm" subroutine. + * @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-500; Review comment: This can be lower. Perhaps 0x1.0p-511? ########## File path: commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Norms.java ########## @@ -0,0 +1,451 @@ +/* + * 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. Primary sources for the included algorithms are 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> and the + * <a href="http://www.netlib.org/minpack";>minpack</a> "enorm" subroutine. + * @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-500; + + /** 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; + if (xabs > yabs) { + 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)) { + if (Double.isNaN(x) || Double.isNaN(y)) { + return Double.NaN; + } + return Double.POSITIVE_INFINITY; + } 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)) { + if (Double.isNaN(x) || Double.isNaN(y) || Double.isNaN(z)) { + return Double.NaN; + } + return Double.POSITIVE_INFINITY; + } + + // 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_i^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|, ... |v_i|)}\), i.e., the maximum of the absolute values of the Review comment: This should be `\(\max{(|v_0| + ... + |v_i|)}\)` to match the style of the other javadocs. I rendered the javadocs and prefer the `\ldots`. Also use the length of the input n so the summation is to n-1. Edit to `\(\max{(|v_0| + \ldots + |v_{n-1}|)}\)`. ``` rm -f target/maven-javadoc-plugin-stale-data.txt mvn javadoc:javadoc open target/site/apidocs/org/apache/commons/numbers/arrays/Norms.html ``` ########## File path: commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/NormsPerformance.java ########## @@ -0,0 +1,194 @@ +/* + * 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.numbers.examples.jmh.DoubleUtils; +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.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 { + + /** Number of samples used in each benchmark. */ + private static final int SAMPLES = 100_000; + + /** Minimum random double exponent. */ + private static final int MIN_EXP = -550; + + /** Maximum random double exponent. */ + private static final int MAX_EXP = +550; + + /** Range of exponents within a single vector. */ + private static final int VECTOR_EXP_RANGE = 26; + + /** Class providing input vectors for benchmarks. + */ + @State(Scope.Benchmark) + public static class VectorArrayInput { + + /** 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(MAX_EXP - MIN_EXP) + MIN_EXP; Review comment: Invalid range again. Add a 1 to it. ########## File path: commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Norms.java ########## @@ -0,0 +1,451 @@ +/* + * 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. Primary sources for the included algorithms are 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> and the + * <a href="http://www.netlib.org/minpack";>minpack</a> "enorm" subroutine. + * @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-500; + + /** 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; + if (xabs > yabs) { + 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)) { + if (Double.isNaN(x) || Double.isNaN(y)) { + return Double.NaN; + } + return Double.POSITIVE_INFINITY; + } 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)) { + if (Double.isNaN(x) || Double.isNaN(y) || Double.isNaN(z)) { Review comment: A numerical fiddle would be to return: ```java return xabs * yabs * zabs; ``` ########## File path: commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/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; + +import org.apache.commons.rng.UniformRandomProvider; + +/** Class containing utility methods for working with doubles. + */ +public final class DoubleUtils { Review comment: Does all this have to be public? I think it can move to the arrays package and be package private. ########## File path: commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Norms.java ########## @@ -0,0 +1,451 @@ +/* + * 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. Primary sources for the included algorithms are 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> and the + * <a href="http://www.netlib.org/minpack";>minpack</a> "enorm" subroutine. Review comment: This class is no longer based on minpack enorm. I think the reference to it can be dropped. ########## File path: commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAlgorithmPerformance.java ########## @@ -0,0 +1,235 @@ +/* + * 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.examples.jmh.DoubleUtils; +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 { + + /** Default number of samples used in each benchmark. */ + private static final int DEFAULT_SAMPLES = 100_000; + + /** Default length of vectors used in the benchmarks. */ + private static final int DEFAULT_VECTOR_LEN = 100; + + /** 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"; + Review comment: Remove the extra line ########## File path: commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/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; + +import org.apache.commons.rng.UniformRandomProvider; + +/** Class containing utility methods for working with doubles. + */ +public final class DoubleUtils { + + /** No instantiation. */ + private DoubleUtils() {} + + /** Create a random double value with exponent in the range {@code [minExp, maxExp]}. + * @param minExp minimum exponent value + * @param maxExp maximum exponent value + * @param rng random number generator + * @return random double + */ + public 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(Math.abs(maxExp - minExp)) + minExp + 1023; Review comment: Again the range should have a +1 to avoid calling nextInt with the wrong limit. ########## File path: commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/ExtendedPrecisionTest.java ########## @@ -102,6 +102,23 @@ void testSubNormalSplit() { 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(-1.4293872661474477E-16, Math.PI); Review comment: Can you compute this result using BigDecimal? ```java final BigDecimal pi = new BigDecimal(Math.PI); final double pi2low = pi.pow(2).subtract(new BigDecimal(Math.PI * Math.PI)).doubleValue(); assertSquareLowUnscaled(pi2low, Math.PI); ``` It is a bit cryptic to have a constant here. ########## File path: commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Norms.java ########## @@ -0,0 +1,451 @@ +/* + * 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. Primary sources for the included algorithms are 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> and the + * <a href="http://www.netlib.org/minpack";>minpack</a> "enorm" subroutine. + * @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-500; + + /** 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; + if (xabs > yabs) { + 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)) { Review comment: If x is NaN and y is finite then (xabs > yabs) is false and max will be yabs. Here nan detection will fail. This can lead to you running the entire calculation with one number as NaN. You still get a NaN result but do a lot of work to get there. These tests pass: ```java Assertions.assertEquals(Double.NaN, Norms.euclidean(Double.NaN, 0d)); Assertions.assertEquals(Double.NaN, Norms.euclidean(Double.NaN, Double.MAX_VALUE)); ``` But the first one returns NaN due to the detection of a large gap in the exponents. The second performs the whole calculation and detects NaN at the end. Since a few lines down you extract the exponents then perhaps you should extract the exponents first. You then find the max based on the exponent. The mantissa is not important for ordering as the max is just for scaling detection. Or extract the entire raw long bits and remove the sign bit: ```java long bitsa = Double.doubleToRawLongBits(x) & Long.MAX_VALUE; long bitsb = Double.doubleToRawLongBits(y) & Long.MAX_VALUE; // Swap so a is max; if (bitsb > bitsa) { final long tmp = bitsa; bitsa = bitsb; bitsb = tmp; } if ((bitsa >>> 52) == 2047) { // NaN/inf - Could be combined with the condition below as we return the same result return Double.longBitsToDouble(bitsa); } else if (bitsa - bitsb > (EXP_DIFF_THRESHOLD_2D << 52)) { return Double.longBitsToDouble(bitsa); } double xabs = Double.longBitsToDouble(bitsa); double yabs = Double.longBitsToDouble(bitsb); // Compute ``` NaN and inf will have the highest exponent (or NaN will be highest if using the raw bits). If this is found then you just return NaN or inf as appropriate. For a result similar to the 3D case then a simpler fix if to change to: ```java // NaN will be highest if (Double.compare(xabs, yabs) > 0) { ``` This is more readable. ########## File path: commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/Norms.java ########## @@ -0,0 +1,451 @@ +/* + * 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. Primary sources for the included algorithms are 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> and the + * <a href="http://www.netlib.org/minpack";>minpack</a> "enorm" subroutine. + * @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-500; + + /** 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; + if (xabs > yabs) { + 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)) { + if (Double.isNaN(x) || Double.isNaN(y)) { + return Double.NaN; + } + return Double.POSITIVE_INFINITY; + } 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); Review comment: Possibly use `maximum(x, y, z)`. I do not think the absolutes are required for the rest of the computation. ########## File path: commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAlgorithmPerformance.java ########## @@ -0,0 +1,235 @@ +/* + * 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.examples.jmh.DoubleUtils; +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 { + + /** Default number of samples used in each benchmark. */ + private static final int DEFAULT_SAMPLES = 100_000; Review comment: This is only used in one place. Just hard code it there and avoid the ugly `int + ""` declaration in the `@Param` tag. Same for DEFAULT_VECTOR_LEN. ########## File path: commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAccuracyTest.java ########## @@ -0,0 +1,142 @@ +/* + * 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.numbers.examples.jmh.DoubleUtils; +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]}.. Review comment: Two periods at the end of the line. ########## File path: commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormEvaluator.java ########## @@ -0,0 +1,248 @@ +/* + * 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. This method assumes that + * the arguments have the same exponent. Review comment: They do not have to have the same exponent. The ulp difference between 1.0 and Math.nextDown(1.0) will be 1. They do have to have the same sign. All your results should be positive so this is OK. If the numbers are very different then this will return a truncated difference by long to int primitive conversion. A return type of `long` would fix this. I do not think that you require `int` during testing as the sums use a double. This comment also applies to `DoubleTestUtils` in the main package. ########## File path: commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/DoubleTestUtils.java ########## @@ -0,0 +1,70 @@ +/* + * 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 the arguments. The arguments are assumed to + * have the same exponent. + * @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 + * @param maxExp maximum exponent + * @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 = Math.abs(maxExp - minExp); Review comment: Should this be Math.abs(maxExp - minExp) + 1? Otherwise the call to rng.nextInt(int) could be called with 0 when max is the same as min. ########## File path: commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/ExtendedPrecisionTest.java ########## @@ -102,6 +102,23 @@ void testSubNormalSplit() { final double lo2 = a - hi2; Assertions.assertEquals(a, hi2); Assertions.assertEquals(0, lo2); + Review comment: Remove extra line ########## File path: commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/NormsPerformance.java ########## @@ -0,0 +1,194 @@ +/* + * 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.numbers.examples.jmh.DoubleUtils; +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.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 { + + /** Number of samples used in each benchmark. */ + private static final int SAMPLES = 100_000; Review comment: Can all these be put into `@Param` annotated attributes in VectorArrayInput. This allows them to be adjusted by passing values to JMH on the command line. ########## File path: commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormAlgorithmPerformance.java ########## @@ -0,0 +1,235 @@ +/* + * 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.examples.jmh.DoubleUtils; +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 { + + /** Default number of samples used in each benchmark. */ + private static final int DEFAULT_SAMPLES = 100_000; + + /** Default length of vectors used in the benchmarks. */ + private static final int DEFAULT_VECTOR_LEN = 100; + + /** 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(DEFAULT_SAMPLES + "") + private int samples; + + /** The length of each vector. */ + @Param(DEFAULT_VECTOR_LEN + "") + 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: + minExp = -550; + maxExp = +550; + break; Review comment: I prefer here to throw in the default case statement. So if you pass parameters to JMH on the command line it will fail: ``` mvn package -Pexamples-jmh java -jar target/examples-jmh.jar EuclideanNormAlgorithmPerformance -p vectorLength=10 -p type=foobar ``` This should error but instead runs with `type=FULL`. ########## File path: commons-numbers-arrays/src/test/java/org/apache/commons/numbers/arrays/NormsTest.java ########## @@ -0,0 +1,564 @@ +/* + * 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 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(-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, 1, fn); + checkScaledEuclideanNorm(ones, 0x1.0p496, fn); + checkScaledEuclideanNorm(ones, 0x1.0p497, fn); + checkScaledEuclideanNorm(ones, 0x1.0p-100, fn); + checkScaledEuclideanNorm(ones, 0x1.0p-101, fn); + checkScaledEuclideanNorm(ones, 0x1.0p-500, fn); + checkScaledEuclideanNorm(ones, 0x1.0p-501, fn); + + + checkScaledEuclideanNorm(multiplesOfTen, 1, fn); + checkScaledEuclideanNorm(multiplesOfTen, 0x1.0p-100, fn); + checkScaledEuclideanNorm(multiplesOfTen, 0x1.0p-101, fn); + checkScaledEuclideanNorm(multiplesOfTen, 0x1.0p495, fn); + checkScaledEuclideanNorm(multiplesOfTen, 0x1.0p-500, 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(+490, +510, samples, rng); + assertEuclidean2dVersusHypot(-510, -490, 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 double hypotResult = computeUlpErrorStdDev(inputs, exactResults, v -> Math.hypot(v[0], v[1])); + final double normResult = computeUlpErrorStdDev(inputs, exactResults, v -> Norms.euclidean(v[0], v[1])); + + // ensure that we are within the ballpark of Math.hypot + Assertions.assertTrue(normResult <= (hypotResult + HYPOT_COMPARE_EPS), Review comment: You should also check the mean. What if Norms computed the exact result but shifted by 20 ULP. It would have a standard deviation of zero but a mean of 20. It would be very precise in its inaccuracy. Effectively you should test accuracy (mean ulp) and precision (std dev of the ulp). ########## File path: commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/EuclideanNormEvaluator.java ########## @@ -0,0 +1,248 @@ +/* + * 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. This method assumes that + * the arguments have the same exponent. + * @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 = Math.sqrt(diffSumSq / (successCount - 1)); Review comment: In the ridiculous case where success count is 1 this will be invalid. zero successes will compute a NaN mean which is a useful marker, and stdDev will be zero. -- This is an automated message from the Apache Git Service. To respond to the message, please log on to GitHub and use the URL above to go to the specific comment. For queries about this service, please contact Infrastructure at: [email protected]
