You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2020/04/06 12:34:03 UTC
[commons-numbers] 01/03: [NUMBERS-142] Added LinearCombination
implementations to the JMH module.
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit 0494360c0eb9f7df75db2752cd4d094e4b880853
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Thu Jan 16 11:32:48 2020 +0000
[NUMBERS-142] Added LinearCombination implementations to the JMH module.
Adds implementations based on the method of Dekker, Ogita et al, and
Shewchuk.
All implementations pass the tests copied from
arrays.LinearCombinationTest.
Tests were added using ill-conditioned dot products to demonstrate the
condition number where each implementation will fail.
The method for computing a double length product of two numbers has been
written to ensure no overflow of intermediates.
All implementations are added to a JMH performance test.
---
commons-numbers-angle/LICENSE | 30 +
commons-numbers-examples/examples-jmh/pom.xml | 10 +
.../examples/jmh/arrays/DoublePrecision.java | 705 +++++++
.../jmh/arrays/DoubleSplitPerformance.java | 615 ++++++
.../examples/jmh/arrays/LinearCombination.java | 98 +
.../jmh/arrays/LinearCombinationPerformance.java | 371 ++++
.../jmh/arrays/LinearCombinationUtils.java | 211 ++
.../examples/jmh/arrays/LinearCombinations.java | 2007 ++++++++++++++++++++
.../examples/jmh/arrays/StickySumPerformance.java | 609 ++++++
.../numbers/examples/jmh/arrays/package-info.java | 21 +
.../examples/jmh/arrays/DoublePrecisionTest.java | 93 +
.../jmh/arrays/LinearCombinationAccuracyTest.java | 241 +++
.../jmh/arrays/LinearCombinationsTest.java | 588 ++++++
.../checkstyle/checkstyle-suppressions.xml | 4 +-
14 files changed, 5602 insertions(+), 1 deletion(-)
diff --git a/commons-numbers-angle/LICENSE b/commons-numbers-angle/LICENSE
index 261eeb9..4cc610a 100644
--- a/commons-numbers-angle/LICENSE
+++ b/commons-numbers-angle/LICENSE
@@ -199,3 +199,33 @@
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.
+
+================================================================================
+
+Class "org.apache.commons.numbers.examples.jmh.arrays.StickySumPerformance" contains
+Java code ported from "Free BSD" in C for benchmarking.
+The source files contain the following notice:
+
+ * Copyright (c) 2005-2011 David Schultz <da...@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
diff --git a/commons-numbers-examples/examples-jmh/pom.xml b/commons-numbers-examples/examples-jmh/pom.xml
index 544c24f..4c7ab41 100644
--- a/commons-numbers-examples/examples-jmh/pom.xml
+++ b/commons-numbers-examples/examples-jmh/pom.xml
@@ -33,6 +33,11 @@
<dependencies>
<dependency>
<groupId>org.apache.commons</groupId>
+ <artifactId>commons-numbers-arrays</artifactId>
+ </dependency>
+
+ <dependency>
+ <groupId>org.apache.commons</groupId>
<artifactId>commons-numbers-complex</artifactId>
</dependency>
@@ -43,6 +48,11 @@
<dependency>
<groupId>org.apache.commons</groupId>
+ <artifactId>commons-numbers-fraction</artifactId>
+ </dependency>
+
+ <dependency>
+ <groupId>org.apache.commons</groupId>
<artifactId>commons-math3</artifactId>
</dependency>
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/DoublePrecision.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/DoublePrecision.java
new file mode 100644
index 0000000..b929b70
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/DoublePrecision.java
@@ -0,0 +1,705 @@
+/*
+ * 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;
+
+/**
+ * Computes double-length precision floating-point operations.
+ *
+ * <p>It is based on the 1971 paper
+ * <a href="https://doi.org/10.1007/BF01397083">
+ * Dekker (1971) A floating-point technique for extending the available precision</a>.
+ */
+final class DoublePrecision {
+ /*
+ * Caveat:
+ *
+ * The code below uses many additions/subtractions that may
+ * appear redundant. However, they should NOT be simplified, as they
+ * do use IEEE754 floating point arithmetic rounding properties.
+ *
+ * Algorithms are based on computing the product or sum of two values x and y in
+ * extended precision. The standard result is stored using a double (high part z) and
+ * the round-off error (or low part zz) is stored in a second double, e.g:
+ * x * y = (z, zz); z + zz = x * y
+ * x + y = (z, zz); z + zz = x + y
+ *
+ * To sum multiple (z, zz) results ideally the parts are sorted in order of
+ * non-decreasing magnitude and summed. This is exact if each number's most significant
+ * bit is below the least significant bit of the next (i.e. does not
+ * overlap). Creating non-overlapping parts requires a rebalancing
+ * of adjacent pairs using a summation z + zz = (z1, zz1) iteratively through the parts
+ * (see Shewchuk (1997) Grow-Expansion and Expansion-Sum [1]).
+ *
+ * [1] Shewchuk (1997): Arbitrary Precision Floating-Point Arithmetic
+ * http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps
+ */
+
+ /**
+ * The multiplier used to split the double value into high and low parts. From
+ * Dekker (1971): "The constant should be chosen equal to 2^(p - p/2) + 1,
+ * where p is the number of binary digits in the mantissa". Here p is 53
+ * and the multiplier is {@code 2^27 + 1}.
+ */
+ private static final double MULTIPLIER = 1.34217729E8;
+
+ /** The upper limit above which a number may overflow during the split into a high part.
+ * Assuming the multiplier is above 2^27 and the maximum exponent is 1023 then a safe
+ * limit is a value with an exponent of (1023 - 27) = 2^996. */
+ private static final double SAFE_UPPER = 0x1.0p996;
+
+ /** The scale to use when down-scaling during a split into a high part.
+ * This must be smaller than the inverse of the multiplier and a power of 2 for exact scaling. */
+ private static final double DOWN_SCALE = 0x1.0p-30;
+
+ /** The scale to use when re-scaling during a split into a high part.
+ * This is the inverse of {@link #DOWN_SCALE}. */
+ private static final double UP_SCALE = 0x1.0p30;
+
+ /** The mask to zero the lower 27-bits of a long . */
+ private static final long ZERO_LOWER_27_BITS = 0xffff_ffff_f800_0000L;
+
+ /** The mask to extract the raw 11-bit exponent.
+ * The value must be shifted 52-bits to remove the mantissa bits. */
+ private static final int EXP_MASK = 0x7ff;
+
+ /** The value 2046 converted for use if using {@link Integer#compareUnsigned(int, int)}.
+ * This requires adding {@link Integer#MIN_VALUE} to 2046. */
+ private static final int CMP_UNSIGNED_2046 = Integer.MIN_VALUE + 2046;
+
+ /** The value -1 converted for use if using {@link Integer#compareUnsigned(int, int)}.
+ * This requires adding {@link Integer#MIN_VALUE} to -1. */
+ private static final int CMP_UNSIGNED_MINUS_1 = Integer.MIN_VALUE - 1;
+
+ /**
+ * Represents a floating-point number with twice the precision of a {@code double}.
+ */
+ static final class Quad {
+ // This is treated as a simple struct.
+ // CHECKSTYLE: stop VisibilityModifier
+ /** The high part of the number. */
+ double hi;
+ /** The low part of the number. */
+ double lo;
+ // CHECKSTYLE: resume VisibilityModifier
+ }
+
+ /** Private constructor. */
+ private DoublePrecision() {
+ // intentionally empty.
+ }
+
+ /**
+ * Multiply the values {@code x} and {@code y} into a double-precision result {@code z}.
+ * It is assumed the numbers are normalized so no over/underflow will occurs.
+ *
+ * <p>Implements Dekker's mul12 method to split the numbers and multiply them
+ * in extended precision.
+ *
+ * <p>Note: The quad satisfies the condition {@code x * y == z.hi + z.lo}. The high
+ * part may be different from {@code x * y} by 1 ulp due to rounding.
+ *
+ * @param x First value
+ * @param y Second value
+ * @param z Result
+ */
+ static void multiplyUnscaled(double x, double y, Quad z) {
+ // Note: The original mul12 algorithm avoids x * y and saves 1 multiplication.
+ double p;
+ p = x * MULTIPLIER;
+ final double hx = x - p + p;
+ final double lx = x - hx;
+ p = y * MULTIPLIER;
+ final double hy = y - p + p;
+ final double ly = y - hy;
+ p = hx * hy;
+ final double q = hx * ly + lx * hy;
+ z.hi = p + q;
+ z.lo = p - z.hi + q + lx * ly;
+ }
+
+ /**
+ * Multiply the values {@code x} and {@code y} into a double-precision result {@code c}.
+ * Scaling is performed on the numbers to protect against intermediate over/underflow.
+ *
+ * <p>The quadruple precision result has the standard double precision result
+ * {@code x * y} in the high part and the round-off in the low part,
+ *
+ * <p>Special cases:
+ *
+ * <ul>
+ * <li>If {@code x * y} is sub-normal or zero then the low part is 0.0.
+ * <li>If {@code x * y} is infinite or NaN then the low part is NaN.
+ * </ul>
+ *
+ * <p>Note: This does not represent the low part of infinity with zero. This is because the
+ * method is intended to be used for extended precision computations. The NaN low part
+ * signals that an extended precision computation using the result is invalid (i.e. the
+ * result of summation/multiplication of the parts will not be finite).
+ *
+ * @param x First value
+ * @param y Second value
+ * @param c Result
+ * @see DoublePrecision#productLowUnscaled(double, double, double)
+ */
+ static void multiply(double x, double y, Quad c) {
+ // Special cases. Check the product.
+ final double xy = x * y;
+ if (isNotNormal(xy)) {
+ c.hi = xy;
+ // Returns 0.0 for sub-normal xy, otherwise NaN for inf/nan
+ c.lo = xy - xy;
+ return;
+ }
+ // Extract biased exponent and normalise.
+ // Sub-normals are scaled by 2^54 and the exponent adjusted.
+ // This is equivalent to the c function frexp which decomposes given floating
+ // point value arg into a normalized fraction and an integral power of two.
+ // Here we use a biased exponent as it is later adjusted when re-scaling.
+ long xb = Double.doubleToRawLongBits(x);
+ int xe = getBiasedExponent(xb);
+ double xs;
+ if (xe == 0) {
+ // Sub-normal. Scale up and extract again
+ xs = x * 0x1.0p54;
+ xb = Double.doubleToRawLongBits(xs);
+ xe = getBiasedExponent(xb) - 54;
+ }
+ xs = getNormalisedFraction(xb);
+
+ long yb = Double.doubleToRawLongBits(y);
+ int ye = getBiasedExponent(yb);
+ double ys;
+ if (ye == 0) {
+ // Sub-normal. Scale up and extract again
+ ys = y * 0x1.0p54;
+ yb = Double.doubleToRawLongBits(ys);
+ ye = getBiasedExponent(yb) - 54;
+ }
+ ys = getNormalisedFraction(yb);
+
+ // Compute hi as x*y.
+ // Thus if the standard precision result is finite (as verified in the initial test
+ // on x * y) then the extended precision result will be.
+ double z = xs * ys;
+ double zz = productLowUnscaled(xs, ys, z);
+
+ // Re-scale. The result is currently in the range [0.25, 1) so no checks for
+ // 0, nan, inf (the result exponent will be -2 or -1).
+ // Both exponents are currently biased so subtract 1023 to get the biased scale.
+ int scale = xe + ye - 1023;
+ // Compute scaling by multiplication so we can scale both together.
+ // If a single multiplication to a normal number then handle here.
+ if (scale <= 2046 && scale > 0) {
+ // Convert to a normalised power of 2
+ final double d = Double.longBitsToDouble(((long) scale) << 52);
+ z *= d;
+ zz *= d;
+ } else {
+ // Delegate to java.util.Math
+ // We have to adjust the biased scale to unbiased using the exponent offset 1023.
+ scale -= 1023;
+ z = Math.scalb(z, scale);
+ zz = Math.scalb(zz, scale);
+ }
+
+ // Final result. The hi part should be same as the IEEE754 result.
+ // assert z == xy;
+ c.hi = z;
+ c.lo = zz;
+ }
+
+ /**
+ * Checks if the number is not normal. This is functionally equivalent to:
+ * <pre>
+ * final double abs = Math.abs(a);
+ * return (abs <= Double.MIN_NORMAL || !(absXy <= Double.MAX_VALUE));
+ * </pre>
+ *
+ * @param a The value.
+ * @return true if the value is not normal
+ */
+ static boolean isNotNormal(double a) {
+ // Sub-normal numbers have a biased exponent of 0.
+ // Inf/NaN numbers have a biased exponent of 2047.
+ // Catch both cases by extracting the raw exponent, subtracting 1
+ // and compare unsigned (so 0 underflows to a large value).
+ final int baisedExponent = ((int) (Double.doubleToRawLongBits(a) >>> 52)) & EXP_MASK;
+ // Pre-compute the additions used by Integer.compareUnsigned
+ return baisedExponent + CMP_UNSIGNED_MINUS_1 >= CMP_UNSIGNED_2046;
+ }
+
+ /**
+ * Gets the exponent.
+ *
+ * @param bits the bits
+ * @return the exponent
+ */
+ private static int getBiasedExponent(long bits) {
+ return (int)(bits >>> 52) & 0x7ff;
+ }
+
+ /**
+ * Gets the normalised fraction in the range [0.5, 1).
+ *
+ * @param bits the bits
+ * @return the exponent
+ */
+ private static double getNormalisedFraction(long bits) {
+ // Mask out the exponent and set it to 1022.
+ return Double.longBitsToDouble((bits & 0x800f_ffff_ffff_ffffL) | 0x3fe0_0000_0000_0000L);
+ }
+
+ /**
+ * Implement Dekker's method to split a value into two parts. Multiplying by (2^s + 1) creates
+ * a big value from which to derive the two split parts.
+ * <pre>
+ * c = (2^s + 1) * a
+ * a_big = c - a
+ * a_hi = c - a_big
+ * a_lo = a - a_hi
+ * a = a_hi + a_lo
+ * </pre>
+ *
+ * <p>The multiplicand allows a p-bit value to be split into
+ * (p-s)-bit value {@code a_hi} and a non-overlapping (s-1)-bit value {@code a_lo}.
+ * Combined they have (p-1) bits of significand but the sign bit of {@code a_lo}
+ * contains a bit of information. The constant is chosen so that s is ceil(p/2) where
+ * the precision p for a double is 53-bits (1-bit of the mantissa is assumed to be
+ * 1 for a non sub-normal number) and s is 27.
+ *
+ * <p>This conversion uses scaling to avoid overflow in intermediate computations.
+ *
+ * <p>Splitting a NaN or infinite value will return NaN. Any finite value will return
+ * a finite value.
+ *
+ * @param value Value.
+ * @return the high part of the value.
+ */
+ static double highPart(double value) {
+ // Avoid overflow
+ if (Math.abs(value) >= SAFE_UPPER) {
+ // Do scaling.
+ final double hi = highPartUnscaled(value * DOWN_SCALE) * UP_SCALE;
+ if (Double.isInfinite(hi)) {
+ // Number is too large.
+ // This occurs if value is infinite or close to Double.MAX_VALUE.
+ // Note that Dekker's split creates an approximating 26-bit number which may
+ // have an exponent 1 greater than the input value. This will overflow if the
+ // exponent is already +1023. Revert to the raw upper 26 bits of the 53-bit
+ // mantissa (including the assumed leading 1 bit). This conversion will result in
+ // the low part being a 27-bit significand and the potential loss of bits during
+ // addition and multiplication. (Contrast to the Dekker split which creates two
+ // 26-bit numbers with a bit of information moved to the sign of low.)
+ // The conversion will maintain Infinite in the high part where the resulting
+ // low part a_lo = a - a_hi = inf - inf = NaN.
+ return highPartSplit(value);
+ }
+ return hi;
+ }
+ // normal conversion
+ return highPartUnscaled(value);
+ }
+
+ /**
+ * Implement Dekker's method to split a value into two parts (see {@link #highPart(double)}).
+ *
+ * <p>This conversion does not use scaling and the result of overflow is NaN. Overflow
+ * may occur when the exponent of the input value is above 996.
+ *
+ * <p>Splitting a NaN or infinite value will return NaN.
+ *
+ * @param value Value.
+ * @return the high part of the value.
+ * @see Math#getExponent(double)
+ */
+ static double highPartUnscaled(double value) {
+ final double c = MULTIPLIER * value;
+ return c - (c - value);
+ }
+
+ /**
+ * Implement a split using the upper and lower raw bits from the value.
+ *
+ * <p>Note: This method will not work for very small sub-normal numbers
+ * ({@code <= 27} bits) as the high part will be zero and the low part will
+ * have all the information. Methods that assume {@code hi > lo} will have
+ * undefined behaviour.
+ *
+ * <p>Splitting a NaN value will return NaN or infinite. Splitting an infinite
+ * value will return infinite. Any finite value will return a finite value.
+ *
+ * @param value Value.
+ * @return the high part of the value.
+ */
+ static double highPartSplit(double value) {
+ return Double.longBitsToDouble(Double.doubleToRawLongBits(value) & ZERO_LOWER_27_BITS);
+ }
+
+ /**
+ * Compute the low part of the double length number {@code (z,zz)} for the exact
+ * product of {@code x} and {@code y}. This is equivalent to computing a {@code double}
+ * containing the magnitude of the rounding error when converting the exact 106-bit
+ * significand of the multiplication result to a 53-bit significand.
+ *
+ * <p>The method is written to be functionally similar to using a fused multiply add (FMA)
+ * operation to compute the low part, for example JDK 9's Math.fma function (note the sign
+ * change in the input argument for the product):
+ * <pre>
+ * double x = ...;
+ * double y = ...;
+ * double xy = x * y;
+ * double low1 = Math.fma(x, y, -xy);
+ * double low2 = DoublePrecision.productLow(x, y, xy);
+ * </pre>
+ *
+ * <p>Special cases:
+ *
+ * <ul>
+ * <li>If {@code x * y} is sub-normal or zero then the result is 0.0.
+ * <li>If {@code x * y} is infinite or NaN then the result is NaN.
+ * </ul>
+ *
+ * @param x First factor.
+ * @param y Second factor.
+ * @param xy Product of the factors (x * y).
+ * @return the low part of the product double length number
+ * @see #highPart(double)
+ * @see #productLow(double, double, double, double, double)
+ */
+ static double productLow(double x, double y, double xy) {
+ // Verify the input. This must be NaN safe.
+ //assert Double.compare(x * y, xy) == 0
+
+ // If the number is sub-normal, inf or nan there is no round-off.
+ if (isNotNormal(xy)) {
+ // Returns 0.0 for sub-normal xy, otherwise NaN for inf/nan:
+ return xy - xy;
+ }
+
+ // The result xy is finite and normal.
+ // Use Dekker's mul12 algorithm that splits the values into high and low parts.
+ // Dekker's split using multiplication will overflow if the value is within 2^27
+ // of double max value. It can also produce 26-bit approximations that are larger
+ // than the input numbers for the high part causing overflow in hx * hy when
+ // x * y does not overflow. So we must scale down big numbers.
+ // We only have to scale the largest number as we know the product does not overflow
+ // (if one is too big then the other cannot be).
+ // We also scale if the product is close to overflow to avoid intermediate overflow.
+ // This could be done at a higher limit (e.g. Math.abs(xy) > Double.MAX_VALUE / 4)
+ // but is included here to have a single low probability branch condition.
+
+ // Add the absolute inputs for a single comparison. The sum will not be more than
+ // 3-fold higher than any component.
+ final double a = Math.abs(x);
+ final double b = Math.abs(y);
+ if (a + b + Math.abs(xy) >= SAFE_UPPER) {
+ // Only required to scale the largest number as x*y does not overflow.
+ if (a > b) {
+ return productLowUnscaled(x * DOWN_SCALE, y, xy * DOWN_SCALE) * UP_SCALE;
+ }
+ return productLowUnscaled(x, y * DOWN_SCALE, xy * DOWN_SCALE) * UP_SCALE;
+ }
+
+ // No scaling required
+ return productLowUnscaled(x, y, xy);
+ }
+
+ /**
+ * Compute the low part of the double length number {@code (z,zz)} for the exact
+ * product of {@code x} and {@code y}. This is equivalent to computing a {@code double}
+ * containing the magnitude of the rounding error when converting the exact 106-bit
+ * significand of the multiplication result to a 53-bit significand.
+ *
+ * <p>Special cases:
+ *
+ * <ul>
+ * <li>If {@code x * y} is sub-normal or zero then the result is 0.0.
+ * <li>If {@code x * y} is infinite, and {@code x} and {@code y} are finite then the
+ * result is the opposite infinity.
+ * <li>If {@code x} or {@code y} are infinite then the result is NaN.
+ * <li>If {@code x * y} is NaN then the result is NaN.
+ * </ul>
+ *
+ * @param x First factor.
+ * @param y Second factor.
+ * @param xy Product of the factors (x * y).
+ * @return the low part of the product double length number
+ * @see #highPart(double)
+ * @see #productLow(double, double, double, double, double)
+ */
+ static double productLow1(double x, double y, double xy) {
+ // Verify the input. This must be NaN safe.
+ //assert Double.compare(x * y, xy) == 0
+
+ // Logic as per productLow but with no check for sub-normal or NaN.
+ final double a = Math.abs(x);
+ final double b = Math.abs(y);
+ if (a + b + Math.abs(xy) >= SAFE_UPPER) {
+ // Only required to scale the largest number as x*y does not overflow.
+ if (a > b) {
+ return productLowUnscaled(x * DOWN_SCALE, y, xy * DOWN_SCALE) * UP_SCALE;
+ }
+ return productLowUnscaled(x, y * DOWN_SCALE, xy * DOWN_SCALE) * UP_SCALE;
+ }
+
+ // No scaling required
+ return productLowUnscaled(x, y, xy);
+ }
+ /**
+ * Compute the low part of the double length number {@code (z,zz)} for the exact
+ * product of {@code x} and {@code y}. This is equivalent to computing a {@code double}
+ * containing the magnitude of the rounding error when converting the exact 106-bit
+ * significand of the multiplication result to a 53-bit significand.
+ *
+ * <p>The method is written to be functionally similar to using a fused multiply add (FMA)
+ * operation to compute the low part, for example JDK 9's Math.fma function (note the sign
+ * change in the input argument for the product):
+ * <pre>
+ * double x = ...;
+ * double y = ...;
+ * double xy = x * y;
+ * double low1 = Math.fma(x, y, -xy);
+ * double low2 = DoublePrecision.productLow(x, y, xy);
+ * </pre>
+ *
+ * <p>Special cases:
+ *
+ * <ul>
+ * <li>If {@code x * y} is sub-normal or zero then the result is 0.0.
+ * <li>If {@code x * y} is infinite or NaN then the result is NaN.
+ * </ul>
+ *
+ * @param x First factor.
+ * @param y Second factor.
+ * @param xy Product of the factors (x * y).
+ * @return the low part of the product double length number
+ * @see #highPart(double)
+ * @see #productLow(double, double, double, double, double)
+ */
+ static double productLow2(double x, double y, double xy) {
+ // Verify the input. This must be NaN safe.
+ //assert Double.compare(x * y, xy) == 0
+
+ // If the number is sub-normal, inf or nan there is no round-off.
+ if (isNotNormal(xy)) {
+ // Returns 0.0 for sub-normal xy, otherwise NaN for inf/nan:
+ return xy - xy;
+ }
+
+ // The result xy is finite and normal.
+ // Use Dekker's mul12 algorithm that splits the values into high and low parts.
+ // Dekker's split using multiplication will overflow if the value is within 2^27
+ // of double max value. It can also produce 26-bit approximations that are larger
+ // than the input numbers for the high part causing overflow in hx * hy when
+ // x * y does not overflow. So we must scale down big numbers.
+ // We only have to scale the largest number as we know the product does not overflow
+ // (if one is too big then the other cannot be).
+ // Also scale if the product is close to max value.
+
+ if (Math.abs(x) >= SAFE_UPPER) {
+ return productLowUnscaled(x * DOWN_SCALE, y, xy * DOWN_SCALE) * UP_SCALE;
+ }
+ if (Math.abs(y) >= SAFE_UPPER || Math.abs(xy) >= Double.MAX_VALUE / 4) {
+ return productLowUnscaled(x, y * DOWN_SCALE, xy * DOWN_SCALE) * UP_SCALE;
+ }
+ // No scaling required
+ return productLowUnscaled(x, y, xy);
+ }
+
+ /**
+ * Compute the low part of the double length number {@code (z,zz)} for the exact
+ * product of {@code x} and {@code y} using Dekker's mult12 algorithm. The standard
+ * precision product {@code x*y} must be provided. The numbers {@code x} and {@code y}
+ * are split into high and low parts using Dekker's algorithm.
+ *
+ * <p>This method performs scaling in Dekker's split for large finite numbers to avoid
+ * overflow when generating the high part of the number.
+ *
+ * <p>Warning: Dekker's split can produce high parts that are larger in magnitude than
+ * the input number as the high part is a 26-bit approximation of the number. Thus it is
+ * possible that the standard product {@code x * y} does not overflow but the extended
+ * precision sub-product {@code hx * hy} does overflow. This method should not be
+ * considered safe for all combinations where {@code Double.isFinite(x * y)} is true.
+ * The method is used for benchmarking.
+ *
+ * @param x First factor.
+ * @param y Second factor.
+ * @param xy Product of the factors (x * y).
+ * @return the low part of the product double length number
+ * @see #highPart(double)
+ * @see #productLow(double, double, double, double, double)
+ */
+ static double productLow3(double x, double y, double xy) {
+ // Split the numbers using Dekker's algorithm
+ final double hx = highPart(x);
+ final double lx = x - hx;
+
+ final double hy = highPart(y);
+ final double ly = y - hy;
+
+ return productLow(hx, lx, hy, ly, xy);
+ }
+
+ /**
+ * Compute the low part of the double length number {@code (z,zz)} for the
+ * product of {@code x} and {@code y} using a Dekker's mult12 algorithm. The
+ * standard precision product {@code x*y} must be provided. The numbers
+ * {@code x} and {@code y} are split into high and low parts by zeroing the
+ * lower 27-bits of the mantissa to create the high part. This may lose 1 bit of
+ * precision in the resulting low part computed by subtraction. The intermediate
+ * computations will not overflow as the split results are always smaller in
+ * magnitude than the input numbers.
+ *
+ * <p>The method is used for benchmarking as results may not be exact due to
+ * loss of a bit during splitting of the input factors.
+ *
+ * @param x First factor.
+ * @param y Second factor.
+ * @param xy Product of the factors (x * y).
+ * @return the low part of the product double length number
+ * @see #highPart(double)
+ * @see #productLow(double, double, double, double, double)
+ */
+ static double productLowSplit(double x, double y, double xy) {
+ // Split the numbers using Dekker's algorithm
+ final double hx = highPartSplit(x);
+ final double lx = x - hx;
+
+ final double hy = highPartSplit(y);
+ final double ly = y - hy;
+
+ return productLow(hx, lx, hy, ly, xy);
+ }
+
+ /**
+ * Compute the low part of the double length number {@code (z,zz)} for the exact
+ * product of {@code x} and {@code y} using Dekker's mult12 algorithm. The standard
+ * precision product {@code x*y} must be provided. The numbers {@code x} and {@code y}
+ * are split into high and low parts using Dekker's algorithm.
+ *
+ * <p>Warning: This method does not perform scaling in Dekker's split and large
+ * finite numbers can create NaN results.
+ *
+ * @param x First factor.
+ * @param y Second factor.
+ * @param xy Product of the factors (x * y).
+ * @return the low part of the product double length number
+ * @see #highPartUnscaled(double)
+ * @see #productLow(double, double, double, double, double)
+ */
+ static double productLowUnscaled(double x, double y, double xy) {
+ // Split the numbers using Dekker's algorithm without scaling
+ final double hx = highPartUnscaled(x);
+ final double lx = x - hx;
+
+ final double hy = highPartUnscaled(y);
+ final double ly = y - hy;
+
+ return productLow(hx, lx, hy, ly, xy);
+ }
+
+ /**
+ * Compute the low part of the double length number {@code (z,zz)} for the exact
+ * product of {@code x} and {@code y} using Dekker's mult12 algorithm. The standard
+ * precision product {@code x*y} must be provided. The numbers {@code x} and {@code y}
+ * should already be split into low and high parts.
+ *
+ * <p>Note: This uses the high part of the result {@code (z,zz)} as {@code x * y} and not
+ * {@code hx * hy + hx * ty + tx * hy} as specified in Dekker's original paper.
+ * See Shewchuk (1997) for working examples.
+ *
+ * @param hx High part of first factor.
+ * @param lx Low part of first factor.
+ * @param hy High part of second factor.
+ * @param ly Low part of second factor.
+ * @param xy Product of the factors.
+ * @return <code>lx * ly - (((xy - hx * hy) - lx * hy) - hx * ly)</code>
+ * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+ * Shewchuk (1997) Theorum 18</a>
+ */
+ static double productLow(double hx, double lx, double hy, double ly, double xy) {
+ // Compute the multiply low part:
+ // err1 = xy - hx * hy
+ // err2 = err1 - lx * hy
+ // err3 = err2 - hx * ly
+ // low = lx * ly - err3
+ return lx * ly - (((xy - hx * hy) - lx * hy) - hx * ly);
+ }
+
+ /**
+ * Compute the round-off {@code s} from the sum of two split numbers {@code (x, xx)}
+ * and {@code (y, yy)} using Dekker's add2 algorithm. The values are not required to be
+ * ordered by magnitude as an absolute comparison is made to determine the summation order.
+ * The sum of the high parts {@code r} must be provided.
+ *
+ * <p>The result {@code (r, s)} must be re-balanced to create the split result {@code (z, zz)}:
+ * <pre>
+ * z = r + s
+ * zz = r - z + s
+ * </pre>
+ *
+ * @param x High part of first number.
+ * @param xx Low part of first number.
+ * @param y High part of second number.
+ * @param yy Low part of second number.
+ * @param r Sum of the parts (x + y) = r
+ * @return The round-off from the sum (x + y) = s
+ */
+ static double sumLow(double x, double xx, double y, double yy, double r) {
+ return Math.abs(x) > Math.abs(y) ?
+ x - r + y + yy + xx :
+ y - r + x + xx + yy;
+ }
+
+ /**
+ * Compute the round-off from the sum of two numbers {@code a} and {@code b} using
+ * Dekker's two-sum algorithm. The values are required to be ordered by magnitude
+ * {@code |a| >= |b|}. The standard precision sum must be provided.
+ *
+ * @param a First part of sum.
+ * @param b Second part of sum.
+ * @param sum Sum of the parts (a + b).
+ * @return <code>b - (sum - a)</code>
+ * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+ * Shewchuk (1997) Theorum 6</a>
+ */
+ static double fastTwoSumLow(double a, double b, double sum) {
+ // bVitual = sum - a
+ // b - bVirtual == b round-off
+ return b - (sum - a);
+ }
+
+ /**
+ * Compute the round-off from the sum of two numbers {@code a} and {@code b} using
+ * Knuth's two-sum algorithm. The values are not required to be ordered by magnitude.
+ * The standard precision sum must be provided.
+ *
+ * @param a First part of sum.
+ * @param b Second part of sum.
+ * @param sum Sum of the parts (a + b).
+ * @return <code>(b - (sum - (sum - b))) + (a - (sum - b))</code>
+ * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+ * Shewchuk (1997) Theorum 7</a>
+ */
+ static double twoSumLow(double a, double b, double sum) {
+ final double bVirtual = sum - a;
+ // sum - bVirtual == aVirtual.
+ // a - aVirtual == a round-off
+ // b - bVirtual == b round-off
+ return (a - (sum - bVirtual)) + (b - bVirtual);
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/DoubleSplitPerformance.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/DoubleSplitPerformance.java
new file mode 100644
index 0000000..4aece34
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/DoubleSplitPerformance.java
@@ -0,0 +1,615 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.numbers.examples.jmh.arrays;
+
+import org.apache.commons.rng.UniformRandomProvider;
+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;
+
+import java.util.concurrent.TimeUnit;
+import java.util.function.DoubleBinaryOperator;
+import java.util.function.DoublePredicate;
+import java.util.function.DoubleUnaryOperator;
+
+/**
+ * Executes a benchmark to measure the speed of operations in the {@link LinearCombination} class.
+ * Benchmarks focus on the split of a double value into high and low parts.
+ */
+@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 DoubleSplitPerformance {
+ /** The mask for the sign bit and the mantissa. */
+ private static final long SIGN_MATISSA_MASK = 0x800f_ffff_ffff_ffffL;
+
+ /**
+ * The multiplier used to split the double value into high and low parts. From
+ * Dekker (1971): "The constant should be chosen equal to 2^(p - p/2) + 1,
+ * where p is the number of binary digits in the mantissa". Here p is 53
+ * and the multiplier is {@code 2^27 + 1}.
+ */
+ private static final double MULTIPLIER = 1.34217729E8;
+
+ /** The upper limit above which a number may overflow during the split into a high part.
+ * Assuming the multiplier is above 2^27 and the maximum exponent is 1023 then a safe
+ * limit is a value with an exponent of (1023 - 27) = 2^996. */
+ private static final double SAFE_UPPER = 0x1.0p996;
+
+ /** The scale to use when down-scaling during a split into a high part.
+ * This must be smaller than the inverse of the multiplier and a power of 2 for exact scaling. */
+ private static final double DOWN_SCALE = 0x1.0p-30;
+
+ /** The scale to use when re-scaling during a split into a high part.
+ * This is the inverse of {@link #DOWN_SCALE}. */
+ private static final double UP_SCALE = 0x1.0p30;
+
+ /** The mask to zero the lower 27-bits of a long . */
+ private static final long ZERO_LOWER_27_BITS = 0xffff_ffff_f800_0000L;
+
+ /** Constant to no method. */
+ private static final String NONE = "none";
+
+ /**
+ * The numbers to split.
+ */
+ @State(Scope.Benchmark)
+ public static class Numbers {
+ /** The exponent for small numbers. */
+ private static final long EXP_SMALL = Double.doubleToRawLongBits(1.0);
+ /** The exponent for big numbers. */
+ private static final long EXP_BIG = Double.doubleToRawLongBits(SAFE_UPPER);
+
+ /**
+ * The count of numbers.
+ */
+ @Param({"10000"})
+ private int size;
+
+ /**
+ * The fraction of small numbers.
+ *
+ * <p>Note: The split method may employ multiplications.
+ * Big numbers are edge cases that would cause overflow in multiplications.
+ */
+ @Param({"1", "0.999", "0.99", "0.9"})
+ private double edge;
+
+ /** Numbers. */
+ private double[] a;
+
+ /**
+ * Gets the factors.
+ *
+ * @return Factors.
+ */
+ public double[] getNumbers() {
+ return a;
+ }
+
+ /**
+ * Create the factors.
+ */
+ @Setup
+ public void setup() {
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP);
+ a = new double[size];
+ for (int i = 0; i < size; i++) {
+ long bits = rng.nextLong() & SIGN_MATISSA_MASK;
+ // The exponent will either be small or big
+ if (rng.nextDouble() < edge) {
+ bits |= EXP_SMALL;
+ } else {
+ bits |= EXP_BIG;
+ }
+ a[i] = Double.longBitsToDouble(bits);
+ }
+ }
+ }
+
+ /**
+ * The factors to multiply.
+ */
+ @State(Scope.Benchmark)
+ public static class BiFactors {
+ /** The exponent for small numbers. */
+ private static final long EXP_SMALL = Double.doubleToRawLongBits(1.0);
+
+ /**
+ * The count of products.
+ */
+ @Param({"5000"})
+ private int size;
+
+ /**
+ * The exponent for big numbers.
+ * Two big numbers multiplied together should overflow.
+ *
+ * <p>Note: If this is set below the safe upper for Dekker's split then a method
+ * that computes the split of the two factors independently and then does the multiply
+ * may create split parts that will overflow during multiplication.
+ */
+ @Param({"600", "1000", "1023"})
+ private int exp;
+
+ /**
+ * The fraction of small numbers.
+ *
+ * <p>Note: The split numbers are used in multiplications.
+ * It is unlikely that many numbers will be larger than the upper limit.
+ * These numbers are edge cases that would cause overflow in multiplications if
+ * the other number is anywhere close to the same magnitude.
+ */
+ @Param({"1", "0.95", "0.9"})
+ private double edge;
+
+ /** Factors. */
+ private double[] a;
+
+ /**
+ * Gets the factors to be multiplied as pairs. The array length will be even.
+ *
+ * @return Factors.
+ */
+ public double[] getFactors() {
+ return a;
+ }
+
+ /**
+ * Create the factors.
+ */
+ @Setup
+ public void setup() {
+ // Validate the big exponent
+ final double d = Math.scalb(1.0, exp);
+ assert Double.isInfinite(d * d) : "Product of big numbers does not overflow";
+ final long expBig = Double.doubleToRawLongBits(d);
+
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP);
+ a = new double[size * 2];
+ for (int i = 0; i < a.length; i++) {
+ long bits = rng.nextLong() & SIGN_MATISSA_MASK;
+ // The exponent will either be small or big
+ if (rng.nextDouble() < edge) {
+ bits |= EXP_SMALL;
+ } else {
+ bits |= expBig;
+ }
+ a[i] = Double.longBitsToDouble(bits);
+ }
+ }
+ }
+
+ /**
+ * The numbers to test to determine if they are not normal.
+ */
+ @State(Scope.Benchmark)
+ public static class NonNormalNumbers {
+ /** Non-normal positive numbers. */
+ private static final double[] NON_NORMAL =
+ {Double.POSITIVE_INFINITY, Double.NaN, Double.MIN_NORMAL};
+
+ /**
+ * The count of numbers.
+ */
+ @Param({"10000"})
+ private int size;
+
+ /**
+ * The fraction of non-normal factors.
+ */
+ @Param({"1", "0.999", "0.99", "0.9"})
+ private double edge;
+
+ /** Numbers. */
+ private double[] a;
+
+ /**
+ * Gets the factors.
+ *
+ * @return Factors.
+ */
+ public double[] getFactors() {
+ return a;
+ }
+
+ /**
+ * Create the factors.
+ */
+ @Setup
+ public void setup() {
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP);
+ a = new double[size];
+ for (int i = 0; i < size; i++) {
+ // Value in (-1, 1)
+ double value = rng.nextDouble() * (rng.nextBoolean() ? -1 : 1);
+ // The number will either be small or non-normal
+ if (rng.nextDouble() < edge) {
+ value *= NON_NORMAL[rng.nextInt(NON_NORMAL.length)];
+ }
+ a[i] = value;
+ }
+ }
+ }
+
+ /**
+ * The split method.
+ */
+ @State(Scope.Benchmark)
+ public static class SplitMethod {
+ /**
+ * The name of the method.
+ */
+ @Param({NONE, "dekker", "dekkerAbs", "dekkerRaw", "bits"})
+ private String name;
+
+ /** The function. */
+ private DoubleUnaryOperator fun;
+
+ /**
+ * Gets the function.
+ *
+ * @return the function
+ */
+ public DoubleUnaryOperator getFunction() {
+ return fun;
+ }
+
+ /**
+ * Create the function.
+ */
+ @Setup
+ public void setup() {
+ if (NONE.equals(name)) {
+ fun = a -> a;
+ } else if ("dekker".equals(name)) {
+ fun = DoubleSplitPerformance::splitDekker;
+ } else if ("dekkerAbs".equals(name)) {
+ fun = DoubleSplitPerformance::splitDekkerAbs;
+ } else if ("dekkerRaw".equals(name)) {
+ fun = DoubleSplitPerformance::splitDekkerRaw;
+ } else if ("bits".equals(name)) {
+ fun = DoubleSplitPerformance::splitBits;
+ } else {
+ throw new IllegalStateException("Unknown split method: " + name);
+ }
+ }
+ }
+
+ /**
+ * The method to test for a non-normal number.
+ */
+ @State(Scope.Benchmark)
+ public static class NonNormalMethod {
+ /**
+ * The name of the method.
+ */
+ @Param({NONE, "if", "exponent", "exponent2"})
+ private String name;
+
+ /** The function. */
+ private DoublePredicate fun;
+
+ /**
+ * Gets the function.
+ *
+ * @return the function
+ */
+ public DoublePredicate getFunction() {
+ return fun;
+ }
+
+ /**
+ * Create the function.
+ */
+ @Setup
+ public void setup() {
+ if (NONE.equals(name)) {
+ fun = a -> false;
+ } else if ("if".equals(name)) {
+ fun = DoubleSplitPerformance::isNotNormalIf;
+ } else if ("exponent".equals(name)) {
+ fun = DoubleSplitPerformance::isNotNormalExponent;
+ } else if ("exponent2".equals(name)) {
+ fun = DoubleSplitPerformance::isNotNormalExponent2;
+ } else {
+ throw new IllegalStateException("Unknown is non-normal method: " + name);
+ }
+ }
+ }
+
+ /**
+ * The method to compute the product round-off.
+ */
+ @State(Scope.Benchmark)
+ public static class RoundoffMethod {
+ /**
+ * The name of the method.
+ */
+ @Param({NONE, "multiply", "multiplyUnscaled",
+ "productLow", "productLow1", "productLow2", "productLow3", "productLowSplit",
+ "productLowUnscaled"})
+ private String name;
+
+ /** The function. */
+ private DoubleBinaryOperator fun;
+
+ /**
+ * Gets the function.
+ *
+ * @return the function
+ */
+ public DoubleBinaryOperator getFunction() {
+ return fun;
+ }
+
+ /**
+ * Create the function.
+ */
+ @Setup
+ public void setup() {
+ if (NONE.equals(name)) {
+ // No actually the round-off but x*y - x*y will be optimised away so this
+ // captures the multiply overhead.
+ fun = (x, y) -> x * y;
+ } else if ("multiply".equals(name)) {
+ final DoublePrecision.Quad result = new DoublePrecision.Quad();
+ fun = (x, y) -> {
+ DoublePrecision.multiply(x, y, result);
+ return result.lo;
+ };
+ } else if ("productLow".equals(name)) {
+ fun = (x, y) -> DoublePrecision.productLow(x, y, x * y);
+ } else if ("productLow1".equals(name)) {
+ fun = (x, y) -> DoublePrecision.productLow1(x, y, x * y);
+ } else if ("productLow2".equals(name)) {
+ fun = (x, y) -> DoublePrecision.productLow2(x, y, x * y);
+ } else if ("productLow3".equals(name)) {
+ fun = (x, y) -> DoublePrecision.productLow3(x, y, x * y);
+ } else if ("productLowSplit".equals(name)) {
+ fun = (x, y) -> DoublePrecision.productLowSplit(x, y, x * y);
+ } else if ("multiplyUnscaled".equals(name)) {
+ final DoublePrecision.Quad result = new DoublePrecision.Quad();
+ fun = (x, y) -> {
+ DoublePrecision.multiplyUnscaled(x, y, result);
+ return result.lo;
+ };
+ } else if ("productLowUnscaled".equals(name)) {
+ fun = (x, y) -> DoublePrecision.productLowUnscaled(x, y, x * y);
+ } else {
+ throw new IllegalStateException("Unknown round-off method: " + name);
+ }
+ }
+ }
+
+ /**
+ * Implement Dekker's method to split a value into two parts. Multiplying by (2^s + 1) creates
+ * a big value from which to derive the two split parts.
+ * <pre>
+ * c = (2^s + 1) * a
+ * a_big = c - a
+ * a_hi = c - a_big
+ * a_lo = a - a_hi
+ * a = a_hi + a_lo
+ * </pre>
+ *
+ * <p>The multiplicand allows a p-bit value to be split into
+ * (p-s)-bit value {@code a_hi} and a non-overlapping (s-1)-bit value {@code a_lo}.
+ * Combined they have (p-1) bits of significand but the sign bit of {@code a_lo}
+ * contains a bit of information. The constant is chosen so that s is ceil(p/2) where
+ * the precision p for a double is 53-bits (1-bit of the mantissa is assumed to be
+ * 1 for a non sub-normal number) and s is 27.
+ *
+ * @param value Value.
+ * @return the high part of the value.
+ */
+ private static double splitDekker(double value) {
+ // Avoid overflow
+ if (value >= SAFE_UPPER || value <= -SAFE_UPPER) {
+ // Do scaling.
+ final double x = value * DOWN_SCALE;
+ final double c = MULTIPLIER * x;
+ final double hi = (c - (c - x)) * UP_SCALE;
+ if (Double.isInfinite(hi)) {
+ // Number is too large.
+ // This occurs if value is infinite or close to Double.MAX_VALUE.
+ // Note that Dekker's split creates an approximating 26-bit number which may
+ // have an exponent 1 greater than the input value. This will overflow if the
+ // exponent is already +1023. Revert to the raw upper 26 bits of the 53-bit
+ // mantissa (including the assumed leading 1 bit). This conversion will result in
+ // the low part being a 27-bit significand and the potential loss of bits during
+ // addition and multiplication. (Contrast to the Dekker split which creates two
+ // 26-bit numbers with a bit of information moved to the sign of low.)
+ // The conversion will maintain Infinite in the high part where the resulting
+ // low part a_lo = a - a_hi = inf - inf = NaN.
+ return Double.longBitsToDouble(Double.doubleToRawLongBits(value) & ZERO_LOWER_27_BITS);
+ }
+ return hi;
+ }
+ // normal conversion
+ final double c = MULTIPLIER * value;
+ return c - (c - value);
+ }
+
+ /**
+ * Implement Dekker's method to split a value into two parts.
+ * This is as per {@link #splitDekker(double)} but uses a {@link Math#abs(double)} in the
+ * condition statement.
+ *
+ * @param value Value.
+ * @return the high part of the value.
+ */
+ private static double splitDekkerAbs(double value) {
+ if (Math.abs(value) >= SAFE_UPPER) {
+ final double x = value * DOWN_SCALE;
+ final double c = MULTIPLIER * x;
+ final double hi = (c - (c - x)) * UP_SCALE;
+ if (Double.isInfinite(hi)) {
+ return Double.longBitsToDouble(Double.doubleToRawLongBits(value) & ZERO_LOWER_27_BITS);
+ }
+ return hi;
+ }
+ final double c = MULTIPLIER * value;
+ return c - (c - value);
+ }
+
+ /**
+ * Implement Dekker's method to split a value into two parts.
+ * This is as per {@link #splitDekker(double)} but has no overflow protection.
+ *
+ * @param value Value.
+ * @return the high part of the value.
+ */
+ private static double splitDekkerRaw(double value) {
+ final double c = MULTIPLIER * value;
+ return c - (c - value);
+ }
+
+ /**
+ * Implement a split using the upper and lower raw bits from the value.
+ *
+ * @param value Value.
+ * @return the high part of the value.
+ */
+ private static double splitBits(double value) {
+ return Double.longBitsToDouble(Double.doubleToRawLongBits(value) & ZERO_LOWER_27_BITS);
+ }
+
+ /**
+ * Checks if the number is not normal. This is functionally equivalent to:
+ * <pre>
+ * </pre>
+ *
+ * @param a The value.
+ * @return true if the value is not normal
+ */
+ private static boolean isNotNormalIf(double a) {
+ final double abs = Math.abs(a);
+ return abs <= Double.MIN_NORMAL || !(abs <= Double.MAX_VALUE);
+ }
+
+ /**
+ * Checks if the number is not normal.
+ *
+ * @param a The value.
+ * @return true if the value is not normal
+ */
+ private static boolean isNotNormalExponent(double a) {
+ // Sub-normal numbers have a biased exponent of 0.
+ // Inf/NaN numbers have a biased exponent of 2047.
+ // Catch both cases by extracting the raw exponent, subtracting 1
+ // and make unsigned. 0 will underflow to a large value.
+ final int baisedExponent = ((int) (Double.doubleToRawLongBits(a) >>> 52)) & 0x7ff;
+ return ((baisedExponent - 1) & 0xffff) >= 2046;
+ }
+
+ /**
+ * Checks if the number is not normal.
+ *
+ * @param a The value.
+ * @return true if the value is not normal
+ */
+ private static boolean isNotNormalExponent2(double a) {
+ // Sub-normal numbers have a biased exponent of 0.
+ // Inf/NaN numbers have a biased exponent of 2047.
+ // Catch both cases by extracting the raw exponent, subtracting 1
+ // and compare unsigned (so 0 underflows to a large value).
+ final int baisedExponent = ((int) (Double.doubleToRawLongBits(a) >>> 52)) & 0x7ff;
+ // Adding int min value is equal to compare unsigned
+ return baisedExponent + Integer.MIN_VALUE - 1 >= 2046 + Integer.MIN_VALUE;
+ }
+
+ // Benchmark methods.
+
+ /**
+ * Benchmark extracting the high part of the split number.
+ *
+ * @param numbers Factors.
+ * @param bh Data sink.
+ * @param method Split method
+ */
+ @Benchmark
+ public void high(Numbers numbers, Blackhole bh, SplitMethod method) {
+ final DoubleUnaryOperator fun = method.getFunction();
+ final double[] a = numbers.getNumbers();
+ for (int i = 0; i < a.length; i++) {
+ bh.consume(fun.applyAsDouble(a[i]));
+ }
+ }
+
+ /**
+ * Benchmark extracting the low part of the split number.
+ *
+ * @param numbers Factors.
+ * @param bh Data sink.
+ * @param method Split method.
+ */
+ @Benchmark
+ public void low(Numbers numbers, Blackhole bh, SplitMethod method) {
+ final DoubleUnaryOperator fun = method.getFunction();
+ final double[] a = numbers.getNumbers();
+ for (int i = 0; i < a.length; i++) {
+ bh.consume(a[i] - fun.applyAsDouble(a[i]));
+ }
+ }
+
+ /**
+ * Benchmark testing if a number is non-normal.
+ *
+ * @param numbers Factors.
+ * @param bh Data sink.
+ * @param method Split method.
+ */
+ @Benchmark
+ public void nonNormal(NonNormalNumbers numbers, Blackhole bh, NonNormalMethod method) {
+ final DoublePredicate fun = method.getFunction();
+ final double[] a = numbers.getFactors();
+ for (int i = 0; i < a.length; i++) {
+ bh.consume(fun.test(a[i]));
+ }
+ }
+
+ /**
+ * Benchmark extracting the round-off from the product of two numbers.
+ *
+ * @param factors Factors.
+ * @param bh Data sink.
+ * @param method Round-off method.
+ */
+ @Benchmark
+ public void productLow(BiFactors factors, Blackhole bh, RoundoffMethod method) {
+ final DoubleBinaryOperator fun = method.getFunction();
+ final double[] a = factors.getFactors();
+ for (int i = 0; i < a.length; i += 2) {
+ bh.consume(fun.applyAsDouble(a[i], a[i + 1]));
+ }
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombination.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombination.java
new file mode 100644
index 0000000..b7bf902
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombination.java
@@ -0,0 +1,98 @@
+/*
+ * 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;
+
+/**
+ * Computes linear combinations as the the sum of the products of two sequences of numbers
+ * <code>a<sub>i</sub> b<sub>i</sub></code>.
+ *
+ * @see <a href="https://en.wikipedia.org/wiki/Dot_product">Dot product</a>
+ */
+public interface LinearCombination {
+ /**
+ * Compute the sum of the products of two sequences of 2 factors.
+ */
+ @FunctionalInterface
+ interface TwoD {
+ /**
+ * Compute the sum of the products of two sequences of factors.
+ *
+ * @param a1 First factor of the first term.
+ * @param b1 Second factor of the first term.
+ * @param a2 First factor of the second term.
+ * @param b2 Second factor of the second term.
+ * @return \( a_1 b_1 + a_2 b_2 \)
+ */
+ double value(double a1, double b1, double a2, double b2);
+ }
+
+ /**
+ * Compute the sum of the products of two sequences of 3 factors.
+ */
+ @FunctionalInterface
+ interface ThreeD {
+ /**
+ * Compute the sum of the products of two sequences of factors.
+ *
+ * @param a1 First factor of the first term.
+ * @param b1 Second factor of the first term.
+ * @param a2 First factor of the second term.
+ * @param b2 Second factor of the second term.
+ * @param a3 First factor of the third term.
+ * @param b3 Second factor of the third term.
+ * @return \( a_1 b_1 + a_2 b_2 + a_3 b_3 \)
+ */
+ double value(double a1, double b1, double a2, double b2, double a3, double b3);
+ }
+
+ /**
+ * Compute the sum of the products of two sequences of 4 factors.
+ */
+ @FunctionalInterface
+ interface FourD {
+ /**
+ * Compute the sum of the products of two sequences of factors.
+ *
+ * @param a1 First factor of the first term.
+ * @param b1 Second factor of the first term.
+ * @param a2 First factor of the second term.
+ * @param b2 Second factor of the second term.
+ * @param a3 First factor of the third term.
+ * @param b3 Second factor of the third term.
+ * @param a4 First factor of the fourth term.
+ * @param b4 Second factor of the fourth term.
+ * @return \( a_1 b_1 + a_2 b_2 + a_3 b_3 + a_4 b_4 \)
+ */
+ double value(double a1, double b1, double a2, double b2, double a3, double b3, double a4, double b4);
+ }
+
+ /**
+ * Compute the sum of the products of two sequences of {@code n} factors.
+ */
+ @FunctionalInterface
+ interface ND {
+ /**
+ * Compute the sum of the products of two sequences of factors.
+ *
+ * @param a Factors.
+ * @param b Factors.
+ * @return \( \sum_i a_i b_i \).
+ * @throws IllegalArgumentException if the sizes of the arrays are different.
+ */
+ double value(double[] a, double[] b);
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinationPerformance.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinationPerformance.java
new file mode 100644
index 0000000..e620a87
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinationPerformance.java
@@ -0,0 +1,371 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.numbers.examples.jmh.arrays;
+
+import org.apache.commons.numbers.arrays.LinearCombination;
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.FourD;
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.ND;
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.ThreeD;
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.TwoD;
+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;
+
+import java.math.MathContext;
+import java.util.Arrays;
+import java.util.concurrent.TimeUnit;
+import java.util.function.IntFunction;
+
+/**
+ * Executes a benchmark to measure the speed of operations in the {@link LinearCombination} 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 LinearCombinationPerformance {
+ /**
+ * The seed to use to create the factors.
+ * Using a fixed seed ensures the same factors are created for the variable
+ * length arrays as for the small fixed size arrays.
+ */
+ private static final long SEED = System.currentTimeMillis();
+
+ /**
+ * The factors to multiply.
+ */
+ @State(Scope.Benchmark)
+ public static class Factors {
+ /**
+ * The condition number of the generated data.
+ */
+ @Param({"1e20"})
+ private double c;
+
+ /**
+ * The number of factors.
+ */
+ @Param({"1000"})
+ private int size;
+
+ /** Factors a. */
+ private double[][] a;
+
+ /** Factors b. */
+ private double[][] b;
+
+ /**
+ * Gets the length of the array of factors.
+ * This exists to be overridden by factors of a specific length.
+ * The default is to create factors of length 4 for use in the inlined
+ * scalar product methods.
+ *
+ * @return the length
+ */
+ public int getLength() {
+ return 4;
+ }
+
+ /**
+ * Gets the number of scalar products to compute.
+ *
+ * @return the size
+ */
+ public int getSize() {
+ return size;
+ }
+
+ /**
+ * Gets the a factors.
+ *
+ * @param index the index
+ * @return Factors b.
+ */
+ public double[] getA(int index) {
+ return a[index];
+ }
+
+ /**
+ * Gets the b factors.
+ *
+ * @param index the index
+ * @return Factors b.
+ */
+ public double[] getB(int index) {
+ return b[index];
+ }
+
+ /**
+ * Create the factors.
+ */
+ @Setup
+ public void setup() {
+ final UniformRandomProvider rng =
+ RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP, SEED);
+ // Use the ill conditioned data generation method.
+ // This requires an array of at least 6.
+ final int n = Math.max(6, getLength());
+ final double[] x = new double[n];
+ final double[] y = new double[n];
+ a = new double[size][];
+ b = new double[size][];
+ // Limit precision to allow large array lengths to be generated.
+ final MathContext mathContext = new MathContext(100);
+ for (int i = 0; i < size; i++) {
+ LinearCombinationUtils.genDot(c, rng, x, y, null, mathContext);
+ a[i] = Arrays.copyOf(x, getLength());
+ b[i] = Arrays.copyOf(y, getLength());
+ }
+ }
+ }
+
+ /**
+ * The factors to multiply of a specific length.
+ */
+ @State(Scope.Benchmark)
+ public static class LengthFactors extends Factors {
+ /**
+ * The length of each factors array.
+ */
+ @Param({"2", "3", "4", "8", "16", "32", "64"})
+ private int length;
+
+ /** {@inheritDoc} */
+ @Override
+ public int getLength() {
+ return length;
+ }
+ }
+
+ /**
+ * The {@link LinearCombination} implementation.
+ */
+ @State(Scope.Benchmark)
+ public static class Calculator {
+ /**
+ * The implementation name.
+ */
+ @Param({"standard",
+ "current",
+ "dekker",
+ "dot2s",
+ "dot2", "dot3", "dot4", "dot5", "dot6", "dot7",
+ "exact",
+ "extended", "extended2", "extended_exact", "extended_exact2",
+ // Cached working double[] array.
+ // Only faster when 'length' is >16. Below this the array
+ // is small enough to be allocated locally
+ // (Search for Thread Local Allocation Buffer (TLAB))
+ "dot3c", "extendedc"})
+ private String name;
+
+ /** The 2D implementation. */
+ private TwoD twod;
+ /** The 3D implementation. */
+ private ThreeD threed;
+ /** The 4D implementation. */
+ private FourD fourd;
+ /** The ND implementation. */
+ private ND nd;
+
+ /**
+ * @return the 2D implementation
+ */
+ public TwoD getTwoD() {
+ return twod;
+ }
+
+ /**
+ * @return the 3D implementation
+ */
+ public ThreeD getThreeD() {
+ return threed;
+ }
+
+ /**
+ * @return the 4D implementation
+ */
+ public FourD getFourD() {
+ return fourd;
+ }
+
+ /**
+ * @return the ND implementation
+ */
+ public ND getND() {
+ return nd;
+ }
+
+ /**
+ * Setup the implementation.
+ */
+ @Setup
+ public void setup() {
+ if ("current".endsWith(name)) {
+ twod = LinearCombination::value;
+ threed = LinearCombination::value;
+ fourd = LinearCombination::value;
+ nd = LinearCombination::value;
+ return;
+ }
+ // All implementations below are expected to implement all the interfaces.
+ if ("standard".endsWith(name)) {
+ nd = LinearCombinations.StandardPrecision.INSTANCE;
+ } else if ("dekker".equals(name)) {
+ nd = LinearCombinations.Dekker.INSTANCE;
+ } else if ("dot2s".equals(name)) {
+ nd = LinearCombinations.Dot2s.INSTANCE;
+ } else if ("dot2".equals(name)) {
+ nd = new LinearCombinations.DotK(2);
+ } else if ("dot3".equals(name)) {
+ nd = LinearCombinations.DotK.DOT_3;
+ } else if ("dot4".equals(name)) {
+ nd = LinearCombinations.DotK.DOT_4;
+ } else if ("dot5".equals(name)) {
+ nd = LinearCombinations.DotK.DOT_5;
+ } else if ("dot6".equals(name)) {
+ nd = LinearCombinations.DotK.DOT_6;
+ } else if ("dot7".equals(name)) {
+ nd = LinearCombinations.DotK.DOT_7;
+ } else if ("exact".equals(name)) {
+ nd = LinearCombinations.Exact.INSTANCE;
+ } else if ("extended".equals(name)) {
+ nd = LinearCombinations.ExtendedPrecision.INSTANCE;
+ } else if ("extended2".equals(name)) {
+ nd = LinearCombinations.ExtendedPrecision.DOUBLE;
+ } else if ("extended_exact".equals(name)) {
+ nd = LinearCombinations.ExtendedPrecision.EXACT;
+ } else if ("extended_exact2".equals(name)) {
+ nd = LinearCombinations.ExtendedPrecision.EXACT2;
+ } else if ("dot3c".equals(name)) {
+ nd = new LinearCombinations.DotK(3, new CachedArrayFactory());
+ } else if ("extendedc".equals(name)) {
+ nd = LinearCombinations.ExtendedPrecision.of(
+ LinearCombinations.ExtendedPrecision.Summation.STANDARD, new CachedArrayFactory());
+ } else {
+ throw new IllegalStateException("Unknown implementation: " + name);
+ }
+ // Possible class-cast exception for partial implementations...
+ twod = (TwoD) nd;
+ threed = (ThreeD) nd;
+ fourd = (FourD) nd;
+ }
+ }
+
+ /**
+ * Create or return a cached array.
+ */
+ static final class CachedArrayFactory implements IntFunction<double[]> {
+ /** An empty double array. */
+ private static final double[] EMPTY = new double[0];
+
+ /** The cached array. */
+ private double[] array = EMPTY;
+
+ @Override
+ public double[] apply(int value) {
+ double[] a = array;
+ if (a.length < value) {
+ array = a = new double[value];
+ }
+ return a;
+ }
+ }
+
+ /**
+ * Compute the 2D scalar product for all the factors.
+ *
+ * @param factors Factors.
+ * @param bh Data sink.
+ * @param calc Scalar product calculator.
+ */
+ @Benchmark
+ public void twoD(Factors factors, Blackhole bh, Calculator calc) {
+ final TwoD fun = calc.getTwoD();
+ for (int i = 0; i < factors.getSize(); i++) {
+ final double[] a = factors.getA(i);
+ final double[] b = factors.getB(i);
+ bh.consume(fun.value(a[0], b[0], a[1], b[1]));
+ }
+ }
+
+ /**
+ * Compute the 3D scalar product for all the factors.
+ *
+ * @param factors Factors.
+ * @param bh Data sink.
+ * @param calc Scalar product calculator.
+ */
+ @Benchmark
+ public void threeD(Factors factors, Blackhole bh, Calculator calc) {
+ final ThreeD fun = calc.getThreeD();
+ for (int i = 0; i < factors.getSize(); i++) {
+ final double[] a = factors.getA(i);
+ final double[] b = factors.getB(i);
+ bh.consume(fun.value(a[0], b[0], a[1], b[1], a[2], b[2]));
+ }
+ }
+
+ /**
+ * Compute the 4D scalar product for all the factors.
+ *
+ * @param factors Factors.
+ * @param bh Data sink.
+ * @param calc Scalar product calculator.
+ */
+ @Benchmark
+ public void fourD(Factors factors, Blackhole bh, Calculator calc) {
+ final FourD fun = calc.getFourD();
+ for (int i = 0; i < factors.getSize(); i++) {
+ final double[] a = factors.getA(i);
+ final double[] b = factors.getB(i);
+ bh.consume(fun.value(a[0], b[0], a[1], b[1], a[2], b[2], a[3], b[3]));
+ }
+ }
+
+ /**
+ * Compute the ND scalar product for all the factors.
+ *
+ * @param factors Factors.
+ * @param bh Data sink.
+ * @param calc Scalar product calculator.
+ */
+ @Benchmark
+ public void nD(LengthFactors factors, Blackhole bh, Calculator calc) {
+ final ND fun = calc.getND();
+ for (int i = 0; i < factors.getSize(); i++) {
+ // These should be pre-computed to the correct length
+ final double[] a = factors.getA(i);
+ final double[] b = factors.getB(i);
+ bh.consume(fun.value(a, b));
+ }
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinationUtils.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinationUtils.java
new file mode 100644
index 0000000..a13a47d
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinationUtils.java
@@ -0,0 +1,211 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.numbers.examples.jmh.arrays;
+
+import org.apache.commons.rng.UniformRandomProvider;
+
+import java.math.BigDecimal;
+import java.math.MathContext;
+
+/**
+ * Utility class to create data for linear combinations.
+ */
+final class LinearCombinationUtils {
+ /** ln(2). */
+ private static final double LN_2 = Math.log(2);
+
+ /** No construction. */
+ private LinearCombinationUtils() {}
+ /**
+ * Generates ill conditioned dot products.
+ * See {@link #genDot(double, UniformRandomProvider, double[], double[], double[], MathContext)}.
+ *
+ * <p>The exact dot product is computed using {@link MathContext#UNLIMITED}. This
+ * will not scale to large lengths.
+ *
+ * @param c anticipated condition number of x’*y
+ * @param rng source of randomness
+ * @param x output array vector (length at least 6)
+ * @param y output array vector (length at least 6)
+ * @param computeC if not null compute the condition number and place at index 0
+ * @return the exact dot product
+ * @throws IllegalArgumentException If the vector length is below 6
+ */
+ static double genDot(double c, UniformRandomProvider rng,
+ double[] x, double[] y, double[] computeC) {
+ return genDot(c, rng, x, y, computeC, MathContext.UNLIMITED);
+ }
+
+ /**
+ * Generates ill conditioned dot products. The length of the vectors should be
+ * {@code n>=6}.
+ *
+ * <p>The condition number is defined as the inverse of the cosine of the angle between
+ * the two vectors:
+ *
+ * <pre>
+ * C = 1 / cos angle(x, y)
+ * = ||x'|| ||y|| / |x'*y|
+ * </pre>
+ * <p>where |x'*y| is the absolute of the dot product and ||.|| is the 2-norm (i.e. the
+ * Euclidean length of each vector).
+ *
+ * <p>A high condition number means that small perturbations in the values result in
+ * a large change in the final result. This occurs when the cosine of the angle approaches
+ * 0 and the vectors are close to orthogonal.
+ *
+ * <p>Computation of the actual dot product requires an exact method to avoid cancellation
+ * floating-point errors. BigDecimal is used to compute the exact result to the accuracy
+ * defined by the {@link MathContext}. The dot product result is created in the range [-1, 1].
+ * The actual condition number can be obtained by passing an array {@code computeC}. The
+ * routine has been tested with anticipated condition numbers up to 1e300 to generate data
+ * where the standard precision dot product will not overflow.
+ *
+ * <p>Uses the GenDot algorithm 6.1 from <a
+ * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
+ * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump, and
+ * Shin'ichi Oishi published in <em>SIAM J. Sci. Comput</em>.
+ *
+ * <p>Note: Ogita et al state "in general the actual condition number is a
+ * little larger than the anticipated one with quite some variation". This
+ * implementation is computing a condition number approximately equal to the
+ * requested for the parameters that have been tested.
+ *
+ * @param c anticipated condition number of x’*y
+ * @param rng source of randomness
+ * @param x output array vector (length at least 6)
+ * @param y output array vector (length at least 6)
+ * @param computeC if not null compute the condition number and place at index 0
+ * @param context the context to use for rounding the exact sum
+ * @return the exact dot product
+ * @throws IllegalArgumentException If the vector length is below 6
+ * @see <a href="https://en.wikipedia.org/wiki/Condition_number">Condition number</a>
+ * @see <a href="https://en.wikipedia.org/wiki/Dot_product">Dot product</a>
+ */
+ static double genDot(double c, UniformRandomProvider rng,
+ double[] x, double[] y, double[] computeC, MathContext context) {
+ // Initialisation
+ if (x.length < 6) {
+ throw new IllegalArgumentException("Incorrect usage: n < 6: " + x.length);
+ }
+ final int n = x.length;
+ final int n2 = n / 2;
+
+ // log2(c)
+ final double b = Math.log(c) / LN_2;
+ final double b2 = b / 2;
+ // e vector of exponents between 0 and b/2
+ int[] e = new int[n2];
+ // make sure exponents b/2 and 0 actually occur in e
+ e[0] = (int) Math.round(b2) + 1;
+ for (int i = 1; i < n2 - 1; i++) {
+ e[i] = (int) Math.round(rng.nextDouble() * b2);
+ }
+ // e[end] = 0;
+
+ // Generate first half vectors.
+ // Maintain the exact dot product for use later
+ BigDecimal exact = BigDecimal.ZERO;
+ for (int i = 0; i < n2; i++) {
+ x[i] = Math.scalb(m1p1(rng), e[i]);
+ y[i] = Math.scalb(m1p1(rng), e[i]);
+ exact = exact.add(new BigDecimal(x[i]).multiply(new BigDecimal(y[i])), context);
+ }
+
+ // for i=n2+1:n and v=1:i,
+ // generate x(i), y(i) such that (*) x(v)’*y(v) ~ 2^e(i-n2)
+ // i.e. the dot product up to position i is a value that will increasingly approach 0
+ e = new int[n - n2];
+ // exponents for second half as a linear vector of exponents from b/2 to 0
+ // linspace(b/2, 0, n-n2)
+ for (int i = 0; i < e.length - 1; i++) {
+ e[i] = (int) Math.round(b2 * (e.length - i - 1) / (e.length - 1));
+ }
+
+ for (int i = n2; i < x.length; i++) {
+ // x(i) random with generated exponent
+ x[i] = Math.scalb(m1p1(rng), e[i - n2]);
+ // y(i) according to (*).
+ // sum(i) = xi * yi + sum(i-1)
+ // yi = (sum(i) - sum(i-1)) / xi
+ // Here the new sum(i) is a random number with the exponent gradually
+ // reducing to e[end] = 0 so the final result is in [-1, 1].
+ y[i] = (Math.scalb(m1p1(rng), e[i - n2]) - exact.doubleValue()) / x[i];
+ // Maintain the exact dot product
+ exact = exact.add(new BigDecimal(x[i]).multiply(new BigDecimal(y[i])), context);
+ }
+
+ // Shuffle x and y. Do a parallel Fisher-Yates shuffle.
+ for (int i = n; i > 1; i--) {
+ final int j = rng.nextInt(i);
+ swap(x, i - 1, j);
+ swap(y, i - 1, j);
+ }
+
+ // Ogita el at.
+ // Compute condition number:
+ // d = DotExact(x’,y); % the true dot product rounded to nearest
+ // C = 2*(abs(x’)*abs(y))/abs(d); % the actual condition number
+
+ // Compare to this:
+ // https://math.stackexchange.com/questions/3147927/condition-number-of-dot-product-of-vectors
+ // Compute the inverse of the cosine between the two vectors.
+ // ||y^t|| ||x|| / | y^t x |
+ // This value is similar to that computed by Ogita et al which effectively
+ // is using the magnitude of the largest component to approximate the vector
+ // length. This works as the vectors are created with matched magnitudes in their
+ // components. Here we compute the actual vector lengths ||y^t|| and ||x||.
+ final double d = exact.doubleValue();
+ if (computeC != null) {
+ // Sum should not overflow as elements are bounded to an exponent half the size
+ // of the input condition number.
+ double s1 = 0;
+ double s2 = 0;
+ for (int i = 0; i < n; i++) {
+ s1 += x[i] * x[i];
+ s2 += y[i] * y[i];
+ }
+ computeC[0] = Math.sqrt(s1) * Math.sqrt(s2) / Math.abs(d);
+ }
+ return d;
+ }
+
+ /**
+ * Create a double in the range [-1, 1).
+ *
+ * @param rng source of randomness
+ * @return the double
+ */
+ private static double m1p1(UniformRandomProvider rng) {
+ // Create in the range [0, 1) then randomly subtract 1.
+ // This samples the 2^54 dyadic rationals in the range.
+ return rng.nextDouble() - rng.nextInt(1);
+ }
+
+ /**
+ * Swaps the two specified elements in the array.
+ *
+ * @param array Array.
+ * @param i First index.
+ * @param j Second index.
+ */
+ private static void swap(double[] array, int i, int j) {
+ final double tmp = array[i];
+ array[i] = array[j];
+ array[j] = tmp;
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinations.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinations.java
new file mode 100644
index 0000000..ad81380
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinations.java
@@ -0,0 +1,2007 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.numbers.examples.jmh.arrays;
+
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.FourD;
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.ND;
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.ThreeD;
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.TwoD;
+
+import java.math.BigDecimal;
+import java.util.function.IntFunction;
+
+/**
+ * Provides implementations to computes linear combinations as the the sum of
+ * the products of two sequences of numbers
+ * <code>a<sub>i</sub> b<sub>i</sub></code>.
+ *
+ * @see LinearCombination
+ */
+public final class LinearCombinations {
+
+ /** No public constructor. */
+ private LinearCombinations() {}
+
+ /**
+ * Base class to compute a linear combination with high accuracy.
+ * Contains common code for computing short combinations and computing
+ * the standard precision sum of products.
+ */
+ public abstract static class BaseLinearCombination implements LinearCombination.ND {
+ /** {@inheritDoc} */
+ @Override
+ public double value(double[] a, double[] b) {
+ if (a.length != b.length) {
+ throw new IllegalArgumentException("Dimension mismatch: " + a.length + " != " + b.length);
+ }
+
+ final int len = a.length;
+
+ if (len == 0) {
+ return 0;
+ }
+ if (len == 1) {
+ // Revert to scalar multiplication.
+ return a[0] * b[0];
+ }
+
+ return computeValue(a, b);
+ }
+
+ /**
+ * Compute the sum of the products of two sequences of factors with high accuracy.
+ * The input arrays will have a length of at least 2; the lengths will
+ * be the same.
+ *
+ * @param a Factors.
+ * @param b Factors.
+ * @return \( \sum_i a_i b_i \).
+ */
+ protected abstract double computeValue(double[] a, double[] b);
+
+ /**
+ * Compute the sum of the products of two sequences of factors.
+ * This is a standard precision summation for the IEEE754 result.
+ *
+ * @param a Factors.
+ * @param b Factors.
+ * @return \( \sum_i a_i b_i \).
+ */
+ static double standardDotProduct(double[] a, double[] b) {
+ double result = 0;
+ for (int i = 0; i < a.length; ++i) {
+ result += a[i] * b[i];
+ }
+ return result;
+ }
+ }
+
+ /**
+ * Computes linear combinations using standard precision multiplication and summation.
+ *
+ * <p>This class is used for a baseline and does not provide high precision results.
+ */
+ static final class StandardPrecision implements TwoD, ThreeD, FourD, ND {
+ /** An instance. */
+ static final StandardPrecision INSTANCE = new StandardPrecision();
+
+ /** Private constructor. */
+ private StandardPrecision() {}
+
+ @Override
+ public double value(double[] a, double[] b) {
+ double sum = 0;
+ for (int i = 0; i < a.length; i++) {
+ sum += a[i] * b[i];
+ }
+ return sum;
+ }
+
+ @Override
+ public double value(double a1, double b1, double a2, double b2) {
+ return a1 * b1 + a2 * b2;
+ }
+
+ @Override
+ public double value(double a1, double b1, double a2, double b2, double a3, double b3) {
+ return a1 * b1 + a2 * b2 + a3 * b3;
+ }
+
+ @Override
+ public double value(double a1, double b1, double a2, double b2, double a3, double b3, double a4, double b4) {
+ return a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4;
+ }
+ }
+
+ /**
+ * Computes linear combinations using the double-length multiplication and summation
+ * algorithms of Dekker.
+ *
+ * @see <a href="https://doi.org/10.1007/BF01397083">
+ * Dekker (1971) A floating-point technique for extending the available precision</a>
+ */
+ public static final class Dekker extends BaseLinearCombination implements TwoD, ThreeD, FourD {
+ /** An instance. */
+ public static final Dekker INSTANCE = new Dekker();
+
+ /** Private constructor. */
+ private Dekker() {}
+
+ @Override
+ protected double computeValue(double[] a, double[] b) {
+ // Dekker's dot product method:
+ // Dekker's multiply and then Dekker's add2 to sum in double length precision.
+
+ double x = a[0] * b[0];
+ double xx = DoublePrecision.productLow(a[0], b[0], x);
+
+ // Remaining split products added to the current sum and round-off stored
+ for (int i = 1; i < a.length; i++) {
+ final double y = a[i] * b[i];
+ final double yy = DoublePrecision.productLow(a[i], b[i], y);
+
+ // sum this to the previous number using Dekker's add2 algorithm
+ final double r = x + y;
+ final double s = DoublePrecision.sumLow(x, xx, y, yy, r);
+
+ x = r + s;
+ // final split product does not require the round-off but included here for brevity
+ xx = r - x + s;
+ }
+
+ if (!Double.isFinite(x)) {
+ // Either we have split infinite numbers or some coefficients were NaNs,
+ // just rely on the naive implementation and let IEEE754 handle this
+ return standardDotProduct(a, b);
+ }
+ return x;
+ }
+
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2) {
+ final double x = a1 * b1;
+ final double xx = DoublePrecision.productLow(a1, b1, x);
+ final double y = a2 * b2;
+ final double yy = DoublePrecision.productLow(a2, b2, y);
+ double r = x + y;
+ final double s = DoublePrecision.sumLow(x, xx, y, yy, r);
+ r = r + s;
+
+ if (!Double.isFinite(r)) {
+ // Either we have split infinite numbers or some coefficients were NaNs,
+ // just rely on the naive implementation and let IEEE754 handle this
+ return x + y;
+ }
+ return r;
+ }
+
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2,
+ double a3, double b3) {
+ double x = a1 * b1;
+ double xx = DoublePrecision.productLow(a1, b1, x);
+ double y = a2 * b2;
+ double yy = DoublePrecision.productLow(a2, b2, y);
+ double r = x + y;
+ double s = DoublePrecision.sumLow(x, xx, y, yy, r);
+ x = r + s;
+ xx = r - x + s;
+ y = a3 * b3;
+ yy = DoublePrecision.productLow(a3, b3, y);
+ r = x + y;
+ s = DoublePrecision.sumLow(x, xx, y, yy, r);
+ x = r + s;
+
+ if (!Double.isFinite(x)) {
+ // Either we have split infinite numbers or some coefficients were NaNs,
+ // just rely on the naive implementation and let IEEE754 handle this
+ return a1 * b1 + a2 * b2 + y;
+ }
+ return x;
+ }
+
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2,
+ double a3, double b3,
+ double a4, double b4) {
+ double x = a1 * b1;
+ double xx = DoublePrecision.productLow(a1, b1, x);
+ double y = a2 * b2;
+ double yy = DoublePrecision.productLow(a2, b2, y);
+ double r = x + y;
+ double s = DoublePrecision.sumLow(x, xx, y, yy, r);
+ x = r + s;
+ xx = r - x + s;
+ y = a3 * b3;
+ yy = DoublePrecision.productLow(a3, b3, y);
+ r = x + y;
+ s = DoublePrecision.sumLow(x, xx, y, yy, r);
+ x = r + s;
+ xx = r - x + s;
+ y = a4 * b4;
+ yy = DoublePrecision.productLow(a4, b4, y);
+ r = x + y;
+ s = DoublePrecision.sumLow(x, xx, y, yy, r);
+ x = r + s;
+
+ if (!Double.isFinite(x)) {
+ // Either we have split infinite numbers or some coefficients were NaNs,
+ // just rely on the naive implementation and let IEEE754 handle this
+ return a1 * b1 + a2 * b2 + a3 * b3 + y;
+ }
+ return x;
+ }
+ }
+
+ /**
+ * Computes linear combinations accurately using the DotK algorithm of Ogita et al
+ * for K-fold precision of the sum.
+ *
+ * <p>It is based on the 2005 paper
+ * <a href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
+ * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump,
+ * and Shin'ichi Oishi published in <em>SIAM J. Sci. Comput</em>.
+ *
+ * <p>Note: It is possible to use this class to compute 2-fold precision. In this case
+ * the round-off parts of the dot-product are stored in an array and summed. The
+ * {@link Dot2s} class provides an alternative faster implementation that sums the
+ * round-off parts during processing to avoid array allocation overhead. The results will
+ * not be identical due to a different order for the summation of the round-off parts.
+ *
+ * <p>The number of operations will scale with {@code k}. When computing the small linear
+ * combinations {@link TwoD}, {@link ThreeD} and {@link FourD} using a {@code k} value too
+ * large will result in more operations than the number used by {@link ExtendedPrecision}.
+ * For a small {@code n}-dimension linear combinations the value {@code k = n+1} will have the
+ * same number of operations as {@link ExtendedPrecision} and users should switch to the
+ * {@link ExtendedPrecision} class. For example 3D combinations should use {@code k} up to
+ * 3 or else switch to using {@link ExtendedPrecision}. This rule does not apply to the
+ * {@link org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.ND ND}
+ * implementations as {@link ExtendedPrecision} eliminates redundant operations
+ * on zeros. In this case equivalent performance will be observed when {@code k <= n+1}.
+ * Choice of implementation must consider performance and accuracy on real-world data.
+ */
+ public static final class DotK extends BaseLinearCombination implements TwoD, ThreeD, FourD {
+ /** An instance computing 3-fold precision. */
+ public static final DotK DOT_3 = new DotK(3);
+ /** An instance computing 4-fold precision. */
+ public static final DotK DOT_4 = new DotK(4);
+ /** An instance computing 5-fold precision. */
+ public static final DotK DOT_5 = new DotK(5);
+ /** An instance computing 6-fold precision. */
+ public static final DotK DOT_6 = new DotK(6);
+ /** An instance computing 7-fold precision. */
+ public static final DotK DOT_7 = new DotK(7);
+
+ /** The k-fold precision to compute. */
+ private final int k;
+ /** The factory to create a working array. */
+ private final IntFunction<double[]> arrayFactory;
+
+ /**
+ * Instantiates a new dot K.
+ *
+ * @param k K-fold precision.
+ */
+ public DotK(int k) {
+ this(k, double[]::new);
+ }
+
+ /**
+ * Instantiates a new dot K.
+ *
+ * @param k K-fold precision.
+ * @param arrayFactory The factory to create a working array.
+ */
+ DotK(int k, IntFunction<double[]> arrayFactory) {
+ this.k = k;
+ this.arrayFactory = arrayFactory;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected double computeValue(double[] a, double[] b) {
+ // Implement dotK (Algorithm 5.10) from Ogita et al (2005).
+ // Store all round-off parts.
+ // Round-off parts of each product are r[0 to (n-1)].
+ // Round-off parts of each sum are r[n to (2n-2)].
+ // The standard precision scalar product is term p which becomes r[2n-1].
+ final int len = a.length;
+ final double[] r = arrayFactory.apply(len * 2);
+
+ // p is the standard scalar product sum initialised with the first product
+ double p = a[0] * b[0];
+ r[0] = DoublePrecision.productLow(a[0], b[0], p);
+
+ // Remaining split products added to the current sum and round-off stored
+ for (int i = 1; i < len; i++) {
+ final double h = a[i] * b[i];
+ r[i] = DoublePrecision.productLow(a[i], b[i], h);
+
+ final double x = p + h;
+ r[i + len - 1] = DoublePrecision.twoSumLow(p, h, x);
+ p = x;
+ }
+
+ // Sum the round-off with the standard sum as the final component.
+ // Here the value passed to sumK is (K-1) for K-fold precision of the sum (dotK).
+ // Increasing K may be of benefit for longer arrays.
+ // The iteration is an error-free transform and higher K never loses precision.
+ r[r.length - 1] = p;
+ return getSum(p, sumK(r, k - 1));
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2) {
+ // Round-off parts of each product are r[0-1].
+ // Round-off parts of each sum are r[2].
+ // Working variables p/q (new/old sum).
+ // The standard precision scalar product is stored in s.
+
+ double p = a1 * b1;
+ double r0 = DoublePrecision.productLow(a1, b1, p);
+ double q = a2 * b2;
+ double r1 = DoublePrecision.productLow(a2, b2, q);
+ final double s = p + q;
+ double r2 = DoublePrecision.twoSumLow(p, q, s);
+ double r3 = s;
+
+ // In-line k-2 rounds of vector sum for k-fold precision
+ for (int i = 2; i < k; i++) {
+ q = r1 + r0;
+ r0 = DoublePrecision.twoSumLow(r1, r0, q);
+ p = r2 + q;
+ r1 = DoublePrecision.twoSumLow(r2, q, p);
+ q = r3 + p;
+ r2 = DoublePrecision.twoSumLow(r3, p, q);
+ r3 = q;
+ }
+
+ // Final summation
+ return getSum(s, r0 + r1 + r2 + r3);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2,
+ double a3, double b3) {
+ // Round-off parts of each product are r[0-2].
+ // Round-off parts of each sum are r[3-4].
+ // Working variables p/q (new/old sum) and h (current product high part).
+ // The standard precision scalar product is stored in s.
+
+ double p = a1 * b1;
+ double r0 = DoublePrecision.productLow(a1, b1, p);
+ double h = a2 * b2;
+ double r1 = DoublePrecision.productLow(a2, b2, h);
+ double q = p + h;
+ double r3 = DoublePrecision.twoSumLow(p, h, q);
+ h = a3 * b3;
+ double r2 = DoublePrecision.productLow(a3, b3, h);
+ final double s = q + h;
+ double r4 = DoublePrecision.twoSumLow(q, h, s);
+ double r5 = s;
+
+ // In-line k-2 rounds of vector sum for k-fold precision
+ for (int i = 2; i < k; i++) {
+ q = r1 + r0;
+ r0 = DoublePrecision.twoSumLow(r1, r0, q);
+ p = r2 + q;
+ r1 = DoublePrecision.twoSumLow(r2, q, p);
+ q = r3 + p;
+ r2 = DoublePrecision.twoSumLow(r3, p, q);
+ p = r4 + q;
+ r3 = DoublePrecision.twoSumLow(r4, q, p);
+ q = r5 + p;
+ r4 = DoublePrecision.twoSumLow(r5, p, q);
+ r5 = q;
+ }
+
+ // Final summation
+ return getSum(s, r0 + r1 + r2 + r3 + r4 + r5);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2,
+ double a3, double b3,
+ double a4, double b4) {
+ // Round-off parts of each product are r[0-3].
+ // Round-off parts of each sum are r[4-6].
+ // Working variables p/q (new/old sum) and h (current product high part).
+ // The standard precision scalar product is stored in s.
+
+ double p = a1 * b1;
+ double r0 = DoublePrecision.productLow(a1, b1, p);
+ double h = a2 * b2;
+ double r1 = DoublePrecision.productLow(a2, b2, h);
+ double q = p + h;
+ double r4 = DoublePrecision.twoSumLow(p, h, q);
+ h = a3 * b3;
+ double r2 = DoublePrecision.productLow(a3, b3, h);
+ p = q + h;
+ double r5 = DoublePrecision.twoSumLow(q, h, p);
+ h = a4 * b4;
+ double r3 = DoublePrecision.productLow(a4, b4, h);
+ final double s = p + h;
+ double r6 = DoublePrecision.twoSumLow(p, h, s);
+ double r7 = s;
+
+ // In-line k-2 rounds of vector sum for k-fold precision
+ for (int i = 2; i < k; i++) {
+ q = r1 + r0;
+ r0 = DoublePrecision.twoSumLow(r1, r0, q);
+ p = r2 + q;
+ r1 = DoublePrecision.twoSumLow(r2, q, p);
+ q = r3 + p;
+ r2 = DoublePrecision.twoSumLow(r3, p, q);
+ p = r4 + q;
+ r3 = DoublePrecision.twoSumLow(r4, q, p);
+ q = r5 + p;
+ r4 = DoublePrecision.twoSumLow(r5, p, q);
+ p = r6 + q;
+ r5 = DoublePrecision.twoSumLow(r6, q, p);
+ q = r7 + p;
+ r6 = DoublePrecision.twoSumLow(r7, p, q);
+ r7 = q;
+ }
+
+ // Final summation
+ return getSum(s, r0 + r1 + r2 + r3 + r4 + r5 + r6 + r7);
+ }
+
+ /**
+ * Sum to K-fold precision.
+ *
+ * @param p Data to sum.
+ * @param km1 The precision (k-1).
+ * @return the sum
+ */
+ private static double sumK(double[] p, int km1) {
+ // (k-1)=1 will skip the vector transformation and sum in standard precision.
+ for (int i = 1; i < km1; i++) {
+ vectorSum(p);
+ }
+ double sum = 0;
+ for (final double pi : p) {
+ sum += pi;
+ }
+ return sum;
+ }
+
+ /**
+ * Error free vector transformation for summation.
+ *
+ * @param p Data.
+ */
+ private static void vectorSum(double[] p) {
+ for (int i = 1; i < p.length; i++) {
+ final double x = p[i] + p[i - 1];
+ p[i - 1] = DoublePrecision.twoSumLow(p[i], p[i - 1], x);
+ p[i] = x;
+ }
+ }
+ }
+
+ /**
+ * Computes linear combinations accurately using the Dot2s algorithm of Ogita et al
+ * for 2-fold precision of the sum.
+ *
+ * <p>It is based on the 2005 paper
+ * <a href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
+ * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump,
+ * and Shin'ichi Oishi published in <em>SIAM J. Sci. Comput</em>.
+ *
+ * <p>This is faster than using {@link DotK} with a {@code k} of 2. The results will
+ * not be identical due to a different summation order of the round-off parts.
+ */
+ public static final class Dot2s extends BaseLinearCombination implements TwoD, ThreeD, FourD {
+ /** An instance computing 2-fold precision. */
+ public static final Dot2s INSTANCE = new Dot2s();
+
+ /** Private constructor. */
+ private Dot2s() {}
+
+ /** {@inheritDoc} */
+ @Override
+ protected double computeValue(double[] a, double[] b) {
+ // Implement dot2s (Algorithm 5.4) from Ogita et al (2005).
+ final int len = a.length;
+
+ // p is the standard scalar product sum.
+ // s is the sum of round-off parts.
+ double p = a[0] * b[0];
+ double s = DoublePrecision.productLow(a[0], b[0], p);
+
+ // Remaining split products added to the current sum and round-off sum.
+ for (int i = 1; i < len; i++) {
+ final double h = a[i] * b[i];
+ final double r = DoublePrecision.productLow(a[i], b[i], h);
+
+ final double x = p + h;
+ // s_i = s_(i-1) + (q_i + r_i)
+ s += DoublePrecision.twoSumLow(p, h, x) + r;
+ p = x;
+ }
+
+ return getSum(p, p + s);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2) {
+ // p/pn are the standard scalar product old/new sum.
+ // s is the sum of round-off parts.
+ final double p = a1 * b1;
+ double s = DoublePrecision.productLow(a1, b1, p);
+ final double h = a2 * b2;
+ final double r = DoublePrecision.productLow(a2, b2, h);
+ final double pn = p + h;
+ s += DoublePrecision.twoSumLow(p, h, pn) + r;
+
+ // Final summation
+ return getSum(pn, pn + s);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2,
+ double a3, double b3) {
+ // Sum round-off parts in s: s_i = s_(i-1) + (q_i + r_i)
+ // The standard precision scalar product is stored in p_n.
+ final double p = a1 * b1;
+ double s = DoublePrecision.productLow(a1, b1, p);
+ double h = a2 * b2;
+ double r = DoublePrecision.productLow(a2, b2, h);
+ final double q = p + h;
+ s += r + DoublePrecision.twoSumLow(p, h, q);
+ h = a3 * b3;
+ r = DoublePrecision.productLow(a3, b3, h);
+ final double pn = q + h;
+ s += r + DoublePrecision.twoSumLow(q, h, pn);
+
+ // Final summation
+ return getSum(pn, pn + s);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2,
+ double a3, double b3,
+ double a4, double b4) {
+ // p/q are the standard scalar product old/new sum (alternating).
+ // s is the sum of round-off parts.
+ // pn is the final scalar product sum.
+ double p = a1 * b1;
+ double s = DoublePrecision.productLow(a1, b1, p);
+ double h = a2 * b2;
+ double r = DoublePrecision.productLow(a2, b2, h);
+ final double q = p + h;
+ s += DoublePrecision.twoSumLow(p, h, q) + r;
+ h = a3 * b3;
+ r = DoublePrecision.productLow(a3, b3, h);
+ p = q + h;
+ s += DoublePrecision.twoSumLow(q, h, p) + r;
+ h = a4 * b4;
+ r = DoublePrecision.productLow(a4, b4, h);
+ final double pn = p + h;
+ s += DoublePrecision.twoSumLow(p, h, pn) + r;
+
+ // Final summation
+ return getSum(pn, pn + s);
+ }
+ }
+
+ /**
+ * Computes linear combinations exactly using BigDecimal.
+ * This computation may be prohibitively slow on large combination sums.
+ */
+ public static final class Exact extends BaseLinearCombination implements TwoD, ThreeD, FourD {
+ /** An instance. */
+ public static final Exact INSTANCE = new Exact();
+
+ /** Private constructor. */
+ private Exact() {}
+
+ @Override
+ protected double computeValue(double[] a, double[] b) {
+ // BigDecimal cannot handle inf/nan so check the arguments and return the IEEE754 result
+ if (!areFinite(a) || !areFinite(b)) {
+ return standardDotProduct(a, b);
+ }
+ BigDecimal sum = multiply(a[0], b[0]);
+ for (int i = 1; i < a.length; i++) {
+ sum = sum.add(multiply(a[i], b[i]));
+ }
+ return sum.doubleValue();
+ }
+
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2) {
+ // BigDecimal cannot handle inf/nan so check the arguments and return the IEEE754 result
+ if (!areFinite(a1, b1, a2, b2)) {
+ return a1 * b1 + a2 * b2;
+ }
+ return multiply(a1, b1)
+ .add(multiply(a2, b2))
+ .doubleValue();
+ }
+
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2,
+ double a3, double b3) {
+ // BigDecimal cannot handle inf/nan so check the arguments and return the IEEE754 result
+ if (!areFinite(a1, b1, a2, b2, a3, b3)) {
+ return a1 * b1 + a2 * b2 + a3 * b3;
+ }
+ return multiply(a1, b1)
+ .add(multiply(a2, b2))
+ .add(multiply(a3, b3))
+ .doubleValue();
+ }
+
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2,
+ double a3, double b3,
+ double a4, double b4) {
+ // BigDecimal cannot handle inf/nan so check the arguments and return the IEEE754 result
+ if (!areFinite(a1, b1, a2, b2, a3, b3, a4, b4)) {
+ return a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4;
+ }
+ return multiply(a1, b1)
+ .add(multiply(a2, b2))
+ .add(multiply(a3, b3))
+ .add(multiply(a4, b4))
+ .doubleValue();
+ }
+
+ /**
+ * Test that all the values are finite.
+ *
+ * @param values the values
+ * @return true if finite
+ */
+ private static boolean areFinite(double... values) {
+ for (final double value : values) {
+ if (!Double.isFinite(value)) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ /**
+ * Multiply the factors to a BigDecimal.
+ *
+ * @param a First factor.
+ * @param b Second factor.
+ * @return the product
+ */
+ private static BigDecimal multiply(double a, double b) {
+ return new BigDecimal(a).multiply(new BigDecimal(b));
+ }
+ }
+
+ /**
+ * Computes linear combinations accurately using extended precision representations of
+ * floating point numbers.
+ *
+ * <p>It is based on the paper by
+ * <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+ * Shewchuk (1997): Arbitrary Precision Floating-Point Arithmetic</a>.
+ */
+ public static final class ExtendedPrecision extends BaseLinearCombination implements TwoD, ThreeD, FourD {
+ /*
+ * Note:
+ *
+ * An expansion representation of a number is a series of non-overlapping floating-point
+ * values where the most significant bit of each value is less than the least significant
+ * bit of the next value. The summation of the expansion is exact (without round-off error)
+ * and is equal to the original number. The largest magnitude value in the expansion is
+ * an approximation of the number.
+ *
+ * Expansions are created from normal numbers by multiplications or additions that
+ * represent the result exactly with extended precision. Addition of expansions f and g
+ * of length m and n creates an expansion of length (m+n). Some parts of the expansion
+ * may be zero and can be eliminated. The size of the expansion is constrained by the
+ * maximum and minimum exponents required to represent the original results of
+ * multiplication and addition. If all products in the linear combination have
+ * approximately the same magnitude the expansion size will be small and the sum
+ * is efficient. If the products have a very large range of magnitudes then the expansion
+ * will require a large size and the sum is computationally more expensive.
+ */
+
+ /**
+ * An instance computing the final summation of the extended precision value
+ * using standard precision. The sum will be within 1 ULP of the exact result.
+ */
+ public static final ExtendedPrecision INSTANCE = new ExtendedPrecision(ApproximationSum.INSTANCE);
+
+ /**
+ * An instance computing the final summation of the extended precision value
+ * using two-fold precision. The sum will be within 1 ULP of the exact result.
+ */
+ public static final ExtendedPrecision DOUBLE = new ExtendedPrecision(TwoSum.INSTANCE);
+
+ /**
+ * An instance computing the final summation of the extended precision value exactly.
+ */
+ public static final ExtendedPrecision EXACT = new ExtendedPrecision(ExactSum.INSTANCE);
+
+ /**
+ * An instance computing the final summation of the extended precision value exactly.
+ */
+ public static final ExtendedPrecision EXACT2 = new ExtendedPrecision(ExactSum2.INSTANCE);
+
+ /**
+ * Define methods to sum an expansion.
+ */
+ private interface ExpansionSum {
+ /**
+ * Sum the expansion of size m. It can be assumed that the size is non-zero.
+ *
+ * @param e Expansion.
+ * @param m Size.
+ * @return the sum
+ */
+ double sum(double[] e, int m);
+
+ /**
+ * Sum the expansion; low parts have smaller magnitudes.
+ *
+ * @param e0 Part 0.
+ * @param e1 Part 1.
+ * @param e2 Part 2.
+ * @param e3 Part 3.
+ * @return the sum
+ */
+ double sum(double e0, double e1, double e2, double e3);
+
+ /**
+ * Sum the expansion; low parts have smaller magnitudes.
+ *
+ * @param e0 Part 0.
+ * @param e1 Part 1.
+ * @param e2 Part 2.
+ * @param e3 Part 3.
+ * @param e4 Part 4.
+ * @param e5 Part 5.
+ * @return the sum
+ */
+ double sum(double e0, double e1, double e2, double e3, double e4, double e5);
+
+ /**
+ * Sum the expansion; low parts have smaller magnitudes.
+ *
+ * @param e0 Part 0.
+ * @param e1 Part 1.
+ * @param e2 Part 2.
+ * @param e3 Part 3.
+ * @param e4 Part 4.
+ * @param e5 Part 5.
+ * @param e6 Part 6.
+ * @param e7 Part 7.
+ * @return the sum
+ */
+ double sum(double e0, double e1, double e2, double e3, double e4, double e5, double e6, double e7);
+ }
+
+ /**
+ * Shewchuk's APPROXIMATION algorithm sums the parts in order of increasing magnitude.
+ * The sum is within 1 ulp of the true expansion value.
+ */
+ private static final class ApproximationSum implements ExpansionSum {
+ /** An instance. */
+ static final ApproximationSum INSTANCE = new ApproximationSum();
+
+ /** No construction. */
+ private ApproximationSum() {}
+
+ @Override
+ public double sum(double[] e, int m) {
+ double sum = 0;
+ for (int i = 0; i < m; i++) {
+ sum += e[i];
+ }
+ return sum;
+ }
+
+ @Override
+ public double sum(double e0, double e1, double e2, double e3) {
+ return e0 + e1 + e2 + e3;
+ }
+
+ @Override
+ public double sum(double e0, double e1, double e2, double e3, double e4, double e5) {
+ return e0 + e1 + e2 + e3 + e4 + e5;
+ }
+
+ @Override
+ public double sum(double e0, double e1, double e2, double e3, double e4, double e5, double e6, double e7) {
+ return e0 + e1 + e2 + e3 + e4 + e5 + e6 + e7;
+ }
+ }
+
+ /**
+ * Perform a two-sum through the expansion.
+ * The single pass with two-sum ensures that the final term e_m is a good approximation
+ * for e: |e - e_m| < ulp(e_m); and the sum of the parts to
+ * e_(m-1) is within 1 ULP of the round-off ulp(|e - e_m|).
+ */
+ private static final class TwoSum implements ExpansionSum {
+ /** An instance. */
+ static final TwoSum INSTANCE = new TwoSum();
+
+ /** No construction. */
+ private TwoSum() {}
+
+ @Override
+ public double sum(double[] e, int m) {
+ double q = e[0];
+ double qq = 0.0;
+ for (int i = 1; i < m; i++) {
+ final double p = q + e[i];
+ qq += DoublePrecision.twoSumLow(q, e[i], p);
+ q = p;
+ }
+ return q + qq;
+ }
+
+ @Override
+ public double sum(double e0, double e1, double e2, double e3) {
+ double q = e0 + e1;
+ double qq = DoublePrecision.fastTwoSumLow(e0, e1, q);
+ final double p = q + e2;
+ qq += DoublePrecision.twoSumLow(q, e2, p);
+ q = p + e3;
+ qq += DoublePrecision.twoSumLow(p, e3, q);
+ return q + qq;
+ }
+
+ @Override
+ public double sum(double e0, double e1, double e2, double e3, double e4, double e5) {
+ double q = e0 + e1;
+ double qq = DoublePrecision.fastTwoSumLow(e0, e1, q);
+ double p = q + e2;
+ qq += DoublePrecision.twoSumLow(q, e2, p);
+ q = p + e3;
+ qq += DoublePrecision.twoSumLow(p, e3, q);
+ p = q + e4;
+ qq += DoublePrecision.twoSumLow(q, e4, p);
+ q = p + e5;
+ qq += DoublePrecision.twoSumLow(p, e5, q);
+ return q + qq;
+ }
+
+ @Override
+ public double sum(double e0, double e1, double e2, double e3, double e4, double e5, double e6, double e7) {
+ double q = e0 + e1;
+ double qq = DoublePrecision.fastTwoSumLow(e0, e1, q);
+ double p = q + e2;
+ qq += DoublePrecision.twoSumLow(q, e2, p);
+ q = p + e3;
+ qq += DoublePrecision.twoSumLow(p, e3, q);
+ p = q + e4;
+ qq += DoublePrecision.twoSumLow(q, e4, p);
+ q = p + e5;
+ qq += DoublePrecision.twoSumLow(p, e5, q);
+ p = q + e6;
+ qq += DoublePrecision.twoSumLow(q, e6, p);
+ q = p + e7;
+ qq += DoublePrecision.twoSumLow(p, e7, q);
+ return q + qq;
+ }
+ }
+
+ /**
+ * Shewchuk's COMPRESS(e) algorithm uses fast-two-sum with a double pass
+ * (big->small then small->big) and zero elimination for up to 6m additions (zero
+ * elimination reduces the size of the second pass). This ensure e_m is within 1
+ * ULP of the expansion value: |e - e_m| < ulp(e_m); and the sum to e_(m-1) is
+ * within 1 ULP of the round-off ulp(|e - e_m|). Summation of the expansion should
+ * then have low error.
+ *
+ * <p>Here we modify the second pass to set a sticky bit on the carried sum
+ * that allows the rounding of the next addition to compute round-to-nearest,
+ * ties-to-even rounding. This simulates the floating point addition of a+b into
+ * an extended register where some extra digits are held beyond the precision of
+ * the result for use in rounding. The final digit is the sticky bit, a bit wise OR
+ * of all the bits from the exact result that are discarded. This works if the value
+ * with the sticky bit is added to a value with a higher exponent. The use of the
+ * COMPRESS algorithm ensures the next addition is to a number larger in magnitude.
+ *
+ * <p>Note: The final addition of the sticky carry must be to a non-zero value.
+ * Otherwise the sticky bit will be left in the result potentially causing a 1
+ * ULP error.
+ *
+ * <p>Details of the sticky bit can be found in:
+ * <blockquote>
+ * Coonen, J.T., "An Implementation Guide to a Proposed Standard for Floating Point
+ * Arithmetic", Computer, Vol. 13, No. 1, Jan. 1980, pp 68-79.
+ * </blockquote>
+ */
+ private static class ExactSum implements ExpansionSum {
+ /** An instance. */
+ static final ExactSum INSTANCE = new ExactSum();
+
+ /** Package-private construction to allow class extension. */
+ ExactSum() {}
+
+ @Override
+ public double sum(double[] e, int size) {
+ if (size == 1) {
+ return e[0];
+ }
+
+ // First traversal (from big to small)
+ // Shewchuk uses (Q,q) for (big,small) from fast-two-sum and carries Q.
+ // Here use (q,qq); p is used for intermediates.
+ final int m = size - 1;
+ double q = e[m];
+ int bottom = m;
+ for (int i = m - 1; i >= 0; i--) {
+ final double p = q + e[i];
+ final double qq = DoublePrecision.fastTwoSumLow(q, e[i], p);
+ if (qq != 0) {
+ // Store larger component in e and carry the smaller
+ e[bottom] = p;
+ bottom--;
+ q = qq;
+ } else {
+ // Compression.
+ q = p;
+ }
+ }
+ // Redundant store
+ // e[bottom] = q
+
+ if (bottom == m) {
+ // Complete compression to size 1
+ return q;
+ }
+
+ // Second traversal (from small to big)
+ // Here we do not store the round-off parts back into the expansion
+ // but summarise them into the carry using a sticky bit.
+ for (int i = bottom + 1; i < m; i++) {
+ q = fastSumWithStickyBit(e[i], q);
+ }
+ // Final sum. No requirement to compute round-off.
+ return e[m] + q;
+ }
+
+ @Override
+ public double sum(double e0, double e1, double e2, double e3) {
+ // compress(e) into g.
+ double g3 = 0.0;
+ double g2 = 0.0;
+ double g1 = 0.0;
+
+ // First traversal (from big to small).
+ double q = e3 + e2;
+ double qq = DoublePrecision.fastTwoSumLow(e3, e2, q);
+ if (qq != 0) {
+ g3 = q;
+ q = qq;
+ }
+ double p = q + e1;
+ qq = DoublePrecision.fastTwoSumLow(q, e1, p);
+ if (qq != 0) {
+ g2 = p;
+ p = qq;
+ }
+ q = p + e0;
+ qq = DoublePrecision.fastTwoSumLow(p, e0, q);
+ if (qq != 0) {
+ g1 = p;
+ q = qq;
+ }
+ // g0 = q
+
+ // Second traversal (from small to big)
+ // Here we do not store the round-off parts back into the expansion
+ // but summarise them into the carry using a sticky bit sum.
+ return stickySum(q, g1, g2, g3);
+ }
+
+ @Override
+ public double sum(double e0, double e1, double e2, double e3, double e4, double e5) {
+ // compress(e) into g.
+ double g5 = 0.0;
+ double g4 = 0.0;
+ double g3 = 0.0;
+ double g2 = 0.0;
+ double g1 = 0.0;
+
+ // First traversal (from big to small).
+ double q = e5 + e4;
+ double qq = DoublePrecision.fastTwoSumLow(e5, e4, q);
+ if (qq != 0) {
+ g5 = q;
+ q = qq;
+ }
+ double p = q + e3;
+ qq = DoublePrecision.fastTwoSumLow(q, e3, p);
+ if (qq != 0) {
+ g4 = p;
+ p = qq;
+ }
+ q = p + e2;
+ qq = DoublePrecision.fastTwoSumLow(p, e2, q);
+ if (qq != 0) {
+ g3 = q;
+ q = qq;
+ }
+ p = q + e1;
+ qq = DoublePrecision.fastTwoSumLow(q, e1, p);
+ if (qq != 0) {
+ g2 = p;
+ p = qq;
+ }
+ q = p + e0;
+ qq = DoublePrecision.fastTwoSumLow(p, e0, q);
+ if (qq != 0) {
+ g1 = q;
+ q = qq;
+ }
+ // g0 = q
+
+ // Second traversal (from small to big)
+ // Here we do not store the round-off parts back into the expansion
+ // but summarise them into the carry using a sticky bit sum.
+ return stickySum(q, g1, g2, g3, g4, g5);
+ }
+
+ @Override
+ public double sum(double e0, double e1, double e2, double e3, double e4, double e5, double e6, double e7) {
+ // compress(e) into g.
+ double g7 = 0.0;
+ double g6 = 0.0;
+ double g5 = 0.0;
+ double g4 = 0.0;
+ double g3 = 0.0;
+ double g2 = 0.0;
+ double g1 = 0.0;
+
+ // First traversal (from big to small).
+ // If round-off is non-zero it is carried down and the high part deposited in g.
+ // If round-off is zero compression has occurred. The high part is carried and
+ // gi remains zero. This does not eliminate zeros from g but spurious zeros are
+ // left to avoid use of array indexing. These are ignored from the sum operations.
+ double q = e7 + e6;
+ double qq = DoublePrecision.fastTwoSumLow(e7, e6, q);
+ if (qq != 0) {
+ g7 = q;
+ q = qq;
+ }
+ double p = q + e5;
+ qq = DoublePrecision.fastTwoSumLow(q, e5, p);
+ if (qq != 0) {
+ g6 = p;
+ p = qq;
+ }
+ q = p + e4;
+ qq = DoublePrecision.fastTwoSumLow(p, e4, q);
+ if (qq != 0) {
+ g5 = q;
+ q = qq;
+ }
+ p = q + e3;
+ qq = DoublePrecision.fastTwoSumLow(q, e3, p);
+ if (qq != 0) {
+ g4 = p;
+ p = qq;
+ }
+ q = p + e2;
+ qq = DoublePrecision.fastTwoSumLow(p, e2, q);
+ if (qq != 0) {
+ g3 = q;
+ q = qq;
+ }
+ p = q + e1;
+ qq = DoublePrecision.fastTwoSumLow(q, e1, p);
+ if (qq != 0) {
+ g2 = p;
+ p = qq;
+ }
+ q = p + e0;
+ qq = DoublePrecision.fastTwoSumLow(p, e0, q);
+ if (qq != 0) {
+ g1 = q;
+ q = qq;
+ }
+ // g0 = q
+
+ // Second traversal (from small to big)
+ // Here we do not store the round-off parts back into the expansion
+ // but summarise them into the carry using a sticky bit sum.
+ return stickySum(q, g1, g2, g3, g4, g5, g6, g7);
+ }
+
+ /**
+ * Compute the summation of the parts using a fast-two-sum. The remainder is
+ * carried into the next part using a sticky bit for the final correctly
+ * rounded result.
+ *
+ * @param g0 Part 0
+ * @param g1 Part 1
+ * @param g2 Part 2
+ * @param g3 Part 3
+ * @param g4 Part 4
+ * @param g5 Part 5
+ * @param g6 Part 6
+ * @param g7 Part 7
+ * @return the sum
+ */
+ private static double stickySum(double g0, double g1, double g2, double g3,
+ double g4, double g5, double g6, double g7) {
+ if (g7 == 0) {
+ return stickySum(g0, g1, g2, g3, g4, g5, g6);
+ }
+ double q = fastSumWithStickyBit(g1, g0);
+ q = fastSumWithStickyBit(g2, q);
+ q = fastSumWithStickyBit(g3, q);
+ q = fastSumWithStickyBit(g4, q);
+ q = fastSumWithStickyBit(g5, q);
+ return g7 + fastSumWithStickyBit(g6, q);
+ }
+
+ /**
+ * Compute the summation of the parts using a fast-two-sum. The remainder is
+ * carried into the next part using a sticky bit for the final correctly
+ * rounded result.
+ *
+ * @param g0 Part 0
+ * @param g1 Part 1
+ * @param g2 Part 2
+ * @param g3 Part 3
+ * @param g4 Part 4
+ * @param g5 Part 5
+ * @param g6 Part 6
+ * @return the sum
+ */
+ private static double stickySum(double g0, double g1, double g2, double g3,
+ double g4, double g5, double g6) {
+ if (g6 == 0) {
+ return stickySum(g0, g1, g2, g3, g4, g5);
+ }
+ double q = fastSumWithStickyBit(g1, g0);
+ q = fastSumWithStickyBit(g2, q);
+ q = fastSumWithStickyBit(g3, q);
+ q = fastSumWithStickyBit(g4, q);
+ return g6 + fastSumWithStickyBit(g5, q);
+ }
+
+ /**
+ * Compute the summation of the parts using a fast-two-sum. The remainder is
+ * carried into the next part using a sticky bit for the final correctly
+ * rounded result.
+ *
+ * @param g0 Part 0
+ * @param g1 Part 1
+ * @param g2 Part 2
+ * @param g3 Part 3
+ * @param g4 Part 4
+ * @param g5 Part 5
+ * @return the sum
+ */
+ private static double stickySum(double g0, double g1, double g2, double g3,
+ double g4, double g5) {
+ if (g5 == 0) {
+ return stickySum(g0, g1, g2, g3, g4);
+ }
+ double q = fastSumWithStickyBit(g1, g0);
+ q = fastSumWithStickyBit(g2, q);
+ q = fastSumWithStickyBit(g3, q);
+ return g5 + fastSumWithStickyBit(g4, q);
+ }
+
+ /**
+ * Compute the summation of the parts using a fast-two-sum. The remainder is
+ * carried into the next part using a sticky bit for the final correctly
+ * rounded result.
+ *
+ * @param g0 Part 0
+ * @param g1 Part 1
+ * @param g2 Part 2
+ * @param g3 Part 3
+ * @param g4 Part 4
+ * @return the sum
+ */
+ private static double stickySum(double g0, double g1, double g2, double g3,
+ double g4) {
+ if (g4 == 0) {
+ return stickySum(g0, g1, g2, g3);
+ }
+ double q = fastSumWithStickyBit(g1, g0);
+ q = fastSumWithStickyBit(g2, q);
+ return g4 + fastSumWithStickyBit(g3, q);
+ }
+
+ /**
+ * Compute the summation of the parts using a fast-two-sum. The remainder is
+ * carried into the next part using a sticky bit for the final correctly
+ * rounded result.
+ *
+ * @param g0 Part 0
+ * @param g1 Part 1
+ * @param g2 Part 2
+ * @param g3 Part 3
+ * @return the sum
+ */
+ private static double stickySum(double g0, double g1, double g2, double g3) {
+ if (g3 == 0) {
+ return stickySum(g0, g1, g2);
+ }
+ final double q = fastSumWithStickyBit(g1, g0);
+ return g3 + fastSumWithStickyBit(g2, q);
+ }
+
+ /**
+ * Compute the summation of the parts using a fast-two-sum. The remainder is
+ * carried into the next part using a sticky bit for the final correctly
+ * rounded result.
+ *
+ * @param g0 Part 0
+ * @param g1 Part 1
+ * @param g2 Part 2
+ * @return the sum
+ */
+ private static double stickySum(double g0, double g1, double g2) {
+ if (g2 == 0) {
+ // No sticky bit needed
+ return g1 + g0;
+ }
+ return g2 + fastSumWithStickyBit(g1, g0);
+ }
+
+ /**
+ * Compute the sum of two numbers {@code a} and {@code b} using
+ * Dekker's two-sum algorithm. The values are required to be ordered by magnitude
+ * {@code |a| >= |b|}. The result is adjusted to set the lowest bit as a sticky
+ * bit that summarises the magnitude of the round-off that were lost. The
+ * result is not the correctly rounded result; it is intended the result is to
+ * be used in an addition with a value with a greater magnitude exponent. This
+ * addition will have exact round-to-nearest, ties-to-even rounding taking account
+ * of bits lots in the previous sum.
+ *
+ * <p>Details of the sticky bit can be found in:
+ * <blockquote>
+ * Coonen, J.T., "An Implementation Guide to a Proposed Standard for Floating Point
+ * Arithmetic", Computer, Vol. 13, No. 1, Jan. 1980, pp 68-79.
+ * </blockquote>
+ *
+ * @param a First part of sum.
+ * @param b Second part of sum.
+ * @return <code>b - (sum - a)</code>
+ * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+ * Shewchuk (1997) Theorum 6</a>
+ */
+ static double fastSumWithStickyBit(double a, double b) {
+ double sum = a + b;
+ // bVitual = sum - a
+ // b - bVirtual == b round-off
+ final double r = b - (sum - a);
+
+ if (r != 0) {
+ // Bits will be lost.
+ // In floating-point arithmetic the sticky bit is the bit-wise OR
+ // of the rest of the binary bits that cannot be stored in the
+ // preliminary representation of the result:
+ //
+ // sgn | exp | V | N | .............. | L | G | R | S
+ //
+ // sgn : sign bit
+ // exp : exponent
+ // V : overflow for significand field
+ // N and L : most and least significant bits of the result
+ // G and R : the two bits beyond
+ // S : sticky bit, bitwise OR of all bits thereafter
+ //
+ // The sticky bit is a flag indicating if there is more magnitude beyond
+ // the last bits. Here the round-off is signed so we have to consider the
+ // sign of the sum and round-off together and either add the sticky or
+ // remove it. The final bit is thus used to push up the next addition using
+ // the sum to a higher value, or down to a lower value, when tie breaking for
+ // the correct round-to-nearest, ties-to-even result.
+ long hi = Double.doubleToRawLongBits(sum);
+ // Can only set a sticky bit if the bit is not set.
+ if ((hi & 0x1) == 0) {
+ // In a standard extended precision result for (a+b) the bits are extra
+ // magnitude lost and the sticky bit is positive.
+ // Here the round-off magnitude (r) can be negative so the sticky
+ // bit should be added (same sign) or subtracted (different sign).
+ if (sum > 0) {
+ hi += (r > 0) ? 1 : -1;
+ } else {
+ hi += (r < 0) ? 1 : -1;
+ }
+ sum = Double.longBitsToDouble(hi);
+ }
+ }
+ return sum;
+ }
+ }
+
+ /**
+ * Shewchuk's COMPRESS(e) algorithm uses fast-two-sum with a double pass
+ * (big->small then small->big) and zero elimination. This is a variant
+ * of {@link ExactSum}. Methods differ for the inlined versions which use an array
+ * to store non-zero values for an efficient second pass summation. This avoids
+ * the cascading function calls used in {@link ExactSum}. This variant exists for
+ * benchmarking.
+ */
+ private static final class ExactSum2 extends ExactSum {
+ /** An instance. */
+ static final ExactSum2 INSTANCE = new ExactSum2();
+
+ /** No construction. */
+ private ExactSum2() {}
+
+ @Override
+ public double sum(double[] e, int size) {
+ if (size == 1) {
+ return e[0];
+ }
+
+ // First traversal (from big to small)
+ // Shewchuk uses (Q,q) for (big,small) from fast-two-sum and carries Q.
+ // Here use (q,qq); p is used for intermediates.
+ final int m = size - 1;
+ double q = e[m];
+ int bottom = m;
+ for (int i = m - 1; i >= 0; i--) {
+ final double p = q + e[i];
+ final double qq = DoublePrecision.fastTwoSumLow(q, e[i], p);
+ if (qq != 0) {
+ // Store larger component in e and carry the smaller
+ e[bottom] = p;
+ bottom--;
+ q = qq;
+ } else {
+ // Compression.
+ q = p;
+ }
+ }
+ // Redundant store
+ // e[bottom] = q
+
+ return stickySum(q, e, bottom, m);
+ }
+
+ /**
+ * Compute the summation of the parts using a fast-two-sum. The remainder is
+ * carried into the next part using a sticky bit for the final correctly
+ * rounded result.
+ *
+ * @param q the lowest term
+ * @param e the remaining terms from in the range (bottom, m]
+ * @param bottom the bottom of the range (exclusive)
+ * @param m the upper of the range (inclusive)
+ * @return the sum
+ */
+ private static double stickySum(double q, double[] e, int bottom, int m) {
+ if (bottom == m) {
+ // Complete compression to size 1
+ return q;
+ }
+
+ // Second traversal (from small to big)
+ // Here we do not store the round-off parts back into the expansion
+ // but summarise them into the carry using a sticky bit.
+ for (int i = bottom + 1; i < m; i++) {
+ q = fastSumWithStickyBit(e[i], q);
+ }
+ // Final sum. No requirement to compute round-off.
+ return e[m] + q;
+ }
+
+ @Override
+ public double sum(double e0, double e1, double e2, double e3) {
+ int bottom = 2;
+ final double[] e = new double[3];
+
+ // First traversal (from big to small).
+ double q = e3 + e2;
+ double qq = DoublePrecision.fastTwoSumLow(e3, e2, q);
+ if (qq != 0) {
+ e[bottom--] = q;
+ q = qq;
+ }
+ double p = q + e1;
+ qq = DoublePrecision.fastTwoSumLow(q, e1, p);
+ if (qq != 0) {
+ e[bottom--] = p;
+ p = qq;
+ }
+ q = p + e0;
+ qq = DoublePrecision.fastTwoSumLow(p, e0, q);
+ if (qq != 0) {
+ e[bottom--] = q;
+ q = qq;
+ }
+
+ return stickySum(q, e, bottom, 2);
+ }
+
+ @Override
+ public double sum(double e0, double e1, double e2, double e3, double e4, double e5) {
+ int bottom = 4;
+ final double[] e = new double[5];
+
+ // First traversal (from big to small).
+ double q = e5 + e4;
+ double qq = DoublePrecision.fastTwoSumLow(e5, e4, q);
+ if (qq != 0) {
+ e[bottom--] = q;
+ q = qq;
+ }
+ double p = q + e3;
+ qq = DoublePrecision.fastTwoSumLow(q, e3, p);
+ if (qq != 0) {
+ e[bottom--] = p;
+ p = qq;
+ }
+ q = p + e2;
+ qq = DoublePrecision.fastTwoSumLow(p, e2, q);
+ if (qq != 0) {
+ e[bottom--] = q;
+ q = qq;
+ }
+ p = q + e1;
+ qq = DoublePrecision.fastTwoSumLow(q, e1, p);
+ if (qq != 0) {
+ e[bottom--] = p;
+ p = qq;
+ }
+ q = p + e0;
+ qq = DoublePrecision.fastTwoSumLow(p, e0, q);
+ if (qq != 0) {
+ e[bottom--] = q;
+ q = qq;
+ }
+
+ return stickySum(q, e, bottom, 4);
+ }
+
+ @Override
+ public double sum(double e0, double e1, double e2, double e3, double e4, double e5, double e6, double e7) {
+ int bottom = 6;
+ final double[] e = new double[7];
+
+ // First traversal (from big to small).
+ // If round-off is non-zero it is carried down and the high part deposited in g.
+ // If round-off is zero compression has occurred. The high part is carried and
+ // gi remains zero. This does not eliminate zeros from g but spurious zeros are
+ // left to avoid use of array indexing. These are ignored from the sum operations.
+ double q = e7 + e6;
+ double qq = DoublePrecision.fastTwoSumLow(e7, e6, q);
+ if (qq != 0) {
+ e[bottom--] = q;
+ q = qq;
+ }
+ double p = q + e5;
+ qq = DoublePrecision.fastTwoSumLow(q, e5, p);
+ if (qq != 0) {
+ e[bottom--] = p;
+ p = qq;
+ }
+ q = p + e4;
+ qq = DoublePrecision.fastTwoSumLow(p, e4, q);
+ if (qq != 0) {
+ e[bottom--] = q;
+ q = qq;
+ }
+ p = q + e3;
+ qq = DoublePrecision.fastTwoSumLow(q, e3, p);
+ if (qq != 0) {
+ e[bottom--] = p;
+ p = qq;
+ }
+ q = p + e2;
+ qq = DoublePrecision.fastTwoSumLow(p, e2, q);
+ if (qq != 0) {
+ e[bottom--] = q;
+ q = qq;
+ }
+ p = q + e1;
+ qq = DoublePrecision.fastTwoSumLow(q, e1, p);
+ if (qq != 0) {
+ e[bottom--] = p;
+ p = qq;
+ }
+ q = p + e0;
+ qq = DoublePrecision.fastTwoSumLow(p, e0, q);
+ if (qq != 0) {
+ e[bottom--] = q;
+ q = qq;
+ }
+
+ return stickySum(q, e, bottom, 6);
+ }
+ }
+
+ /**
+ * Define the summation algorithm for the expansion.
+ */
+ enum Summation {
+ /** Standard precision sum of the parts. */
+ STANDARD,
+ /** Double-length precision sum of the parts. */
+ DOUBLE,
+ /** Exact sum of the parts. */
+ EXACT
+ }
+
+ /** The summation algorithm. */
+ private final ExpansionSum summation;
+ /** The factory to create a working array. */
+ private final IntFunction<double[]> arrayFactory;
+
+ /**
+ * Private constructor.
+ *
+ * @param summation The summation algorithm.
+ */
+ private ExtendedPrecision(ExpansionSum summation) {
+ this(summation, double[]::new);
+ }
+
+ /**
+ * Private constructor.
+ *
+ * @param summation The summation algorithm.
+ * @param arrayFactory The factory to create a working array.
+ */
+ private ExtendedPrecision(ExpansionSum summation, IntFunction<double[]> arrayFactory) {
+ this.summation = summation;
+ this.arrayFactory = arrayFactory;
+ }
+
+ /**
+ * Create an instance with the specified summation algorithm and array factory.
+ * Thread-safety is dependent on the array factory.
+ *
+ * @param type Summation algorithm
+ * @param arrayFactory The factory to create a working array.
+ * @return the instance
+ */
+ static ExtendedPrecision of(Summation type, IntFunction<double[]> arrayFactory) {
+ switch (type) {
+ case DOUBLE:
+ return new ExtendedPrecision(TwoSum.INSTANCE);
+ case EXACT:
+ return new ExtendedPrecision(ExactSum.INSTANCE);
+ case STANDARD:
+ default:
+ return new ExtendedPrecision(ApproximationSum.INSTANCE);
+ }
+ }
+
+ @Override
+ protected double computeValue(double[] a, double[] b) {
+ // This method uses an optimised two-product to create an expansion for a1*b1 + a2*b2.
+ // This is added to the current sum using an expansion sum.
+ // The initial two-product creates the expansion e.
+ // This may grow to contain 'length' split numbers.
+ // The size is kept smaller using zero elimination.
+ final int len = a.length;
+ final double[] e = arrayFactory.apply(len * 2);
+ int size = sumProduct(a[0], b[0], a[1], b[1], e);
+
+ // Remaining split products added to the current sum.
+ // Index i is the second of the pair
+ for (int i = 3; i < len; i += 2) {
+ // Create the expansion to add inline.
+ // Do not re-use sumProduct to limit array read/writes.
+
+ final double a1 = a[i - 1];
+ final double b1 = b[i - 1];
+ final double a2 = a[i];
+ final double b2 = b[i];
+
+ // Expansion e
+ double e1 = a1 * b1;
+ double e0 = DoublePrecision.productLow(a1, b1, e1);
+
+ // Expansion f
+ final double f1 = a2 * b2;
+ final double f0 = DoublePrecision.productLow(a2, b2, f1);
+
+ // Inline an expansion sum to avoid sorting e and f into a sequence g.
+ // f0 into e
+ double q = e0 + f0;
+ e0 = DoublePrecision.twoSumLow(e0, f0, q);
+ double e2 = e1 + q;
+ e1 = DoublePrecision.twoSumLow(e1, q, e2);
+ // f1 into e
+ q = e1 + f1;
+ e1 = DoublePrecision.twoSumLow(e1, f1, q);
+ final double e3 = e2 + q;
+ e2 = DoublePrecision.twoSumLow(e2, q, e3);
+
+ // Add the round-off parts if non-zero.
+ int n = 0;
+ if (e0 != 0) {
+ growExpansion(e, size++, n++, e0);
+ }
+ if (e1 != 0) {
+ growExpansion(e, size++, n++, e1);
+ }
+ if (e2 != 0) {
+ growExpansion(e, size++, n++, e2);
+ }
+ // Unlikely that the overall representation of the two-product is zero
+ // so no check for non-zero here.
+ growExpansion(e, size++, n, e3);
+
+ size = zeroElimination(e, size);
+ }
+ // Add a trailing final product
+ if ((len & 0x1) == 0x1) {
+ // Create the expansion f.
+ final int i = len - 1;
+ final double a1 = a[i];
+ final double b1 = b[i];
+ final double f1 = a1 * b1;
+ final double f0 = DoublePrecision.productLow(a1, b1, f1);
+ if (f0 != 0) {
+ growExpansion(e, size++, 0, f0);
+ growExpansion(e, size++, 1, f1);
+ } else {
+ growExpansion(e, size++, 0, f1);
+ }
+ // Ignore zero elimination as the result is now summed.
+ }
+
+ // Final summation.
+ if (size == 0) {
+ return 0.0;
+ }
+
+ final double result = summation.sum(e, size);
+ if (!Double.isFinite(result)) {
+ // Either we have split infinite numbers or some coefficients were NaNs,
+ // just rely on the naive implementation and let IEEE754 handle this
+ return standardDotProduct(a, b);
+ }
+ return result;
+ }
+
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2) {
+ // s is the running scalar product used for a final edge-case check
+ double s;
+
+ // Initial product creates the expansion e[0-1]
+ double e1 = a1 * b1;
+ double e0 = DoublePrecision.productLow(a1, b1, e1);
+ s = e1;
+
+ // Second product creates expansion f[0-1]
+ final double f1 = a2 * b2;
+ final double f0 = DoublePrecision.productLow(a2, b2, f1);
+ s += f1;
+ // Expansion sum f into e to create e[0-3]
+ // f0 into e
+ double q = e0 + f0;
+ e0 = DoublePrecision.twoSumLow(e0, f0, q);
+ double e2 = e1 + q;
+ e1 = DoublePrecision.twoSumLow(e1, q, e2);
+ // f1 into e
+ q = e1 + f1;
+ e1 = DoublePrecision.twoSumLow(e1, f1, q);
+ final double e3 = e2 + q;
+ e2 = DoublePrecision.twoSumLow(e2, q, e3);
+
+ // Final summation
+ return getSum(s, summation.sum(e0, e1, e2, e3));
+ }
+
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2,
+ double a3, double b3) {
+ // s is the running scalar product used for a final edge-case check
+ double s;
+
+ // Initial product creates the expansion e[0-1]
+ double e1 = a1 * b1;
+ double e0 = DoublePrecision.productLow(a1, b1, e1);
+ s = e1;
+
+ // Second product creates expansion f[0-1]
+ double f1 = a2 * b2;
+ double f0 = DoublePrecision.productLow(a2, b2, f1);
+ s += f1;
+ // Expansion sum f into e to create e[0-3]
+ // f0 into e
+ double q = e0 + f0;
+ e0 = DoublePrecision.twoSumLow(e0, f0, q);
+ double e2 = e1 + q;
+ e1 = DoublePrecision.twoSumLow(e1, q, e2);
+ // f1 into e
+ q = e1 + f1;
+ e1 = DoublePrecision.twoSumLow(e1, f1, q);
+ double e3 = e2 + q;
+ e2 = DoublePrecision.twoSumLow(e2, q, e3);
+
+ // Third product creates the expansion f[0-1]
+ f1 = a3 * b3;
+ f0 = DoublePrecision.productLow(a3, b3, f1);
+ s += f1;
+
+ // Expansion sum f into e.
+ // f0 into e
+ q = e0 + f0;
+ e0 = DoublePrecision.twoSumLow(e0, f0, q);
+ double p = e1 + q;
+ e1 = DoublePrecision.twoSumLow(e1, q, p);
+ q = e2 + p;
+ e2 = DoublePrecision.twoSumLow(e2, p, q);
+ double e4 = e3 + q;
+ e3 = DoublePrecision.twoSumLow(e3, q, e4);
+
+ // f1 into e
+ q = e1 + f1;
+ e1 = DoublePrecision.twoSumLow(e1, f1, q);
+ p = e2 + q;
+ e2 = DoublePrecision.twoSumLow(e2, q, p);
+ q = e3 + p;
+ e3 = DoublePrecision.twoSumLow(e3, p, q);
+ final double e5 = e4 + q;
+ e4 = DoublePrecision.twoSumLow(e4, q, e5);
+
+ // Final summation
+ return getSum(s, summation.sum(e0, e1, e2, e3, e4, e5));
+ }
+
+ @Override
+ public double value(double a1, double b1,
+ double a2, double b2,
+ double a3, double b3,
+ double a4, double b4) {
+ // s is the running scalar product used for a final edge-case check
+ double s;
+
+ // Initial product creates the expansion e[0-1]
+ double e1 = a1 * b1;
+ double e0 = DoublePrecision.productLow(a1, b1, e1);
+ s = e1;
+
+ // Second product creates expansion f[0-1]
+ double f1 = a2 * b2;
+ double f0 = DoublePrecision.productLow(a2, b2, f1);
+ s += f1;
+ // Expansion sum f into e to create e[0-3]
+ // f0 into e
+ double q = e0 + f0;
+ e0 = DoublePrecision.twoSumLow(e0, f0, q);
+ double e2 = e1 + q;
+ e1 = DoublePrecision.twoSumLow(e1, q, e2);
+ // f1 into e
+ q = e1 + f1;
+ e1 = DoublePrecision.twoSumLow(e1, f1, q);
+ double e3 = e2 + q;
+ e2 = DoublePrecision.twoSumLow(e2, q, e3);
+
+ // Third product creates the expansion f[0-1]
+ f1 = a3 * b3;
+ f0 = DoublePrecision.productLow(a3, b3, f1);
+ s += f1;
+
+ // Second product creates expansion g[0-1]
+ final double g1 = a4 * b4;
+ final double g0 = DoublePrecision.productLow(a4, b4, g1);
+ s += g1;
+ // Expansion sum g into f to create f[0-3]
+ // g0 into f
+ q = f0 + g0;
+ f0 = DoublePrecision.twoSumLow(f0, g0, q);
+ double f2 = f1 + q;
+ f1 = DoublePrecision.twoSumLow(f1, q, f2);
+ // g1 into f
+ q = f1 + g1;
+ f1 = DoublePrecision.twoSumLow(f1, g1, q);
+ final double f3 = f2 + q;
+ f2 = DoublePrecision.twoSumLow(f2, q, f3);
+
+ // Expansion sum f into e.
+ // f0 into e
+ q = e0 + f0;
+ e0 = DoublePrecision.twoSumLow(e0, f0, q);
+ double p = e1 + q;
+ e1 = DoublePrecision.twoSumLow(e1, q, p);
+ q = e2 + p;
+ e2 = DoublePrecision.twoSumLow(e2, p, q);
+ double e4 = e3 + q;
+ e3 = DoublePrecision.twoSumLow(e3, q, e4);
+
+ // f1 into e
+ q = e1 + f1;
+ e1 = DoublePrecision.twoSumLow(e1, f1, q);
+ p = e2 + q;
+ e2 = DoublePrecision.twoSumLow(e2, q, p);
+ q = e3 + p;
+ e3 = DoublePrecision.twoSumLow(e3, p, q);
+ double e5 = e4 + q;
+ e4 = DoublePrecision.twoSumLow(e4, q, e5);
+
+ // f2 into e
+ q = e2 + f2;
+ e2 = DoublePrecision.twoSumLow(e2, f2, q);
+ p = e3 + q;
+ e3 = DoublePrecision.twoSumLow(e3, q, p);
+ q = e4 + p;
+ e4 = DoublePrecision.twoSumLow(e4, p, q);
+ double e6 = e5 + q;
+ e5 = DoublePrecision.twoSumLow(e5, q, e6);
+
+ // f3 into e
+ q = e3 + f3;
+ e3 = DoublePrecision.twoSumLow(e3, f3, q);
+ p = e4 + q;
+ e4 = DoublePrecision.twoSumLow(e4, q, p);
+ q = e5 + p;
+ e5 = DoublePrecision.twoSumLow(e5, p, q);
+ final double e7 = e6 + q;
+ e6 = DoublePrecision.twoSumLow(e6, q, e7);
+
+ // Final summation
+ return getSum(s, summation.sum(e0, e1, e2, e3, e4, e5, e6, e7));
+ }
+
+ /**
+ * Compute the sum of the two products and store the result in the expansion.
+ * Interspersed zeros are removed and the length returned.
+ *
+ * @param a1 First factor of the first term.
+ * @param b1 Second factor of the first term.
+ * @param a2 First factor of the second term.
+ * @param b2 Second factor of the second term.
+ * @param e Expansion
+ * @return the length of the new expansion (can be zero)
+ */
+ private static int sumProduct(double a1, double b1, double a2, double b2, double[] e) {
+ // Expansion e
+ double e1 = a1 * b1;
+ double e0 = DoublePrecision.productLow(a1, b1, e1);
+
+ // Expansion f
+ final double f1 = a2 * b2;
+ final double f0 = DoublePrecision.productLow(a2, b2, f1);
+
+ // Inline an expansion sum to avoid sorting e and f into a sequence g.
+ // f0 into e
+ double q = e0 + f0;
+ e0 = DoublePrecision.twoSumLow(e0, f0, q);
+ double e2 = e1 + q;
+ e1 = DoublePrecision.twoSumLow(e1, q, e2);
+ // f1 into e
+ q = e1 + f1;
+ e1 = DoublePrecision.twoSumLow(e1, f1, q);
+ final double e3 = e2 + q;
+ e2 = DoublePrecision.twoSumLow(e2, q, e3);
+
+ // Store but remove interspersed zeros
+ int ei = 0;
+ if (e0 != 0) {
+ e[ei++] = e0;
+ }
+ if (e1 != 0) {
+ e[ei++] = e1;
+ }
+ if (e2 != 0) {
+ e[ei++] = e2;
+ }
+ // Unlikely that the overall representation of the two-product is zero
+ // so no check for non-zero here.
+ e[ei++] = e3;
+ return ei;
+ }
+
+ /**
+ * Grow the expansion. This maintains the increasing non-overlapping expansion
+ * by two-summing the new value through the entire expansion from the given
+ * start.
+ *
+ * @param expansion Expansion.
+ * @param length Expansion size.
+ * @param start Start point to begin the merge. To be used for optimised
+ * expansion sum.
+ * @param value Value to add.
+ */
+ private static void growExpansion(double[] expansion, int length, int start, double value) {
+ double p = value;
+ for (int i = start; i < length; i++) {
+ final double ei = expansion[i];
+ final double q = ei + p;
+ expansion[i] = DoublePrecision.twoSumLow(ei, p, q);
+ // Carry the larger magnitude up to the next iteration.
+ p = q;
+ }
+ expansion[length] = p;
+ }
+
+ /**
+ * Perform zero elimination on the expansion.
+ * The new size can be zero.
+ *
+ * @param e Expansion.
+ * @param size Expansion size.
+ * @return the new size
+ */
+ private static int zeroElimination(double[] e, int size) {
+ int newSize = 0;
+ // Skip to the first zero
+ while (newSize < size && e[newSize] != 0) {
+ newSize++;
+ }
+ if (newSize != size) {
+ // Skip the zero and copy remaining non-zeros.
+ // This avoids building blocks of non-zeros to bulk-copy using
+ // System.arraycopy. This extra complexity may be useful
+ // if the number of zeros is small compared to the
+ // length of the expansion.
+ for (int i = newSize + 1; i < size; i++) {
+ if (e[i] != 0) {
+ e[newSize++] = e[i];
+ }
+ }
+ }
+ return newSize;
+ }
+ }
+
+ /**
+ * Gets the final sum. This checks the high precision sum is finite, otherwise
+ * returns the standard precision sum for the IEEE754 result.
+ *
+ * <p>The high precision sum may be non-finite due to input infinite
+ * or NaN numbers or overflow in the summation. However the high precision sum
+ * can also be non-finite when the standard sum is finite. This occurs when
+ * the split product had a component high-part that overflowed during
+ * computation of the hx * hy partial result. In all cases returning the
+ * standard sum ensures the IEEE754 result.
+ *
+ * @param sum Standard sum.
+ * @param hpSum High precision sum.
+ * @return the sum
+ */
+ static double getSum(double sum, double hpSum) {
+ if (!Double.isFinite(hpSum)) {
+ // Either we have split infinite numbers, some coefficients were NaNs,
+ // or the split product overflowed.
+ // Return the naive implementation for the IEEE754 result.
+ return sum;
+ }
+ return hpSum;
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/StickySumPerformance.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/StickySumPerformance.java
new file mode 100644
index 0000000..1821b75
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/StickySumPerformance.java
@@ -0,0 +1,609 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.numbers.examples.jmh.arrays;
+
+import org.apache.commons.rng.UniformRandomProvider;
+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;
+
+import java.util.concurrent.TimeUnit;
+import java.util.function.DoubleBinaryOperator;
+
+/**
+ * Executes a benchmark to measure the speed of operations in the {@link LinearCombination} class.
+ * Benchmarks focus on the sticky summation of two double values.
+ *
+ * <p>Details of the sticky bit can be found in:
+ * <blockquote>
+ * Coonen, J.T., "An Implementation Guide to a Proposed Standard for Floating Point
+ * Arithmetic", Computer, Vol. 13, No. 1, Jan. 1980, pp 68-79.
+ * </blockquote>
+ */
+@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 StickySumPerformance {
+ /** The mask for the sign bit and the mantissa. */
+ private static final long SIGN_MATISSA_MASK = 0x800f_ffff_ffff_ffffL;
+
+ /** Constant for no method. */
+ private static final String NONE = "none";
+ /** Constant for branched method. */
+ private static final String BRANCHED = "branched";
+ /** Constant for branchless method. */
+ private static final String BRANCHLESS = "branchless";
+ /** Constant for single branch method based on the low bit of the high part. */
+ private static final String BRANCH_ON_HI = "branch_on_hi";
+ /** Constant for single branch method based on the low part. */
+ private static final String BRANCH_ON_LO = "branch_on_lo";
+
+ /**
+ * The factors to sum.
+ */
+ @State(Scope.Benchmark)
+ public static class BiFactors {
+ /** The exponent for small numbers. */
+ private static final long EXP = Double.doubleToRawLongBits(1.0);
+
+ /**
+ * The count of sums.
+ */
+ @Param({"10000"})
+ private int size;
+
+ /**
+ * The fraction of numbers that have a zero round-off.
+ */
+ @Param({"0", "0.1"})
+ private double zeroRoundoff;
+
+ /** Factors. */
+ private double[] a;
+
+ /**
+ * Gets the factors to be summed as pairs. The array length will be even.
+ *
+ * @return Factors.
+ */
+ public double[] getFactors() {
+ return a;
+ }
+
+ /**
+ * Create the factors.
+ */
+ @Setup
+ public void setup() {
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP);
+ a = new double[size * 2];
+ // Report on the dataset
+ int nonZero = 0;
+ int unsetSticky = 0;
+ for (int i = 0; i < a.length; i += 2) {
+ // Create data similar to the summation of an expansion.
+ // E.g. Generate 2 numbers with large differences in exponents.
+ // Their sum should have a round-off term containing many bits.
+ final double x = nextDouble(rng);
+ // The round-off can conditionally be forced to zero.
+ final double y = (rng.nextDouble() < zeroRoundoff) ?
+ 0.0 :
+ nextDouble(rng) * 0x1.0p52;
+
+ // Initial expansion
+ double e1 = x + y;
+ double e0 = DoublePrecision.twoSumLow(x, y, e1);
+
+ // Validate methods
+ final double expected = fastSumWithStickyBitBranched(e0, e1);
+ assertEqual(expected, fastSumWithStickyBitBranchless(e0, e1), BRANCHLESS);
+ assertEqual(expected, fastSumWithStickyBitBranchedOnHigh(e0, e1), BRANCH_ON_HI);
+ assertEqual(expected, fastSumWithStickyBitBranchedOnLow(e0, e1), BRANCH_ON_LO);
+
+ // Lower parts of expansion for use in sticky sum
+ a[i] = e0;
+ a[i + 1] = e1;
+
+ // Check the sum and round-off
+ final double sum = e1 + e0;
+ final double r = e0 - (sum - e1);
+ if (r != 0) {
+ nonZero++;
+ }
+ if ((Double.doubleToRawLongBits(sum) & 0x1) == 0) {
+ unsetSticky++;
+ }
+ }
+ // CHECKSTYLE: stop Regexp
+ System.out.printf("%n%nNon-zero %d/%d (%.3f) : Unset sticky %d/%d (%.3f)%n%n",
+ nonZero, size,
+ (double) nonZero / size, unsetSticky, size, (double) unsetSticky / size);
+ // CHECKSTYLE: resume Regexp
+ }
+
+ /**
+ * Create the next double in the range [1, 2). All mantissa bits have an equal
+ * probability of being set.
+ *
+ * @param rng Generator of random numbers.
+ * @return the double
+ */
+ private static double nextDouble(UniformRandomProvider rng) {
+ final long bits = rng.nextLong() & SIGN_MATISSA_MASK;
+ return Double.longBitsToDouble(bits | EXP);
+ }
+
+ /**
+ * Assert the two values are equal.
+ *
+ * @param expected Expected value
+ * @param actual Actual value
+ * @param msg Error message
+ */
+ private static void assertEqual(double expected, double actual, String msg) {
+ if (Double.compare(expected, actual) != 0) {
+ throw new IllegalStateException("Methods do not match: " + msg);
+ }
+ }
+ }
+
+ /**
+ * The summation method.
+ */
+ @State(Scope.Benchmark)
+ public static class SumMethod {
+ /**
+ * The name of the method.
+ */
+ @Param({NONE, BRANCHED, BRANCHLESS, BRANCH_ON_HI, BRANCH_ON_LO})
+ private String name;
+
+ /** The function. */
+ private DoubleBinaryOperator fun;
+
+ /**
+ * Gets the function.
+ *
+ * @return the function
+ */
+ public DoubleBinaryOperator getFunction() {
+ return fun;
+ }
+
+ /**
+ * Create the function.
+ */
+ @Setup
+ public void setup() {
+ if (NONE.equals(name)) {
+ fun = (a, b) -> a + b;
+ } else if (BRANCHED.equals(name)) {
+ fun = StickySumPerformance::fastSumWithStickyBitBranched;
+ } else if (BRANCHLESS.equals(name)) {
+ fun = StickySumPerformance::fastSumWithStickyBitBranchless;
+ } else if (BRANCH_ON_HI.equals(name)) {
+ fun = StickySumPerformance::fastSumWithStickyBitBranchedOnHigh;
+ } else if (BRANCH_ON_LO.equals(name)) {
+ fun = StickySumPerformance::fastSumWithStickyBitBranchedOnLow;
+ } else {
+ throw new IllegalStateException("Unknown sum method: " + name);
+ }
+ }
+ }
+
+ /**
+ * Compute the sum of two numbers {@code a} and {@code b} using
+ * Dekker's two-sum algorithm. The values are required to be ordered by magnitude
+ * {@code |a| >= |b|}. The result is adjusted to set the lowest bit as a sticky
+ * bit that summarises the magnitude of the round-off that were lost. The
+ * result is not the correctly rounded result; it is intended the result is to
+ * be used in an addition with a value with a greater magnitude exponent. This
+ * addition will have exact round-to-nearest, ties-to-even rounding taking account
+ * of bits lots in the previous sum.
+ *
+ * @param a First part of sum.
+ * @param b Second part of sum.
+ * @return <code>b - (sum - a)</code>
+ * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+ * Shewchuk (1997) Theorum 6</a>
+ */
+ private static double fastSumWithStickyBitBranched(double a, double b) {
+ double sum = a + b;
+ // bVitual = sum - a
+ // b - bVirtual == b round-off
+ final double r = b - (sum - a);
+
+ if (r != 0) {
+ // Bits will be lost.
+ // In floating-point arithmetic the sticky bit is the bit-wise OR
+ // of the rest of the binary bits that cannot be stored in the
+ // preliminary representation of the result:
+ //
+ // sgn | exp | V | N | .............. | L | G | R | S
+ //
+ // sgn : sign bit
+ // exp : exponent
+ // V : overflow for significand field
+ // N and L : most and least significant bits of the result
+ // G and R : the two bits beyond
+ // S : sticky bit, bitwise OR of all bits thereafter
+ //
+ // The sticky bit is a flag indicating if there is more magnitude beyond
+ // the last bits. Here the round-off is signed so we have to consider the
+ // sign of the sum and round-off together and either add the sticky or
+ // remove it. The final bit is thus used to push up the next addition using
+ // the sum to a higher value, or down to a lower value, when tie breaking for
+ // the correct round-to-nearest, ties-to-even result.
+ long hi = Double.doubleToRawLongBits(sum);
+ // Can only set a sticky bit if the bit is not set.
+ if ((hi & 0x1) == 0) {
+ // In a standard extended precision result for (a+b) the bits are extra
+ // magnitude lost and the sticky bit is positive.
+ // Here the round-off magnitude (r) can be negative so the sticky
+ // bit should be added (same sign) or subtracted (different sign).
+ if (sum > 0) {
+ hi += (r > 0) ? 1 : -1;
+ } else {
+ hi += (r < 0) ? 1 : -1;
+ }
+ sum = Double.longBitsToDouble(hi);
+ }
+ }
+ return sum;
+ }
+
+ /**
+ * Compute the sum of two numbers {@code a} and {@code b} using
+ * Dekker's two-sum algorithm. The values are required to be ordered by magnitude
+ * {@code |a| >= |b|}. The result is adjusted to set the lowest bit as a sticky
+ * bit that summarises the magnitude of the round-off that were lost. The
+ * result is not the correctly rounded result; it is intended the result is to
+ * be used in an addition with a value with a greater magnitude exponent. This
+ * addition will have exact round-to-nearest, ties-to-even rounding taking account
+ * of bits lots in the previous sum.
+ *
+ * @param a First part of sum.
+ * @param b Second part of sum.
+ * @return <code>b - (sum - a)</code>
+ * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+ * Shewchuk (1997) Theorum 6</a>
+ */
+ private static double fastSumWithStickyBitBranchless(double a, double b) {
+ final double sum = a + b;
+ // bVitual = sum - a
+ // b - bVirtual == b round-off
+ final double r = b - (sum - a);
+
+ // In floating-point arithmetic the sticky bit is the bit-wise OR
+ // of the rest of the binary bits that cannot be stored in the
+ // preliminary representation of the result:
+ //
+ // sgn | exp | V | N | .............. | L | G | R | S
+ //
+ // sgn : sign bit
+ // exp : exponent
+ // V : overflow for significand field
+ // N and L : most and least significant bits of the result
+ // G and R : the two bits beyond
+ // S : sticky bit, bitwise OR of all bits thereafter
+ //
+ // The sticky bit is a flag indicating if there is more magnitude beyond
+ // the last bits.
+ //
+ // Here the round-off is signed so we have to consider the
+ // sign of the sum and round-off together and either add the sticky or
+ // remove it. The final bit is thus used to push up the next addition using
+ // the sum to a higher value, or down to a lower value, when tie breaking for
+ // the correct round-to-nearest, ties-to-even result.
+ //
+ // One extra consideration: sum is already rounded. Since we are using the
+ // last bit to store a sticky bit then if the final bit is 1 then this was
+ // not created by a ties-to-even rounding and is already a sticky bit.
+
+ // Compute the sticky bit addition:
+ // sign sum last bit sum sign r magnitude r sticky
+ // x 1 x x +0
+ //
+ // 1 0 1 1 +1
+ // 1 0 1 0 +0
+ // 1 0 0 1 -1
+ // 1 0 0 0 +0
+ // 0 0 1 1 -1
+ // 0 0 1 0 +0
+ // 0 0 0 1 +1
+ // 0 0 0 0 +0
+ //
+ // Magnitude of r is computed by bitwise OR of the 63-bits from exponent+mantissa
+ // Sign of sum and r is the sign-bit of sum or r
+
+ final long hi = Double.doubleToRawLongBits(sum);
+
+ // Note: >50% of the time all code below here is redundant
+ //if ((hi & 0x1) == 0x1) {
+ // // Already sticky
+ // return sum;
+ //}
+
+ final long lo = Double.doubleToRawLongBits(r);
+
+ // OR compress least significant 63-bits into lowest bit
+ long sticky = lo;
+ sticky |= sticky >>> 31; // Discard sign bit
+ sticky |= sticky >>> 16;
+ sticky |= sticky >>> 8;
+ sticky |= sticky >>> 4;
+ sticky |= sticky >>> 2;
+ sticky |= sticky >>> 1; // final sticky bit is in position 0
+
+ // AND with the inverse of the trailing bit from hi to set it to zero
+ // if the last bit in hi is 1 (already sticky).
+ sticky = sticky & ~hi;
+
+ // Clear the rest. Sticky is now 0 if r was 0.0; or 1 if r was non-zero.
+ sticky = sticky & 0x1;
+
+ // The sign bit is created as + or - using the XOR of hi and lo.
+ // Signed shift will create a flag: -1 to negate, else 0.
+ final long fNegate = (hi ^ lo) >> 63;
+
+ // Conditionally negate a value without branching:
+ // http://graphics.stanford.edu/~seander/bithacks.html#ConditionalNegate
+ // (Logic updated since fNegate is already negative.)
+ // (1 ^ 0) - 0 = 1 - 0 = 1
+ // (0 ^ 0) - 0 = 0 - 0 = 0
+ // (1 ^ -1) - -1 = -2 - -1 = -1
+ // (0 ^ -1) - -1 = -1 - -1 = 0
+ sticky = (sticky ^ fNegate) - fNegate;
+
+ return Double.longBitsToDouble(hi + sticky);
+ }
+
+ /**
+ * Compute the sum of two numbers {@code a} and {@code b} using
+ * Dekker's two-sum algorithm. The values are required to be ordered by magnitude
+ * {@code |a| >= |b|}. The result is adjusted to set the lowest bit as a sticky
+ * bit that summarises the magnitude of the round-off that were lost. The
+ * result is not the correctly rounded result; it is intended the result is to
+ * be used in an addition with a value with a greater magnitude exponent. This
+ * addition will have exact round-to-nearest, ties-to-even rounding taking account
+ * of bits lots in the previous sum.
+ *
+ * @param a First part of sum.
+ * @param b Second part of sum.
+ * @return <code>b - (sum - a)</code>
+ * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+ * Shewchuk (1997) Theorum 6</a>
+ */
+ private static double fastSumWithStickyBitBranchedOnHigh(double a, double b) {
+ final double sum = a + b;
+ // bVitual = sum - a
+ // b - bVirtual == b round-off
+ final double r = b - (sum - a);
+
+ // In floating-point arithmetic the sticky bit is the bit-wise OR
+ // of the rest of the binary bits that cannot be stored in the
+ // preliminary representation of the result:
+ //
+ // sgn | exp | V | N | .............. | L | G | R | S
+ //
+ // sgn : sign bit
+ // exp : exponent
+ // V : overflow for significand field
+ // N and L : most and least significant bits of the result
+ // G and R : the two bits beyond
+ // S : sticky bit, bitwise OR of all bits thereafter
+ //
+ // The sticky bit is a flag indicating if there is more magnitude beyond
+ // the last bits.
+ //
+ // Here the round-off is signed so we have to consider the
+ // sign of the sum and round-off together and either add the sticky or
+ // remove it. The final bit is thus used to push up the next addition using
+ // the sum to a higher value, or down to a lower value, when tie breaking for
+ // the correct round-to-nearest, ties-to-even result.
+ //
+ // One extra consideration: sum is already rounded. Since we are using the
+ // last bit to store a sticky bit then if the final bit is 1 then this was
+ // not created by a ties-to-even rounding and is already a sticky bit.
+
+ // Compute the sticky bit addition:
+ // sign sum last bit sum sign r magnitude r sticky
+ // x 1 x x +0
+ //
+ // 1 0 1 1 +1
+ // 1 0 1 0 +0
+ // 1 0 0 1 -1
+ // 1 0 0 0 +0
+ // 0 0 1 1 -1
+ // 0 0 1 0 +0
+ // 0 0 0 1 +1
+ // 0 0 0 0 +0
+ //
+ // Magnitude of r is computed by bitwise OR of the 63-bits from exponent+mantissa
+ // Sign of sum and r is the sign-bit of sum or r
+
+ final long hi = Double.doubleToRawLongBits(sum);
+
+ if ((hi & 0x1) == 0x1) {
+ // Already sticky
+ return sum;
+ }
+
+ final long lo = Double.doubleToRawLongBits(r);
+
+ // OR compress least significant 63-bits into lowest bit
+ long sticky = lo;
+ sticky |= sticky >>> 31; // Discard sign bit
+ sticky |= sticky >>> 16;
+ sticky |= sticky >>> 8;
+ sticky |= sticky >>> 4;
+ sticky |= sticky >>> 2;
+ sticky |= sticky >>> 1; // final sticky bit is in position 0
+
+ // No requirement for AND with the inverse of the trailing bit from hi as we
+ // have eliminated that condition.
+
+ // Clear the rest. Sticky is now 0 if r was 0.0; or 1 if r was non-zero.
+ sticky = sticky & 0x1;
+
+ // The sign bit is created as + or - using the XOR of hi and lo.
+ // Signed shift will create a flag: -1 to negate, else 0.
+ final long fNegate = (hi ^ lo) >> 63;
+
+ // Conditionally negate a value without branching:
+ // http://graphics.stanford.edu/~seander/bithacks.html#ConditionalNegate
+ // (Logic updated since fNegate is already negative.)
+ // (1 ^ 0) - 0 = 1 - 0 = 1
+ // (0 ^ 0) - 0 = 0 - 0 = 0
+ // (1 ^ -1) - -1 = -2 - -1 = -1
+ // (0 ^ -1) - -1 = -1 - -1 = 0
+ sticky = (sticky ^ fNegate) - fNegate;
+
+ return Double.longBitsToDouble(hi + sticky);
+
+ }
+
+ /**
+ * Compute the sum of two numbers {@code a} and {@code b} using
+ * Dekker's two-sum algorithm. The values are required to be ordered by magnitude
+ * {@code |a| >= |b|}. The result is adjusted to set the lowest bit as a sticky
+ * bit that summarises the magnitude of the round-off that were lost. The
+ * result is not the correctly rounded result; it is intended the result is to
+ * be used in an addition with a value with a greater magnitude exponent. This
+ * addition will have exact round-to-nearest, ties-to-even rounding taking account
+ * of bits lots in the previous sum.
+ *
+ * @param a First part of sum.
+ * @param b Second part of sum.
+ * @return <code>b - (sum - a)</code>
+ * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+ * Shewchuk (1997) Theorum 6</a>
+ */
+ private static double fastSumWithStickyBitBranchedOnLow(double a, double b) {
+ final double sum = a + b;
+ // bVitual = sum - a
+ // b - bVirtual == b round-off
+ final double r = b - (sum - a);
+
+ // In floating-point arithmetic the sticky bit is the bit-wise OR
+ // of the rest of the binary bits that cannot be stored in the
+ // preliminary representation of the result:
+ //
+ // sgn | exp | V | N | .............. | L | G | R | S
+ //
+ // sgn : sign bit
+ // exp : exponent
+ // V : overflow for significand field
+ // N and L : most and least significant bits of the result
+ // G and R : the two bits beyond
+ // S : sticky bit, bitwise OR of all bits thereafter
+ //
+ // The sticky bit is a flag indicating if there is more magnitude beyond
+ // the last bits.
+ //
+ // Here the round-off is signed so we have to consider the
+ // sign of the sum and round-off together and either add the sticky or
+ // remove it. The final bit is thus used to push up the next addition using
+ // the sum to a higher value, or down to a lower value, when tie breaking for
+ // the correct round-to-nearest, ties-to-even result.
+ //
+ // One extra consideration: sum is already rounded. Since we are using the
+ // last bit to store a sticky bit then if the final bit is 1 then this was
+ // not created by a ties-to-even rounding and is already a sticky bit.
+
+ // Compute the sticky bit addition.
+ // This is only done when r is non-zero:
+ // sign sum last bit sum sign r sticky
+ // x 1 x +0
+ //
+ // 1 0 1 +1
+ // 1 0 0 -1
+ // 0 0 1 -1
+ // 0 0 0 +1
+ //
+ // Sign of sum and r is the sign-bit of sum or r
+
+ // In the majority of cases there is some round-off.
+ // Testing for non-zero allows the branch to assume the sticky bit magnitude
+ // is 1 (unless the final bit of hi is already set).
+ if (r != 0) {
+ // Bits will be lost.
+ final long hi = Double.doubleToRawLongBits(sum);
+ final long lo = Double.doubleToRawLongBits(r);
+
+ // Can only set a sticky bit if the bit is not already set.
+ // Flip the bits and extract the lowest bit. This is 1 if
+ // the sticky bit is not currently set. If already set then
+ // 'sticky' is zero and the rest of this execution path has
+ // no effect but we have eliminated the requirement to check the
+ // 50/50 branch:
+ // if ((hi & 0x1) == 0) {
+ // // set sticky ...
+ // }
+ int sticky = ~((int) hi) & 0x1;
+
+ // The sign bit is created as + or - using the XOR of hi and lo.
+ // Signed shift will create a flag: -1 to negate, else 0.
+ final int fNegate = (int) ((hi ^ lo) >> 63);
+
+ // Conditionally negate a value without branching:
+ // http://graphics.stanford.edu/~seander/bithacks.html#ConditionalNegate
+ // (Logic updated since fNegate is already negative.)
+ // (1 ^ 0) - 0 = 1 - 0 = 1
+ // (0 ^ 0) - 0 = 0 - 0 = 0
+ // (1 ^ -1) - -1 = -2 - -1 = -1
+ // (0 ^ -1) - -1 = -1 - -1 = 0
+ sticky = (sticky ^ fNegate) - fNegate;
+
+ return Double.longBitsToDouble(hi + sticky);
+ }
+
+ return sum;
+ }
+
+ // Benchmark methods.
+
+ /**
+ * Benchmark the sticky summation of two numbers.
+ *
+ * @param factors Factors.
+ * @param bh Data sink.
+ * @param method Summation method.
+ */
+ @Benchmark
+ public void stickySum(BiFactors factors, Blackhole bh, SumMethod method) {
+ final DoubleBinaryOperator fun = method.getFunction();
+ final double[] a = factors.getFactors();
+ for (int i = 0; i < a.length; i += 2) {
+ bh.consume(fun.applyAsDouble(a[i], a[i + 1]));
+ }
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/package-info.java b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/package-info.java
new file mode 100644
index 0000000..63e94de
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/arrays/package-info.java
@@ -0,0 +1,21 @@
+/*
+ * 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.
+ */
+
+/**
+ * Benchmarks for the {@code org.apache.commons.numbers.arrays} components.
+ */
+package org.apache.commons.numbers.examples.jmh.arrays;
diff --git a/commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/DoublePrecisionTest.java b/commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/DoublePrecisionTest.java
new file mode 100644
index 0000000..deda670
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/DoublePrecisionTest.java
@@ -0,0 +1,93 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.numbers.examples.jmh.arrays;
+
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
+
+/**
+ * Tests for {@link DoublePrecision}.
+ */
+public class DoublePrecisionTest {
+ @Test
+ public void testSplitAssumptions() {
+ // The multiplier used to split the double value into high and low parts.
+ final double scale = (1 << 27) + 1;
+ // The upper limit above which a number may overflow during the split into a high part.
+ final double limit = 0x1.0p996;
+ Assertions.assertTrue(Double.isFinite(limit * scale));
+ Assertions.assertTrue(Double.isFinite(-limit * scale));
+ // Cannot make the limit the next power up
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, limit * 2 * scale);
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, -limit * 2 * scale);
+ // Check the level for the safe upper limit of the exponent of the sum of the absolute
+ // components of the product
+ Assertions.assertTrue(Math.getExponent(2 * Math.sqrt(Double.MAX_VALUE)) - 2 > 508);
+ }
+
+ @Test
+ public void testHighPart() {
+ Assertions.assertEquals(Double.NaN, DoublePrecision.highPart(Double.POSITIVE_INFINITY));
+ Assertions.assertEquals(Double.NaN, DoublePrecision.highPart(Double.NEGATIVE_INFINITY));
+ Assertions.assertEquals(Double.NaN, DoublePrecision.highPart(Double.NaN));
+ // Any finite number should be split to a finite number
+ Assertions.assertTrue(Double.isFinite(DoublePrecision.highPart(Double.MAX_VALUE)));
+ Assertions.assertTrue(Double.isFinite(DoublePrecision.highPart(-Double.MAX_VALUE)));
+ }
+
+ @Test
+ public void testHighPartUnscaled() {
+ Assertions.assertEquals(Double.NaN, DoublePrecision.highPartUnscaled(Double.POSITIVE_INFINITY));
+ Assertions.assertEquals(Double.NaN, DoublePrecision.highPartUnscaled(Double.NEGATIVE_INFINITY));
+ Assertions.assertEquals(Double.NaN, DoublePrecision.highPartUnscaled(Double.NaN));
+ // Large finite numbers will overflow during the split
+ Assertions.assertEquals(Double.NaN, DoublePrecision.highPartUnscaled(Double.MAX_VALUE));
+ Assertions.assertEquals(Double.NaN, DoublePrecision.highPartUnscaled(-Double.MAX_VALUE));
+ }
+
+ /**
+ * Test {@link DoublePrecision#productLow(double, double, double)} computes the same
+ * result as JDK 9 Math.fma(x, y, -x * y) for edge cases.
+ */
+ @Test
+ public void testProductLow() {
+ assertProductLow(0.0, 1.0, Math.nextDown(Double.MIN_NORMAL));
+ assertProductLow(0.0, -1.0, Math.nextDown(Double.MIN_NORMAL));
+ assertProductLow(Double.NaN, 1.0, Double.POSITIVE_INFINITY);
+ assertProductLow(Double.NaN, 1.0, Double.NEGATIVE_INFINITY);
+ assertProductLow(Double.NaN, 1.0, Double.NaN);
+ assertProductLow(0.0, 1.0, Double.MAX_VALUE);
+ assertProductLow(Double.NaN, 2.0, Double.MAX_VALUE);
+ }
+
+ private static void assertProductLow(double expected, double x, double y) {
+ Assertions.assertEquals(expected, DoublePrecision.productLow(x, y, x * y), 0.0);
+ }
+
+ @Test
+ public void testIsNotNormal() {
+ for (double a : new double[] {Double.MAX_VALUE, 1.0, Double.MIN_NORMAL}) {
+ Assertions.assertFalse(DoublePrecision.isNotNormal(a));
+ Assertions.assertFalse(DoublePrecision.isNotNormal(-a));
+ }
+ for (double a : new double[] {Double.POSITIVE_INFINITY, 0.0,
+ Math.nextDown(Double.MIN_NORMAL), Double.NaN}) {
+ Assertions.assertTrue(DoublePrecision.isNotNormal(a));
+ Assertions.assertTrue(DoublePrecision.isNotNormal(-a));
+ }
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinationAccuracyTest.java b/commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinationAccuracyTest.java
new file mode 100644
index 0000000..befef80
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinationAccuracyTest.java
@@ -0,0 +1,241 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.numbers.examples.jmh.arrays;
+
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.ND;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.simple.RandomSource;
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Disabled;
+import org.junit.jupiter.api.Test;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.Arguments;
+import org.junit.jupiter.params.provider.MethodSource;
+
+import java.io.BufferedWriter;
+import java.io.IOException;
+import java.nio.file.Files;
+import java.nio.file.Paths;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.Formatter;
+import java.util.stream.Collectors;
+import java.util.stream.Stream;
+
+/**
+ * Test the accuracy of each implementation of {@link LinearCombination}.
+ */
+public class LinearCombinationAccuracyTest {
+ /**
+ * The dot product length.
+ * This must be the same for the accuracy report as the assertions.
+ * The accuracy report is used to set the approximate bounds for the pass/fail
+ * condition numbers.
+ */
+ private static final int LENGTH = 100;
+
+ /** The number of samples to test the dot product. */
+ private static final int SAMPLES = 10;
+
+ /**
+ * Provide instances of the LinearCombination interface as arguments
+ * along with the condition numbers where instance will pass and fail.
+ * A pass is a mean relative error to the true dot product of {@code < 1e-3}.
+ *
+ * @return the stream
+ */
+ static Stream<Arguments> provideLinearCombination() {
+ return Stream.of(
+ Arguments.of(LinearCombinations.Dekker.INSTANCE, 1e20, 1e30),
+ Arguments.of(LinearCombinations.Dot2s.INSTANCE, 1e20, 1e30),
+ Arguments.of(LinearCombinations.DotK.DOT_3, 1e35, 1e45),
+ Arguments.of(LinearCombinations.DotK.DOT_4, 1e50, 1e65),
+ Arguments.of(LinearCombinations.DotK.DOT_5, 1e65, 1e85),
+ Arguments.of(LinearCombinations.DotK.DOT_6, 1e80, 1e100),
+ Arguments.of(LinearCombinations.DotK.DOT_7, 1e95, 1e115),
+ Arguments.of(LinearCombinations.ExtendedPrecision.INSTANCE, 1e300, -1),
+ Arguments.of(LinearCombinations.ExtendedPrecision.DOUBLE, 1e300, -1),
+ Arguments.of(LinearCombinations.ExtendedPrecision.EXACT, 1e300, -1),
+ Arguments.of(LinearCombinations.Exact.INSTANCE, 1e300, -1)
+ );
+ }
+
+ /**
+ * Assert the dot product function passes and fails at the specified condition numbers.
+ * The average relative error must be below 1e-3 to pass.
+ *
+ * @param fun the function
+ * @param passC the pass condition number
+ * @param failC the fail condition number (set to negative to ignore)
+ */
+ @ParameterizedTest
+ @MethodSource("provideLinearCombination")
+ public void testDotProduct(ND fun, double passC, double failC) {
+ final double[] x = new double[LENGTH];
+ final double[] y = new double[LENGTH];
+ // Fixed seed to consistency
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP, 9283746);
+
+ // Use an average as the actual condition number of the generated dot product
+ // may not be the requested condition number. It will average out at the desired
+ // level and the pass/fail condition bounds should be suitably broad.
+ double sum = 0;
+ for (int i = 0; i < SAMPLES; i++) {
+ final double expected = LinearCombinationUtils.genDot(passC, rng, x, y, null);
+ final double observed = fun.value(x, y);
+ sum += relativeError(expected, observed);
+ }
+ final double error = sum / SAMPLES;
+ Assertions.assertTrue(error < 1e-3, () -> "Expected to pass at C=" + passC + ". Error = " + error);
+ if (failC < 0) {
+ return;
+ }
+ sum = 0;
+ for (int i = 0; i < SAMPLES; i++) {
+ final double expected = LinearCombinationUtils.genDot(failC, rng, x, y, null);
+ final double observed = fun.value(x, y);
+ sum += relativeError(expected, observed);
+ }
+ final double error2 = sum / SAMPLES;
+ Assertions.assertFalse(error2 < 1e-3, () -> "Expected to fail at C=" + failC + ". Error = " + error2);
+ }
+
+ /**
+ * Report the relative error of various implementations. This is not a test.
+ * The method generates dot products with a random condition number over a suitable range.
+ * The relative error of various dot product implementations are saved to a result file.
+ * The file can be used to set the assertion pass/fail levels for
+ * {@link #testDotProduct(ND, double, double)}.
+ *
+ * @throws IOException Signals that an I/O exception has occurred.
+ */
+ @Test
+ @Disabled("This method is used to output a report of the accuracy of implementations.")
+ public void reportRelativeError() throws IOException {
+ // Ogita et al (2005) Figure 6.2 used length=100, samples=1000 with c in 1 to 1e120.
+ // Ogita et al (2005) Figure 6.4 used length=1000, samples=2000 with c in 1 to 1e120.
+ // Here the condition number is in 1e10 to 1e120 as low condition numbers are
+ // computed by all methods with relative error below 1e-16. This uses the shorter
+ // length for speed.
+ final int samples = 2000;
+ // Sample from the log-uniform distribution:
+ final double logA = Math.log(1e10);
+ final double logB = Math.log(1e120);
+
+ final ArrayList<double[]> data = new ArrayList<>(samples);
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_1024_PP);
+ final double[] x = new double[LENGTH];
+ final double[] y = new double[LENGTH];
+ final double[] cc = new double[1];
+
+ // Instances to test and their names (for the report file)
+ final ArrayList<ND> methods = new ArrayList<>();
+ final ArrayList<String> names = new ArrayList<>();
+ addMethod(methods, names, LinearCombinations.Dekker.INSTANCE, "dekker");
+ addMethod(methods, names, LinearCombinations.Dot2s.INSTANCE, "dot2s");
+ addMethod(methods, names, LinearCombinations.DotK.DOT_3, "dot3");
+ addMethod(methods, names, LinearCombinations.DotK.DOT_4, "dot4");
+ addMethod(methods, names, LinearCombinations.DotK.DOT_5, "dot5");
+ addMethod(methods, names, LinearCombinations.DotK.DOT_6, "dot6");
+ addMethod(methods, names, LinearCombinations.DotK.DOT_7, "dot7");
+ addMethod(methods, names, LinearCombinations.ExtendedPrecision.INSTANCE, "extended");
+ addMethod(methods, names, LinearCombinations.ExtendedPrecision.DOUBLE, "extended2");
+ addMethod(methods, names, LinearCombinations.ExtendedPrecision.EXACT, "extended_exact");
+ addMethod(methods, names, LinearCombinations.Exact.INSTANCE, "exact");
+
+ for (int i = 0; i < samples; i++) {
+ // Random condition number.
+ // This should be approximately uniform over the logarithm of the range.
+ final double u = rng.nextDouble();
+ final double c = Math.exp(u * logB + (1 - u) * logA);
+ final double dot = LinearCombinationUtils.genDot(c, rng, x, y, cc);
+ // Compute with different methods.
+ // Store the condition number of the data first.
+ final double[] e = new double[methods.size() + 1];
+ int j = 0;
+ e[j++] = cc[0];
+ for (final ND method : methods) {
+ e[j++] = compute(x, y, dot, method);
+ }
+ data.add(e);
+ }
+
+ // Sort by condition number
+ Collections.sort(data, (a, b) -> Double.compare(a[0], b[0]));
+
+ // Write to file in the Maven build directory
+ try (BufferedWriter out = Files.newBufferedWriter(Paths.get("target/dot.csv"))) {
+ out.append("Condition no," + names.stream().collect(Collectors.joining(",")));
+ out.newLine();
+ final StringBuilder sb = new StringBuilder(200);
+ try (Formatter formatter = new Formatter(sb)) {
+ for (final double[] e : data) {
+ sb.setLength(0);
+ formatter.format("%.3g", e[0]);
+ for (int i = 1; i < e.length; i++) {
+ formatter.format(",%.3g", e[i]);
+ }
+ out.append(sb);
+ out.newLine();
+ }
+ }
+ }
+ }
+
+ /**
+ * Adds the method to the lists
+ *
+ * @param methods List of methods.
+ * @param names List of method names.
+ * @param method Method implementation.
+ * @param name Method name.
+ */
+ private static void addMethod(ArrayList<ND> methods, ArrayList<String> names,
+ ND method, String name) {
+ methods.add(method);
+ names.add(name);
+ }
+
+ /**
+ * Compute the dot product using the given function and return the relative error.
+ *
+ * @param x the x
+ * @param y the y
+ * @param expected the expected value
+ * @param fun the function
+ * @return the relative error
+ */
+ private static double compute(double[] x, double[] y, double expected, ND fun) {
+ return relativeError(expected, fun.value(x, y));
+ }
+
+ /**
+ * Compute the relative error:
+ * <pre>
+ * |expected − observed| / |expected|
+ * </pre>
+ *
+ * <p>The value is clipped to a maximum of 2.
+ *
+ * @param observed the observed
+ * @param expected the expected
+ * @return the double
+ */
+ private static double relativeError(double expected, double observed) {
+ return Math.min(2, Math.abs(observed - expected) / Math.abs(expected));
+ }
+}
diff --git a/commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinationsTest.java b/commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinationsTest.java
new file mode 100644
index 0000000..f2c99d1
--- /dev/null
+++ b/commons-numbers-examples/examples-jmh/src/test/java/org/apache/commons/numbers/examples/jmh/arrays/LinearCombinationsTest.java
@@ -0,0 +1,588 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.numbers.examples.jmh.arrays;
+
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.FourD;
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.ND;
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.ThreeD;
+import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.TwoD;
+import org.apache.commons.numbers.fraction.BigFraction;
+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;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.Arguments;
+import org.junit.jupiter.params.provider.MethodSource;
+
+import java.math.BigDecimal;
+import java.math.MathContext;
+import java.math.RoundingMode;
+import java.util.SplittableRandom;
+import java.util.stream.Stream;
+
+/**
+ * Test each implementation of the LinearCombination interface.
+ */
+public class LinearCombinationsTest {
+ /** Double.MIN_VALUE as a BigDecimal. Use string constructor to truncate precision to 4.9e-324. */
+ private static final BigDecimal MIN = BigDecimal.valueOf(Double.MIN_VALUE);
+
+ /**
+ * Provide instances of the LinearCombination interface as arguments.
+ *
+ * @return the stream
+ */
+ static Stream<Arguments> provideLinearCombination() {
+ return Stream.of(
+ Arguments.of(LinearCombinations.Dekker.INSTANCE),
+ Arguments.of(LinearCombinations.Dot2s.INSTANCE),
+ Arguments.of(LinearCombinations.DotK.DOT_3),
+ Arguments.of(LinearCombinations.DotK.DOT_4),
+ Arguments.of(LinearCombinations.DotK.DOT_5),
+ Arguments.of(LinearCombinations.DotK.DOT_6),
+ Arguments.of(LinearCombinations.DotK.DOT_7),
+ Arguments.of(LinearCombinations.ExtendedPrecision.INSTANCE),
+ Arguments.of(LinearCombinations.ExtendedPrecision.DOUBLE),
+ Arguments.of(LinearCombinations.ExtendedPrecision.EXACT),
+ Arguments.of(LinearCombinations.ExtendedPrecision.EXACT2),
+ Arguments.of(LinearCombinations.Exact.INSTANCE)
+ );
+ }
+
+ @ParameterizedTest
+ @MethodSource("provideLinearCombination")
+ public void testSingleElementArray(ND fun) {
+ final double[] a = {1.23456789};
+ final double[] b = {98765432.1};
+
+ Assertions.assertEquals(a[0] * b[0], fun.value(a, b));
+ }
+
+ @ParameterizedTest
+ @MethodSource("provideLinearCombination")
+ public void testTwoSums(ND fun) {
+ final BigFraction[] aF = new BigFraction[] {
+ BigFraction.of(-1321008684645961L, 268435456L),
+ BigFraction.of(-5774608829631843L, 268435456L),
+ BigFraction.of(-7645843051051357L, 8589934592L)
+ };
+ final BigFraction[] bF = new BigFraction[] {
+ BigFraction.of(-5712344449280879L, 2097152L),
+ BigFraction.of(-4550117129121957L, 2097152L),
+ BigFraction.of(8846951984510141L, 131072L)
+ };
+ final ThreeD fun3 = (ThreeD) fun;
+
+ final int len = aF.length;
+ final double[] a = new double[len];
+ final double[] b = new double[len];
+ for (int i = 0; i < len; i++) {
+ a[i] = aF[i].getNumerator().doubleValue() / aF[i].getDenominator().doubleValue();
+ b[i] = bF[i].getNumerator().doubleValue() / bF[i].getDenominator().doubleValue();
+ }
+
+ // Ensure "array" and "inline" implementations give the same result.
+ final double abSumInline = fun3.value(a[0], b[0],
+ a[1], b[1],
+ a[2], b[2]);
+ final double abSumArray = fun.value(a, b);
+ Assertions.assertEquals(abSumInline, abSumArray);
+
+ // Compare with arbitrary precision computation.
+ BigFraction result = BigFraction.ZERO;
+ for (int i = 0; i < a.length; i++) {
+ result = result.add(aF[i].multiply(bF[i]));
+ }
+ final double expected = result.doubleValue();
+ Assertions.assertEquals(expected, abSumInline, "Expecting exact result");
+
+ final double naive = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+ Assertions.assertTrue(Math.abs(naive - abSumInline) > 1.5);
+ }
+
+ @ParameterizedTest
+ @MethodSource("provideLinearCombination")
+ public void testHuge(ND fun) {
+ final int scale = 971;
+ final double[] a = new double[] {
+ -1321008684645961.0 / 268435456.0,
+ -5774608829631843.0 / 268435456.0,
+ -7645843051051357.0 / 8589934592.0
+ };
+ final double[] b = new double[] {
+ -5712344449280879.0 / 2097152.0,
+ -4550117129121957.0 / 2097152.0,
+ 8846951984510141.0 / 131072.0
+ };
+ final ThreeD fun3 = (ThreeD) fun;
+
+ final int len = a.length;
+ final double[] scaledA = new double[len];
+ final double[] scaledB = new double[len];
+ for (int i = 0; i < len; ++i) {
+ scaledA[i] = Math.scalb(a[i], -scale);
+ scaledB[i] = Math.scalb(b[i], scale);
+ }
+ final double abSumInline = fun3.value(scaledA[0], scaledB[0],
+ scaledA[1], scaledB[1],
+ scaledA[2], scaledB[2]);
+ final double abSumArray = fun.value(scaledA, scaledB);
+
+ Assertions.assertEquals(abSumInline, abSumArray);
+ Assertions.assertEquals(-1.8551294182586248737720779899, abSumInline, "Expecting exact result");
+
+ final double naive = scaledA[0] * scaledB[0] + scaledA[1] * scaledB[1] + scaledA[2] * scaledB[2];
+ Assertions.assertTrue(Math.abs(naive - abSumInline) > 1.5);
+ }
+
+ @ParameterizedTest
+ @MethodSource("provideLinearCombination")
+ public void testArrayVsInline(ND fun) {
+ // Assume the instance implements the inline functions
+ final TwoD fun2 = (TwoD) fun;
+ final ThreeD fun3 = (ThreeD) fun;
+ final FourD fun4 = (FourD) fun;
+
+ final SplittableRandom rng = new SplittableRandom();
+
+ double sInline;
+ double sArray;
+ final double scale = 1e17;
+ for (int i = 0; i < 1000; ++i) {
+ final double u1 = scale * rng.nextDouble();
+ final double u2 = scale * rng.nextDouble();
+ final double u3 = scale * rng.nextDouble();
+ final double u4 = scale * rng.nextDouble();
+ final double v1 = scale * rng.nextDouble();
+ final double v2 = scale * rng.nextDouble();
+ final double v3 = scale * rng.nextDouble();
+ final double v4 = scale * rng.nextDouble();
+
+ // One sum.
+ sInline = fun2.value(u1, v1, u2, v2);
+ sArray = fun.value(new double[] {u1, u2},
+ new double[] {v1, v2});
+ Assertions.assertEquals(sInline, sArray);
+
+ // Two sums.
+ sInline = fun3.value(u1, v1, u2, v2, u3, v3);
+ sArray = fun.value(new double[] {u1, u2, u3},
+ new double[] {v1, v2, v3});
+ Assertions.assertEquals(sInline, sArray);
+
+ // Three sums.
+ sInline = fun4.value(u1, v1, u2, v2, u3, v3, u4, v4);
+ sArray = fun.value(new double[] {u1, u2, u3, u4},
+ new double[] {v1, v2, v3, v4});
+ Assertions.assertEquals(sInline, sArray);
+ }
+ }
+
+ @ParameterizedTest
+ @MethodSource("provideLinearCombination")
+ public void testNonFinite(ND fun) {
+ final double[][] a = new double[][] {
+ {1, 2, 3, 4},
+ {1, Double.POSITIVE_INFINITY, 3, 4},
+ {1, 2, Double.POSITIVE_INFINITY, 4},
+ {1, Double.POSITIVE_INFINITY, 3, Double.NEGATIVE_INFINITY},
+ {1, 2, 3, 4},
+ {1, 2, 3, 4},
+ {1, 2, 3, 4},
+ {1, 2, 3, 4},
+ {1, Double.MAX_VALUE, 3, 4},
+ {1, 2, Double.MAX_VALUE, 4},
+ {1, Double.MAX_VALUE / 2, 3, -Double.MAX_VALUE / 4},
+ };
+ final double[][] b = new double[][] {
+ {1, -2, 3, 4},
+ {1, -2, 3, 4},
+ {1, -2, 3, 4},
+ {1, -2, 3, 4},
+ {1, Double.POSITIVE_INFINITY, 3, 4},
+ {1, -2, Double.POSITIVE_INFINITY, 4},
+ {1, Double.POSITIVE_INFINITY, 3, Double.NEGATIVE_INFINITY},
+ {Double.NaN, -2, 3, 4},
+ {1, -2, 3, 4},
+ {1, -2, 3, 4},
+ {1, -2, 3, 4},
+ };
+
+ final TwoD fun2 = (TwoD) fun;
+ final ThreeD fun3 = (ThreeD) fun;
+ final FourD fun4 = (FourD) fun;
+
+ Assertions.assertEquals(-3, fun2.value(a[0][0], b[0][0],
+ a[0][1], b[0][1]));
+ Assertions.assertEquals(6, fun3.value(a[0][0], b[0][0],
+ a[0][1], b[0][1],
+ a[0][2], b[0][2]));
+ Assertions.assertEquals(22, fun4.value(a[0][0], b[0][0],
+ a[0][1], b[0][1],
+ a[0][2], b[0][2],
+ a[0][3], b[0][3]));
+ Assertions.assertEquals(22, fun.value(a[0], b[0]));
+
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun2.value(a[1][0], b[1][0],
+ a[1][1], b[1][1]));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun3.value(a[1][0], b[1][0],
+ a[1][1], b[1][1],
+ a[1][2], b[1][2]));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun4.value(a[1][0], b[1][0],
+ a[1][1], b[1][1],
+ a[1][2], b[1][2],
+ a[1][3], b[1][3]));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun.value(a[1], b[1]));
+
+ Assertions.assertEquals(-3, fun2.value(a[2][0], b[2][0],
+ a[2][1], b[2][1]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun3.value(a[2][0], b[2][0],
+ a[2][1], b[2][1],
+ a[2][2], b[2][2]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun4.value(a[2][0], b[2][0],
+ a[2][1], b[2][1],
+ a[2][2], b[2][2],
+ a[2][3], b[2][3]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun.value(a[2], b[2]));
+
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun2.value(a[3][0], b[3][0],
+ a[3][1], b[3][1]));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun3.value(a[3][0], b[3][0],
+ a[3][1], b[3][1],
+ a[3][2], b[3][2]));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun4.value(a[3][0], b[3][0],
+ a[3][1], b[3][1],
+ a[3][2], b[3][2],
+ a[3][3], b[3][3]));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun.value(a[3], b[3]));
+
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun2.value(a[4][0], b[4][0],
+ a[4][1], b[4][1]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun3.value(a[4][0], b[4][0],
+ a[4][1], b[4][1],
+ a[4][2], b[4][2]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun4.value(a[4][0], b[4][0],
+ a[4][1], b[4][1],
+ a[4][2], b[4][2],
+ a[4][3], b[4][3]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun.value(a[4], b[4]));
+
+ Assertions.assertEquals(-3, fun2.value(a[5][0], b[5][0],
+ a[5][1], b[5][1]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun3.value(a[5][0], b[5][0],
+ a[5][1], b[5][1],
+ a[5][2], b[5][2]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun4.value(a[5][0], b[5][0],
+ a[5][1], b[5][1],
+ a[5][2], b[5][2],
+ a[5][3], b[5][3]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun.value(a[5], b[5]));
+
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun2.value(a[6][0], b[6][0],
+ a[6][1], b[6][1]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun3.value(a[6][0], b[6][0],
+ a[6][1], b[6][1],
+ a[6][2], b[6][2]));
+ Assertions.assertEquals(Double.NaN, fun4.value(a[6][0], b[6][0],
+ a[6][1], b[6][1],
+ a[6][2], b[6][2],
+ a[6][3], b[6][3]));
+ Assertions.assertEquals(Double.NaN, fun.value(a[6], b[6]));
+
+ Assertions.assertEquals(Double.NaN, fun2.value(a[7][0], b[7][0],
+ a[7][1], b[7][1]));
+ Assertions.assertEquals(Double.NaN, fun3.value(a[7][0], b[7][0],
+ a[7][1], b[7][1],
+ a[7][2], b[7][2]));
+ Assertions.assertEquals(Double.NaN, fun4.value(a[7][0], b[7][0],
+ a[7][1], b[7][1],
+ a[7][2], b[7][2],
+ a[7][3], b[7][3]));
+ Assertions.assertEquals(Double.NaN, fun.value(a[7], b[7]));
+
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun2.value(a[8][0], b[8][0],
+ a[8][1], b[8][1]));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun3.value(a[8][0], b[8][0],
+ a[8][1], b[8][1],
+ a[8][2], b[8][2]));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun4.value(a[8][0], b[8][0],
+ a[8][1], b[8][1],
+ a[8][2], b[8][2],
+ a[8][3], b[8][3]));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun.value(a[8], b[8]));
+
+ Assertions.assertEquals(-3, fun2.value(a[9][0], b[9][0],
+ a[9][1], b[9][1]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun3.value(a[9][0], b[9][0],
+ a[9][1], b[9][1],
+ a[9][2], b[9][2]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun4.value(a[9][0], b[9][0],
+ a[9][1], b[9][1],
+ a[9][2], b[9][2],
+ a[9][3], b[9][3]));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, fun.value(a[9], b[9]));
+
+ Assertions.assertEquals(-Double.MAX_VALUE, fun2.value(a[10][0], b[10][0],
+ a[10][1], b[10][1]));
+ Assertions.assertEquals(-Double.MAX_VALUE, fun3.value(a[10][0], b[10][0],
+ a[10][1], b[10][1],
+ a[10][2], b[10][2]));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun4.value(a[10][0], b[10][0],
+ a[10][1], b[10][1],
+ a[10][2], b[10][2],
+ a[10][3], b[10][3]));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun.value(a[10], b[10]));
+ }
+
+ /**
+ * This creates a scenario where the split product will overflow but the standard
+ * precision computation will not. The result is expected to be in extended precision,
+ * i.e. the method correctly detects and handles intermediate overflow.
+ *
+ * <p>Note: This test assumes that LinearCombination computes a split number
+ * using Dekker's method. This can result in the high part of the number being
+ * greater in magnitude than the the original number due to round-off in the split.
+ */
+ @Test
+ public void testOverflow() {
+ // Create a simple dot product that is different in high precision and has
+ // values that create a high part above the original number. This can be done using
+ // a mantissa with almost all bits set to 1.
+ final double x = Math.nextDown(2.0);
+ final double y = -Math.nextDown(x);
+ final double xxMxy = x * x + x * y;
+ final double xxMxyHighPrecision = LinearCombinations.DotK.DOT_3.value(x, x, x, y);
+ Assertions.assertNotEquals(xxMxy, xxMxyHighPrecision, "High precision result should be different");
+
+ // Scale it close to max value.
+ // The current exponent is 0 so the combined scale must be 1023-1 as the
+ // partial product x*x and x*y have an exponent 1 higher
+ Assertions.assertEquals(0, Math.getExponent(x));
+ Assertions.assertEquals(0, Math.getExponent(y));
+
+ final double a1 = Math.scalb(x, 1022 - 30);
+ final double b1 = Math.scalb(x, 30);
+ final double a2 = a1;
+ final double b2 = Math.scalb(y, 30);
+ // Verify low precision result is scaled and finite
+ final double sxxMxy = Math.scalb(xxMxy, 1022);
+ Assertions.assertEquals(sxxMxy, a1 * b1 + a2 * b2);
+ Assertions.assertTrue(Double.isFinite(sxxMxy));
+
+ // High precision result.
+ // First demonstrate that Dekker's split will create overflow in the high part.
+ final double m = (1 << 27) + 1;
+ double c;
+ c = a1 * m;
+ final double ha1 = c - (c - a1);
+ c = b1 * m;
+ final double hb1 = c - (c - b1);
+ c = a2 * m;
+ final double ha2 = c - (c - a2);
+ c = b2 * m;
+ final double hb2 = c - (c - b2);
+ Assertions.assertTrue(Double.isFinite(ha1));
+ Assertions.assertTrue(Double.isFinite(hb1));
+ Assertions.assertTrue(Double.isFinite(ha2));
+ Assertions.assertTrue(Double.isFinite(hb2));
+ // High part should be bigger in magnitude
+ Assertions.assertTrue(Math.abs(ha1) > Math.abs(a1));
+ Assertions.assertTrue(Math.abs(hb1) > Math.abs(b1));
+ Assertions.assertTrue(Math.abs(ha2) > Math.abs(a2));
+ Assertions.assertTrue(Math.abs(hb2) > Math.abs(b2));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, ha1 * hb1, "Expected split high part to overflow");
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, ha2 * hb2, "Expected split high part to overflow");
+
+ // LinearCombination should detect and handle intermediate overflow and return the
+ // high precision result.
+ final double expected = Math.scalb(xxMxyHighPrecision, 1022);
+ Assertions.assertEquals(expected, LinearCombinations.DotK.DOT_3.value(a1, b1, a2, b2));
+ Assertions.assertEquals(expected, LinearCombinations.DotK.DOT_3.value(a1, b1, a2, b2, 0, 0));
+ Assertions.assertEquals(expected, LinearCombinations.DotK.DOT_3.value(a1, b1, a2, b2, 0, 0, 0, 0));
+ Assertions.assertEquals(expected, LinearCombinations.DotK.DOT_3.value(new double[] {a1, a2}, new double[] {b1, b2}));
+ }
+
+ /**
+ * This is an extreme case of the sum x^2 + y^2 - 1 when x^2 + y^2 are 1.0 within
+ * floating-point error but if performed using high precision subtracting 1.0 is not 0.0.
+ * This case is derived from computations on a complex cis number.
+ */
+ @Test
+ public void testCisNumber() {
+ final double theta = 5.992112452678286E-7;
+ final double x = Math.cos(theta);
+ final double y = Math.sin(theta);
+ assertValue(LinearCombinations.DotK.DOT_3.value(x, x, y, y, 1, -1),
+ new double[] {x, y, 1},
+ new double[] {x, y, -1});
+ }
+
+ /**
+ * Test the sum of vectors composed of sub-vectors of [a1, a2, a3, a4] * [a1, a2, -a3, -a4]
+ * where a1^2 + a2^2 = 1 and a3^2 + a4^2 = 1 such that the sum is approximately 0 every
+ * 4 products. This is a test that is failed by various implementations that accumulate the
+ * round-off sum in single or 2-fold precision.
+ */
+ @Test
+ public void testSumZero() {
+ // Fixed seed for stability
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PP, 876543L);
+ final int size = 10;
+ // Create random doublets of pairs of numbers that sum to 1 or -1.
+ for (int length = 4; length <= 12; length += 4) {
+ final double[] a = new double[length];
+ final double[] b = new double[length];
+ for (int i = 0; i < size; i++) {
+ // Flip-flop above and below zero
+ double sign = 1;
+ for (int k = 0; k < length; k += 4) {
+ // Create 2 complex cis numbers
+ final double theta1 = rng.nextDouble() * Math.PI / 2;
+ final double theta2 = rng.nextDouble() * Math.PI / 2;
+ a[k + 0] = b[k + 0] = Math.cos(theta1);
+ a[k + 1] = b[k + 1] = Math.sin(theta1);
+ a[k + 2] = b[k + 2] = Math.cos(theta2);
+ a[k + 3] = b[k + 3] = Math.sin(theta2);
+ a[k + 0] *= sign;
+ a[k + 1] *= sign;
+ a[k + 2] *= sign;
+ a[k + 3] *= sign;
+ // Invert second pair.
+ // The sum of the pairs should be zero +/- floating point error.
+ a[k + 2] = -a[k + 2];
+ a[k + 3] = -a[k + 3];
+ sign = -sign;
+ }
+ assertValue(LinearCombinations.DotK.DOT_3.value(a, b), a, b);
+ }
+ }
+ }
+
+ /**
+ * Compute the sum of the product of factors in arbitrary precision and compare it to the
+ * given value.
+ *
+ * @param value the value
+ * @param a factors
+ * @param b factors
+ */
+ private static void assertValue(double value, double[] a, double[] b) {
+ final double expected = computeValue(a, b);
+ Assertions.assertEquals(expected, value, Math.ulp(expected),
+ () -> "Difference in Ulps = " + ulps(expected, value));
+ }
+
+ /**
+ * Compute the sum of the product of pairs of input data using BigDecimal.
+ * The BigDecimal is not allowed to underflow Double.MIN_VALUE.
+ *
+ * @param data the data
+ * @return the sum of products
+ */
+ private static double computeValue(double[] a, double[] b) {
+ BigDecimal sum = new BigDecimal(a[0]).multiply(new BigDecimal(b[0]));
+ for (int i = 1; i < a.length; i++) {
+ sum = clip(sum.add(clip(new BigDecimal(a[i]).multiply(new BigDecimal(b[i])))));
+ }
+ return sum.doubleValue();
+ }
+
+ /**
+ * Compute the units of least precision (ulps) between the two numbers.
+ *
+ * @param a first number
+ * @param b second number
+ * @return the ulps
+ */
+ private static long ulps(double a, double b) {
+ long x = Double.doubleToLongBits(a);
+ long y = Double.doubleToLongBits(b);
+ if (x != y) {
+ if ((x ^ y) < 0L) {
+ // Opposite signs. Measure the combined distance to zero.
+ if (x < 0) {
+ final long tmp = x;
+ x = y;
+ y = tmp;
+ }
+ return (x - Double.doubleToLongBits(0.0)) + (y - Double.doubleToLongBits(-0.0)) + 1;
+ }
+ return Math.abs(x - y);
+ }
+ return 0;
+ }
+
+ /**
+ * Clip the value to the minimum value that can be stored by a double.
+ * Ideally this should round BigDecimal to values occupied by sub-normal numbers.
+ * That is non-trivial so this just removes excess precision in the significand and
+ * clips it to Double.MIN_VALUE or zero if the value is very small. The ultimate use for
+ * the BigDecimal is rounded to the closest double so this method is adequate. It would
+ * take many summations of extended precision sub-normal numbers to create more
+ * than a few ULP difference to the final double value
+ *
+ * <p>In data output by the various tests the values have never been known to require
+ * clipping so this is just a safety threshold.
+ *
+ * @param a the value
+ * @return the clipped value
+ */
+ private static BigDecimal clip(BigDecimal a) {
+ // Min value is approx 4.9e-324. Anything with fewer decimal digits to the right of the
+ // decimal point is OK.
+ if (a.scale() < 324) {
+ return a;
+ }
+ // Reduce the scale
+ final BigDecimal b = a.setScale(MIN.scale(), RoundingMode.HALF_UP);
+ // Clip to min value
+ final BigDecimal bb = b.abs();
+ if (bb.compareTo(MIN) < 0) {
+ // Note the number may be closer to MIN than zero so do rounding
+ if (MIN.subtract(bb).compareTo(bb) < 0) {
+ // Closer to MIN
+ return a.signum() == -1 ? MIN.negate() : MIN;
+ }
+ // Closer to zero
+ return BigDecimal.ZERO;
+ }
+ // Anything above min is allowed.
+ return b;
+ }
+
+ /**
+ * Test the clip method does what it specifies.
+ */
+ @Test
+ public void testClip() {
+ // min value is not affected
+ Assertions.assertEquals(Double.MIN_VALUE, clip(MIN).doubleValue());
+ Assertions.assertEquals(-Double.MIN_VALUE, clip(MIN.negate()).doubleValue());
+ // Round-up to min
+ Assertions.assertEquals(Double.MIN_VALUE, clip(MIN.divide(new BigDecimal(2))).doubleValue());
+ Assertions.assertEquals(-Double.MIN_VALUE, clip(MIN.negate().divide(new BigDecimal(2))).doubleValue());
+ // Round down to zero
+ Assertions.assertEquals(0, clip(MIN.divide(new BigDecimal(2.1), MathContext.DECIMAL64)).doubleValue());
+ Assertions.assertEquals(0, clip(MIN.negate().divide(new BigDecimal(2.1), MathContext.DECIMAL64)).doubleValue());
+ // It does not matter if BigDecimal is more precise than a sub-normal number
+ // when the output is ultimately rounded to a double.
+ Assertions.assertEquals(Double.MIN_VALUE, clip(MIN.multiply(new BigDecimal(1.1))).doubleValue());
+ Assertions.assertEquals(Double.MIN_VALUE, clip(MIN.multiply(new BigDecimal(1.5))).doubleValue());
+ Assertions.assertEquals(Double.MIN_VALUE * 2, clip(MIN.multiply(new BigDecimal(1.6))).doubleValue());
+ }
+}
diff --git a/src/main/resources/checkstyle/checkstyle-suppressions.xml b/src/main/resources/checkstyle/checkstyle-suppressions.xml
index ee631c4..d50f34a 100644
--- a/src/main/resources/checkstyle/checkstyle-suppressions.xml
+++ b/src/main/resources/checkstyle/checkstyle-suppressions.xml
@@ -20,8 +20,9 @@
"https://checkstyle.org/dtds/suppressions_1_2.dtd">
<suppressions>
<suppress checks="Indentation" files=".*[/\\]combinatorics[/\\]Factorial\.java" />
- <suppress checks="ParameterNumber" files=".*[/\\]arrays[/\\]LinearCombination\.java" />
+ <suppress checks="ParameterNumber" files=".*[/\\]arrays[/\\]LinearCombination[s]?\.java" />
<suppress checks="FileLengthCheck" files=".*[/\\]Complex(Test)?\.java" />
+ <suppress checks="FileLengthCheck" files=".*jmh[/\\]arrays[/\\]LinearCombinations.java" />
<suppress checks="MethodLength" files=".*[/\\]Complex\.java" />
<!-- Be more lenient on tests. -->
@@ -32,4 +33,5 @@
<suppress checks="IllegalCatch" files=".*[/\\]test[/\\].*" />
<suppress checks="MethodName" files=".*[/\\]test[/\\].*" />
<suppress checks="ConstantName" files=".*[/\\]test[/\\].*" />
+ <suppress checks="MethodLength" files=".*/LinearCombination.*Test.java" />
</suppressions>