8077587: BigInteger Roots

src/java.base/share/classes/java/math/BigInteger.java

@@ -2743,6 +2743,78 @@ public BigInteger[] sqrtAndRemainder() {
         return new BigInteger[] { sqrtRem[0].toBigInteger(), sqrtRem[1].toBigInteger() };
     }
 
+    /**
+     * Returns the integer {@code n}th root of this BigInteger. The integer
+     * {@code n}th root of the corresponding mathematical integer {@code x} has the
+     * same sign of {@code x}, and its magnitude is the largest integer {@code r}
+     * such that {@code r**n <= abs(x)}. It is equal to the value of
+     * {@code (x.signum() * floor(abs(nthRoot(x, n))))}, where {@code nthRoot(x, n)}
+     * denotes the real {@code n}th root of {@code x} treated as a real. If {@code n}
+     * is even and this BigInteger is negative, an {@code ArithmeticException} will be
+     * thrown.
+     *
+     * <p>Note that the magnitude of the integer {@code n}th root will be less than
+     * the magnitude of the real {@code n}th root if the latter is not representable
+     * as an integral value.
+     *
+     * @param n the root degree
+     * @return the integer {@code n}th root of {@code this}
+     * @throws ArithmeticException if {@code n == 0} (Zeroth roots are not
+     *                             defined.)
+     * @throws ArithmeticException if {@code n} is negative. (This would cause the
+     *                             operation to yield a non-integer value.)
+     * @throws ArithmeticException if {@code n} is even and {@code this} is
+     *                             negative. (This would cause the operation to
+     *                             yield non-real roots.)
+     * @see #sqrt()
+     * @since 26
+     */
+    public BigInteger nthRoot(int n) {
+        return n == 1 ? this : (n == 2 ? sqrt() : nthRootAndRemainder(n, false)[0]);
+    }
+
+    /**
+     * Returns an array of two BigIntegers containing the integer {@code n}th root
+     * {@code r} of {@code this} and its remainder {@code this - r^n},
+     * respectively.
+     *
+     * @param n the root degree
+     * @return an array of two BigIntegers with the integer {@code n}th root at
+     *         offset 0 and the remainder at offset 1
+     * @throws ArithmeticException if {@code n == 0} (Zeroth roots are not
+     *                             defined.)
+     * @throws ArithmeticException if {@code n} is negative. (This would cause the
+     *                             operation to yield a non-integer value.)
+     * @throws ArithmeticException if {@code n} is even and {@code this} is
+     *                             negative. (This would cause the operation to
+     *                             yield non-real roots.)
+     * @see #sqrt()
+     * @see #sqrtAndRemainder()
+     * @see #nthRoot(int)
+     * @since 26
+     */
+    public BigInteger[] nthRootAndRemainder(int n) {
+        return n == 1 ? new BigInteger[] { this, ZERO }
+                      : (n == 2 ? sqrtAndRemainder() : nthRootAndRemainder(n, true));
+    }
+
+    /**
+     * Assume {@code n != 1 && n != 2}
+     */
+    private BigInteger[] nthRootAndRemainder(int n, boolean needRemainder) {
+        if (n <= 0)
+            throw new ArithmeticException("Non-positive root degree");
+
+        if ((n & 1) == 0 && this.signum < 0)
+            throw new ArithmeticException("Negative radicand with even root degree");
+
+        MutableBigInteger[] rootRem = new MutableBigInteger(this.mag).nthRootRem(n);
+        return new BigInteger[] {
+                rootRem[0].toBigInteger(signum),
+                needRemainder ? rootRem[1].toBigInteger(signum) : null
+        };
+    }
+
     /**
      * Returns a BigInteger whose value is the greatest common divisor of
      * {@code abs(this)} and {@code abs(val)}.  Returns 0 if

src/java.base/share/classes/java/math/MutableBigInteger.java

@@ -158,6 +158,27 @@ private void init(int val) {
         value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen);
     }
 
+    /**
+     * Returns a MutableBigInteger with a magnitude specified by
+     * the absolute value of the double val. Any fractional part is discarded.
+     *
+     * Assume val is in the finite double range.
+     */
+    static MutableBigInteger valueOf(double val) {
+        val = Math.abs(val);
+        if (val < 0x1p63)
+            return new MutableBigInteger((long) val);
+        // Translate the double into exponent and significand, according
+        // to the formulae in JLS, Section 20.10.22.
+        long valBits = Double.doubleToRawLongBits(val);
+        int exponent = (int) ((valBits >> 52) & 0x7ffL) - 1075;
+        long significand = (valBits & ((1L << 52) - 1)) | (1L << 52);
+        // At this point, val == significand * 2^exponent, with exponent > 0
+        MutableBigInteger result = new MutableBigInteger(significand);
+        result.leftShift(exponent);
+        return result;
+    }
+
     /**
      * Makes this number an {@code n}-int number all of whose bits are ones.
      * Used by Burnikel-Ziegler division.
@@ -1892,6 +1913,136 @@ private boolean unsignedLongCompare(long one, long two) {
         return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
     }
 
+    /**
+     * Calculate the integer {@code n}th root {@code floor(nthRoot(this, n))} and the remainder,
+     * where {@code nthRoot(., n)} denotes the mathematical {@code n}th root.
+     * The contents of {@code this} are <em>not</em> changed. The value of {@code this}
+     * is assumed to be non-negative and the root degree {@code n >= 3}.
+     *
+     * @implNote The implementation is based on the material in Richard P. Brent
+     * and Paul Zimmermann, <a href="https://maths-people.anu.edu.au/~brent/pd/mca-cup-0.5.9.pdf">
+     * Modern Computer Arithmetic</a>, 27-28.
+     *
+     * @return the integer {@code n}th root of {@code this} and the remainder
+     */
+    MutableBigInteger[] nthRootRem(int n) {
+        // Special cases.
+        if (this.isZero() || this.isOne())
+            return new MutableBigInteger[] { this, new MutableBigInteger() };
+
+        final int bitLength = (int) this.bitLength();
+        // if this < 2^n, result is unity
+        if (bitLength <= n) {
+            MutableBigInteger rem = new MutableBigInteger(this);
+            rem.subtract(ONE);
+            return new MutableBigInteger[] { new MutableBigInteger(1), rem };
+        }
+
+        MutableBigInteger r;
+        if (bitLength <= Long.SIZE) {
+            // Initial estimate is the root of the unsigned long value.
+            final long x = this.toLong();
+            // Use fp arithmetic to get an upper bound of the root
+            final double rad = Math.nextUp(x >= 0 ? x : x + 0x1p64);
+            final double approx = n == 3 ? Math.cbrt(rad) : Math.pow(rad, Math.nextUp(1.0 / n));
+            long rLong = (long) Math.ceil(Math.nextUp(approx));
+
+            if (BigInteger.bitLengthForLong(rLong) * (n - 1) <= Long.SIZE) {
+                // Refine the estimate.
+                long r1 = rLong, rToN1;
+                do {
+                    rLong = r1;
+                    rToN1 = BigInteger.unsignedLongPow(rLong, n - 1);
+                    r1 = ((n - 1) * rLong + Long.divideUnsigned(x, rToN1)) / n;
+                } while (r1 < rLong); // Terminate when non-decreasing.
+
+                return new MutableBigInteger[] {
+                        new MutableBigInteger(rLong), new MutableBigInteger(x - rToN1 * rLong)
+                };
+            } else { // r^(n - 1) could overflow long range, use MutableBigInteger loop instead
+                r = new MutableBigInteger(rLong);
+            }
+        } else {
+            // Set up the initial estimate of the iteration.
+            // Determine a right shift that is a multiple of n into finite double range.
+            long shift;
+            double rad;
+            if (bitLength > Double.MAX_EXPONENT) {
+                shift = bitLength - Double.MAX_EXPONENT;
+                int shiftExcess = (int) (shift % n);
+
+                // Shift the value into finite double range
+                rad = this.toBigInteger().shiftRight((int) shift).doubleValue();
+                // Complete the shift to a multiple of n,
+                // avoiding to lose more bits than necessary.
+                if (shiftExcess != 0) {
+                    int shiftLack = n - shiftExcess;
+                    shift += shiftLack; // shift is long, no overflow
+                    rad /= Double.parseDouble("0x1p" + shiftLack);
+                }
+            } else {
+                shift = 0L;
+                rad = this.toBigInteger().doubleValue();
+            }
+
+            // Use the root of the shifted value as an estimate.
+            rad = Math.nextUp(rad);
+            double approx = n == 3 ? Math.cbrt(rad) : Math.pow(rad, Math.nextUp(1.0 / n));
+            approx = Math.ceil(Math.nextUp(approx));
+            if (shift == 0L) {
+                r = valueOf(approx);
+            } else {
+                // Allocate sufficient space to store the final root
+                r = new MutableBigInteger(new int[(intLen - 1) / n + 1]);
+                r.copyValue(valueOf(approx));
+
+                // Refine the estimate, avoiding to compute non-significant bits
+                final int trailingZeros = this.getLowestSetBit();
+                int rootShift = (int) (shift / n);
+                for (int rootBits = (int) r.bitLength(); rootShift >= rootBits; rootBits <<= 1) {
+                    r.leftShift(rootBits);
+                    rootShift -= rootBits;
+
+                    // Remove useless bits from the radicand
+                    MutableBigInteger x = new MutableBigInteger(this);
+                    int removedBits = rootShift * n;
+                    x.rightShift(removedBits);
+                    if (removedBits > trailingZeros)
+                        x.add(ONE); // round up to ensure r is an upper bound of the root
+
+                    newtonRecurrenceNthRoot(x, r, n, r.toBigInteger().pow(n - 1));
+                }
+
+                // Shift the approximate root back into the original range.
+                r.safeLeftShift(rootShift);
+            }
+        }
+
+        // Refine the estimate.
+        do {
+            BigInteger rBig = r.toBigInteger();
+            BigInteger rToN1 = rBig.pow(n - 1);
+            MutableBigInteger rem = new MutableBigInteger(rToN1.multiply(rBig).mag);
+            if (rem.subtract(this) <= 0)
+                return new MutableBigInteger[] { r, rem };
+
+            newtonRecurrenceNthRoot(this, r, n, rToN1);
+        } while (true);
+    }
+
+    /**
+     * Computes {@code ((n-1)*r + x/rToN1)/n} and places the result in {@code r}.
+     */
+    private static void newtonRecurrenceNthRoot(
+            MutableBigInteger x, MutableBigInteger r, int n, BigInteger rToN1) {
+        MutableBigInteger dividend = new MutableBigInteger();
+        r.mul(n - 1, dividend);
+        MutableBigInteger xDivRToN1 = new MutableBigInteger();
+        x.divide(new MutableBigInteger(rToN1.mag), xDivRToN1, false);
+        dividend.add(xDivRToN1);
+        dividend.divideOneWord(n, r);
+    }
+
     /**
      * Calculate the integer square root {@code floor(sqrt(this))} and the remainder
      * if needed, where {@code sqrt(.)} denotes the mathematical square root.