/*
    * 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.core;
  
  /**
   * 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) {
         // bVirtual = 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);
     }
 }