IntMath.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.statistics.descriptive;

  19. import java.math.BigDecimal;
  20. import java.math.BigInteger;
  21. import org.apache.commons.numbers.core.DD;

  23. /**
  24.  * Support class for integer math.
  25.  *
  26.  * @since 1.1
  27.  */
  28. final class IntMath {
  29.     /** Mask for the lower 32-bits of a long. */
  30.     private static final long MASK32 = 0xffff_ffffL;
  31.     /** Mask for the lower 52-bits of a long. */
  32.     private static final long MASK52 = 0xf_ffff_ffff_ffffL;
  33.     /** Bias offset for the exponent of a double. */
  34.     private static final int EXP_BIAS = 1023;
  35.     /** Shift for the exponent of a double. */
  36.     private static final int EXP_SHIFT = 52;
  37.     /** 0.5. */
  38.     private static final double HALF = 0.5;
  39.     /** 2^53. */
  40.     private static final long TWO_POW_53 = 1L << 53;

  42.     /** No instances. */
  43.     private IntMath() {}

  45.     /**
  46.      * Square the values as if an unsigned 64-bit long to produce the high 64-bits
  47.      * of the 128-bit unsigned result.
  48.      *a
  49.      * <p>This method computes the equivalent of:
  50.      * <pre>{@code
  51.      * Math.multiplyHigh(x, x)
  52.      * Math.unsignedMultiplyHigh(x, x) - (((x >> 63) & x) << 1)
  53.      * }</pre>
  54.      *
  55.      * <p>Note: The method {@code Math.multiplyHigh} was added in JDK 9
  56.      * and should be used as above when the source code targets Java 11
  57.      * to exploit the intrinsic method.
  58.      *
  59.      * <p>Note: The method uses the unsigned multiplication. When the input is negative
  60.      * it can be adjusted to the signed result by subtracting the argument twice from the
  61.      * result.
  62.      *
  63.      * @param x Value
  64.      * @return the high 64-bits of the 128-bit result
  65.      */
  66.     static long squareHigh(long x) {
  67.         // Computation is based on the following observation about the upper (a and x)
  68.         // and lower (b and y) bits of unsigned big-endian integers:
  69.         //   ab * xy
  70.         // =  b *  y
  71.         // +  b * x0
  72.         // + a0 *  y
  73.         // + a0 * x0
  74.         // = b * y
  75.         // + b * x * 2^32
  76.         // + a * y * 2^32
  77.         // + a * x * 2^64
  78.         //
  79.         // Summation using a character for each byte:
  80.         //
  81.         //             byby byby
  82.         // +      bxbx bxbx 0000
  83.         // +      ayay ayay 0000
  84.         // + axax axax 0000 0000
  85.         //
  86.         // The summation can be rearranged to ensure no overflow given
  87.         // that the result of two unsigned 32-bit integers multiplied together
  88.         // plus two full 32-bit integers cannot overflow 64 bits:
  89.         // > long x = (1L << 32) - 1
  90.         // > x * x + x + x == -1 (all bits set, no overflow)
  91.         //
  92.         // The carry is a composed intermediate which will never overflow:
  93.         //
  94.         //             byby byby
  95.         // +           bxbx 0000
  96.         // +      ayay ayay 0000
  97.         //
  98.         // +      bxbx 0000 0000
  99.         // + axax axax 0000 0000

  101.         final long a = x >>> 32;
  102.         final long b = x & MASK32;

  104.         final long aa = a * a;
  105.         final long ab = a * b;
  106.         final long bb = b * b;

  108.         // Cannot overflow
  109.         final long carry = (bb >>> 32) +
  110.                            (ab & MASK32) +
  111.                             ab;
  112.         // Note:
  113.         // low = (carry << 32) | (bb & MASK32)
  114.         // Benchmarking shows outputting low to a long[] output argument
  115.         // has no benefit over computing 'low = value * value' separately.

  117.         final long hi = (ab >>> 32) + (carry >>> 32) + aa;
  118.         // Adjust to the signed result:
  119.         // if x < 0:
  120.         //    hi - 2 * x
  121.         return hi - (((x >> 63) & x) << 1);
  122.     }

  124.     /**
  125.      * Multiply the two values as if unsigned 64-bit longs to produce the high 64-bits
  126.      * of the 128-bit unsigned result.
  127.      *
  128.      * <p>This method computes the equivalent of:
  129.      * <pre>{@code
  130.      * Math.multiplyHigh(a, b) + ((a >> 63) & b) + ((b >> 63) & a)
  131.      * }</pre>
  132.      *
  133.      * <p>Note: The method {@code Math.multiplyHigh} was added in JDK 9
  134.      * and should be used as above when the source code targets Java 11
  135.      * to exploit the intrinsic method.
  136.      *
  137.      * <p>Note: The method {@code Math.unsignedMultiplyHigh} was added in JDK 18
  138.      * and should be used when the source code target allows.
  139.      *
  140.      * <p>Taken from {@code o.a.c.rng.core.source64.LXMSupport}.
  141.      *
  142.      * @param value1 the first value
  143.      * @param value2 the second value
  144.      * @return the high 64-bits of the 128-bit result
  145.      */
  146.     static long unsignedMultiplyHigh(long value1, long value2) {
  147.         // Computation is based on the following observation about the upper (a and x)
  148.         // and lower (b and y) bits of unsigned big-endian integers:
  149.         //   ab * xy
  150.         // =  b *  y
  151.         // +  b * x0
  152.         // + a0 *  y
  153.         // + a0 * x0
  154.         // = b * y
  155.         // + b * x * 2^32
  156.         // + a * y * 2^32
  157.         // + a * x * 2^64
  158.         //
  159.         // Summation using a character for each byte:
  160.         //
  161.         //             byby byby
  162.         // +      bxbx bxbx 0000
  163.         // +      ayay ayay 0000
  164.         // + axax axax 0000 0000
  165.         //
  166.         // The summation can be rearranged to ensure no overflow given
  167.         // that the result of two unsigned 32-bit integers multiplied together
  168.         // plus two full 32-bit integers cannot overflow 64 bits:
  169.         // > long x = (1L << 32) - 1
  170.         // > x * x + x + x == -1 (all bits set, no overflow)
  171.         //
  172.         // The carry is a composed intermediate which will never overflow:
  173.         //
  174.         //             byby byby
  175.         // +           bxbx 0000
  176.         // +      ayay ayay 0000
  177.         //
  178.         // +      bxbx 0000 0000
  179.         // + axax axax 0000 0000

  181.         final long a = value1 >>> 32;
  182.         final long b = value1 & MASK32;
  183.         final long x = value2 >>> 32;
  184.         final long y = value2 & MASK32;

  186.         final long by = b * y;
  187.         final long bx = b * x;
  188.         final long ay = a * y;
  189.         final long ax = a * x;

  191.         // Cannot overflow
  192.         final long carry = (by >>> 32) +
  193.                            (bx & MASK32) +
  194.                             ay;
  195.         // Note:
  196.         // low = (carry << 32) | (by & INT_TO_UNSIGNED_BYTE_MASK)
  197.         // Benchmarking shows outputting low to a long[] output argument
  198.         // has no benefit over computing 'low = value1 * value2' separately.

  200.         return (bx >>> 32) + (carry >>> 32) + ax;
  201.     }

  203.     /**
  204.      * Multiply the arguments as if unsigned integers to a {@code double} result.
  205.      *
  206.      * @param x Value.
  207.      * @param y Value.
  208.      * @return the double
  209.      */
  210.     static double unsignedMultiplyToDouble(long x, long y) {
  211.         final long lo = x * y;
  212.         // Fast case: check the arguments cannot overflow a long.
  213.         // This is true if neither has the upper 33-bits set.
  214.         if (((x | y) >>> 31) == 0) {
  215.             // Implicit conversion to a double.
  216.             return lo;
  217.         }
  218.         return uint128ToDouble(unsignedMultiplyHigh(x, y), lo);
  219.     }

  221.     /**
  222.      * Convert an unsigned 128-bit integer to a {@code double}.
  223.      *
  224.      * @param hi High 64-bits.
  225.      * @param lo Low 64-bits.
  226.      * @return the double
  227.      */
  228.     static double uint128ToDouble(long hi, long lo) {
  229.         // Require the representation:
  230.         // 2^exp * mantissa / 2^53
  231.         // The mantissa has an implied leading 1-bit.

  233.         // We have the mantissa final bit as xxx0 or xxx1.
  234.         // To perform correct rounding we maintain the 54-th bit (a) and
  235.         // a check bit (b) of remaining bits.
  236.         // Cases:
  237.         // xxx0 00 - round-down              [1]
  238.         // xxx0 0b - round-down              [1]
  239.         // xxx0 a0 - half-even, round-down   [4]
  240.         // xxx0 ab - round-up                [2]
  241.         // xxx1 00 - round-down              [1]
  242.         // xxx1 0b - round-down              [1]
  243.         // xxx1 a0 - half-even, round-up     [3]
  244.         // xxx1 ab - round-up                [2]
  245.         // [1] If the 54-th bit is 0 always round-down.
  246.         // [2] Otherwise round-up if the check bit is set or
  247.         // [3] the final bit is odd (half-even rounding up).
  248.         // [4] half-even rounding down.

  250.         if (hi == 0) {
  251.             // If lo is a 63-bit result then we are done
  252.             if (lo >= 0) {
  253.                 return lo;
  254.             }
  255.             // Create a 63-bit number with a sticky bit for rounding, rescale the result
  256.             return 2 * (double) ((lo >>> 1) | (lo & 0x1));
  257.         }

  259.         // Initially we create the most significant 64-bits.
  260.         final int shift = Long.numberOfLeadingZeros(hi);
  261.         // Shift the high bits and add trailing low bits.
  262.         // The mask is for the bits from low that are *not* used.
  263.         // Flipping the mask obtains the bits we concatenate
  264.         // after shifting (64 - shift).
  265.         final long maskLow = -1L >>> shift;
  266.         long bits64 = (hi << shift) | ((lo & ~maskLow) >>> -shift);
  267.         // exponent for 2^exp is the index of the highest bit in the 128 bit integer
  268.         final int exp = 127 - shift;
  269.         // Some of the low bits are lost. If non-zero set
  270.         // a sticky bit for rounding.
  271.         bits64 |= (lo & maskLow) == 0 ? 0 : 1;

  273.         // We have a 64-bit unsigned fraction magnitude and an exponent.
  274.         // This must be converted to a IEEE double by mapping the fraction to a base of 2^53.

  276.         // Create the 53-bit mantissa without the implicit 1-bit
  277.         long bits = (bits64 >>> 11) & MASK52;
  278.         // Extract 54-th bit and a sticky bit
  279.         final long a = (bits64 >>> 10) & 0x1;
  280.         final long b = (bits64 << 54) == 0 ? 0 : 1;
  281.         // Perform half-even rounding.
  282.         bits += a & (b | (bits & 0x1));

  284.         // Add the exponent.
  285.         // No worry about overflow to the sign bit as the max exponent is 127.
  286.         bits += (long) (exp + EXP_BIAS) << EXP_SHIFT;

  288.         return Double.longBitsToDouble(bits);
  289.     }

  291.     /**
  292.      * Return the whole number that is nearest to the {@code double} argument {@code x}
  293.      * as an {@code int}, with ties rounding towards positive infinity.
  294.      *
  295.      * <p>This will raise an {@link ArithmeticException} if the closest
  296.      * integer result is not within the range {@code [-2^31, 2^31)},
  297.      * i.e. it overflows an {@code int}; or the argument {@code x}
  298.      * is not finite.
  299.      *
  300.      * <p>Note: This method is equivalent to:
  301.      * <pre>
  302.      * Math.toIntExact(Math.round(x))
  303.      * </pre>
  304.      *
  305.      * <p>The behaviour has been re-implemented for consistent error handling
  306.      * for {@code int}, {@code long} and {@code BigInteger} types.
  307.      *
  308.      * @param x Value.
  309.      * @return rounded value
  310.      * @throws ArithmeticException if the {@code result} overflows an {@code int},
  311.      * or {@code x} is not finite
  312.      * @see Math#round(double)
  313.      * @see Math#toIntExact(long)
  314.      */
  315.     static int toIntExact(double x) {
  316.         final double r = roundToInteger(x);
  317.         if (r >= -0x1.0p31 && r < 0x1.0p31) {
  318.             return (int) r;
  319.         }
  320.         throw new ArithmeticException("integer overflow: " + x);
  321.     }

  323.     /**
  324.      * Return the whole number that is nearest to the {@code double} argument {@code x}
  325.      * as an {@code long}, with ties rounding towards positive infinity.
  326.      *
  327.      * <p>This will raise an {@link ArithmeticException} if the closest
  328.      * integer result is not within the range {@code [-2^63, 2^63)},
  329.      * i.e. it overflows a {@code long}; or the argument {@code x}
  330.      * is not finite.
  331.      *
  332.      * @param x Value.
  333.      * @return rounded value
  334.      * @throws ArithmeticException if the {@code result} overflows a {@code long},
  335.      * or {@code x} is not finite
  336.      */
  337.     static long toLongExact(double x) {
  338.         final double r = roundToInteger(x);
  339.         if (r >= -0x1.0p63 && r < 0x1.0p63) {
  340.             return (long) r;
  341.         }
  342.         throw new ArithmeticException("long integer overflow: " + x);
  343.     }

  345.     /**
  346.      * Return the whole number that is nearest to the {@code double} argument {@code x}
  347.      * as an {@code int}, with ties rounding towards positive infinity.
  348.      *
  349.      * <p>This will raise an {@link ArithmeticException} if the argument {@code x}
  350.      * is not finite.
  351.      *
  352.      * @param x Value.
  353.      * @return rounded value
  354.      * @throws ArithmeticException if {@code x} is not finite
  355.      */
  356.     static BigInteger toBigIntegerExact(double x) {
  357.         if (!Double.isFinite(x)) {
  358.             throw new ArithmeticException("BigInteger overflow: " + x);
  359.         }
  360.         final double r = roundToInteger(x);
  361.         if (r >= -0x1.0p63 && r < 0x1.0p63) {
  362.             // Representable as a long
  363.             return BigInteger.valueOf((long) r);
  364.         }
  365.         // Large result
  366.         return new BigDecimal(r).toBigInteger();
  367.     }

  369.     /**
  370.      * Get the whole number that is the nearest to x, with ties rounding towards positive infinity.
  371.      *
  372.      * <p>This method is intended to perform the equivalent of
  373.      * {@link Math#round(double)} without converting to a {@code long} primitive type.
  374.      * This allows the domain of the result to be checked against the range {@code [-2^63, 2^63)}.
  375.      *
  376.      * <p>Note: Adapted from {@code o.a.c.math4.AccurateMath.rint} and
  377.      * modified to perform rounding towards positive infinity.
  378.      *
  379.      * @param x Number from which nearest whole number is requested.
  380.      * @return a double number r such that r is an integer {@code r - 0.5 <= x < r + 0.5}
  381.      */
  382.     private static double roundToInteger(double x) {
  383.         final double y = Math.floor(x);
  384.         final double d = x - y;
  385.         if (d >= HALF) {
  386.             // Here we do not preserve the sign of the operand in the case
  387.             // of -0.5 < x <= -0.0 since the rounded result is required as an integer.
  388.             // if y == -1.0:
  389.             //    return -0.0
  390.             return y + 1.0;
  391.         }
  392.         return y;
  393.     }

  395.     /**
  396.      * Divide value {@code x} by the count {@code n}.
  397.      *
  398.      * @param x Value.
  399.      * @param n Count.
  400.      * @return the quotient
  401.      */
  402.     static double divide(Int128 x, long n) {
  403.         final DD a = x.toDD();
  404.         if (n < TWO_POW_53) {
  405.             // n is a representable double
  406.             return a.divide(n).doubleValue();
  407.         }
  408.         // Extended precision divide when n > 2^53
  409.         return a.divide(DD.of(n)).doubleValue();
  410.     }
  411. }
Created with JaCoCo 0.8.12.202403310830