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 <das@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>