From ab12fbfda2c364bb16ddf03b923989639f437f6a Mon Sep 17 00:00:00 2001 From: Fabio Romano Date: Mon, 8 Sep 2025 16:10:22 +0000 Subject: [PATCH] 8077587: BigInteger Roots Reviewed-by: rgiulietti --- .../share/classes/java/math/BigInteger.java | 79 +++++++ .../classes/java/math/MutableBigInteger.java | 200 +++++++++++++++++- .../java/math/BigInteger/BigIntegerTest.java | 157 +++++++++++++- 3 files changed, 433 insertions(+), 3 deletions(-) diff --git a/src/java.base/share/classes/java/math/BigInteger.java b/src/java.base/share/classes/java/math/BigInteger.java index 51d935f10c1..21f8598266f 100644 --- a/src/java.base/share/classes/java/math/BigInteger.java +++ b/src/java.base/share/classes/java/math/BigInteger.java @@ -2742,6 +2742,85 @@ public class BigInteger extends Number implements Comparable { 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 {@code r} of the corresponding mathematical integer {@code x} + * is defined as follows: + * + * If the root is defined, it is equal to the value of + * {@code x.signum()}⋅ ⌊{@code |nthRoot(x, n)|}⌋, + * where {@code nthRoot(x, n)} denotes the real {@code n}th root of {@code x} + * treated as a real. + * Otherwise, the method throws an {@code ArithmeticException}. + * + *

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}. + * @throws ArithmeticException if {@code n} is even and {@code this} is negative. + * @see #sqrt() + * @since 26 + * @apiNote Note that calling {@code nthRoot(2)} is equivalent to calling {@code sqrt()}. + */ + public BigInteger nthRoot(int n) { + if (n == 1) + return this; + + if (n == 2) + return sqrt(); + + checkRootDegree(n); + return new MutableBigInteger(this.mag).nthRootRem(n)[0].toBigInteger(signum); + } + + /** + * Returns an array of two BigIntegers containing the integer {@code n}th root + * {@code r} of {@code this} and its remainder {@code this - r}{@code 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}. + * @throws ArithmeticException if {@code n} is even and {@code this} is negative. + * @see #sqrt() + * @see #sqrtAndRemainder() + * @see #nthRoot(int) + * @since 26 + * @apiNote Note that calling {@code nthRootAndRemainder(2)} is equivalent to calling + * {@code sqrtAndRemainder()}. + */ + public BigInteger[] nthRootAndRemainder(int n) { + if (n == 1) + return new BigInteger[] { this, ZERO }; + + if (n == 2) + return sqrtAndRemainder(); + + checkRootDegree(n); + MutableBigInteger[] rootRem = new MutableBigInteger(this.mag).nthRootRem(n); + return new BigInteger[] { + rootRem[0].toBigInteger(signum), + rootRem[1].toBigInteger(signum) + }; + } + + private void checkRootDegree(int n) { + 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"); + } + /** * Returns a BigInteger whose value is the greatest common divisor of * {@code abs(this)} and {@code abs(val)}. Returns 0 if diff --git a/src/java.base/share/classes/java/math/MutableBigInteger.java b/src/java.base/share/classes/java/math/MutableBigInteger.java index 6ff435ba1ed..dd1da29ddd2 100644 --- a/src/java.base/share/classes/java/math/MutableBigInteger.java +++ b/src/java.base/share/classes/java/math/MutableBigInteger.java @@ -1,5 +1,5 @@ /* - * Copyright (c) 1999, 2024, Oracle and/or its affiliates. All rights reserved. + * Copyright (c) 1999, 2025, Oracle and/or its affiliates. All rights reserved. * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. * * This code is free software; you can redistribute it and/or modify it @@ -1892,6 +1892,204 @@ class MutableBigInteger { 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 not changed. The value of {@code this} + * is assumed to be non-negative and the root degree {@code n >= 3}. + * Assumes {@code bitLength() <= Integer.MAX_VALUE}. + * + * @implNote The implementation is based on the material in Richard P. Brent + * and Paul Zimmermann, + * Modern Computer Arithmetic, p. 27-28. + * + * @param n the root degree + * @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 s; + if (bitLength <= Long.SIZE) { + // Initial estimate is the root of the unsigned long value. + final long x = this.toLong(); + long sLong = (long) nthRootApprox(Math.nextUp(x >= 0 ? x : x + 0x1p64), n) + 1L; + /* The integer-valued recurrence formula in the algorithm of Brent&Zimmermann + * simply discards the fraction part of the real-valued Newton recurrence + * on the function f discussed in the referenced work. + * Indeed, for real x and integer n > 0, the equality ⌊x/n⌋ == ⌊⌊x⌋/n⌋ holds, + * from which the claim follows. + * As a consequence, an initial underestimate (not discussed in BZ) + * will immediately lead to a (weak) overestimate during the 1st iteration, + * thus meeting BZ requirements for termination and correctness. + */ + if (BigInteger.bitLengthForLong(sLong) * (n - 1) <= Long.SIZE) { + // Do the 1st iteration outside the loop to ensure an overestimate + long sToN1 = BigInteger.unsignedLongPow(sLong, n - 1); + sLong = ((n - 1) * sLong + Long.divideUnsigned(x, sToN1)) / n; + + if (BigInteger.bitLengthForLong(sLong) * (n - 1) <= Long.SIZE) { + // Refine the estimate. + long u = sLong; + do { + sLong = u; + sToN1 = BigInteger.unsignedLongPow(sLong, n - 1); + u = ((n - 1) * sLong + Long.divideUnsigned(x, sToN1)) / n; + } while (u < sLong); // Terminate when non-decreasing. + + return new MutableBigInteger[] { + new MutableBigInteger(sLong), new MutableBigInteger(x - sToN1 * sLong) + }; + } + } + // s^(n - 1) could overflow long range, use MutableBigInteger loop instead + s = new MutableBigInteger(sLong); + } else { + final int rootLen = (bitLength - 1) / n + 1; // ⌈bitLength / n⌉ + int rootSh; + double rad = 0.0, approx = 0.0; + if (n < Double.PRECISION) { + // Set up the initial estimate of the iteration. + /* Since the following equality holds: + * nthRoot(x, n) == nthRoot(x/2^sh, n) * 2^(sh/n), + * + * to get an upper bound of the root of x, it suffices to find an integer sh + * and a real s such that s >= nthRoot(x/2^sh, n) and sh % n == 0. + * The upper bound will be s * 2^(sh/n), indeed: + * s * 2^(sh/n) >= nthRoot(x/2^sh, n) * 2^(sh/n) == nthRoot(x, n). + * To achieve this, we right shift the input of sh bits into finite double range, + * rounding up the result. + * + * The value of the shift sh is chosen in order to have the smallest number of + * trailing zeros in the double value of s after the significand (minimizing + * non-significant bits), to avoid losing bits in the significand. + */ + // Determine a right shift that is a multiple of n into finite double range. + rootSh = (bitLength - Double.PRECISION) / n; // rootSh < rootLen + /* Let x = this, P = Double.PRECISION, ME = Double.MAX_EXPONENT, + * bl = bitLength, sh = rootSh * n, ex = (bl - P) % n + * + * We have bl-sh = bl-((bl-P)-ex) = P + ex + * Since ex < n < P, we get P + ex ≤ ME, and so bl-sh ≤ ME. + * + * Recalling x < 2^bl: + * x >> sh < 2^(bl-sh) ≤ 2^ME < Double.MAX_VALUE + * Thus, rad ≤ 2^ME is in the range of finite doubles. + * + * Noting that ex ≥ 0, we get bl-sh = P + ex ≥ P + * which shows that x >> sh has at least P bits of precision, + * since bl-sh is its bit length. + */ + // Shift the value into finite double range + rad = this.toBigInteger().shiftRight(rootSh * n).doubleValue(); + + // Use the root of the shifted value as an estimate. + // rad ≤ 2^ME, so Math.nextUp(rad) < Double.MAX_VALUE + rad = Math.nextUp(rad); + approx = nthRootApprox(rad, n); + } else { // fp arithmetic gives too few correct bits + // Set the root shift to the root's bit length minus 1 + // The initial estimate will be 2^rootLen == 2 << (rootLen - 1) + rootSh = rootLen - 1; + } + + if (rootSh == 0) { + // approx has at most ⌈Double.PRECISION / n⌉ + 1 ≤ 19 integer bits + s = new MutableBigInteger((int) approx + 1); + } else { + // Allocate ⌈intLen / n⌉ ints to store the final root + s = new MutableBigInteger(new int[(intLen - 1) / n + 1]); + + if (n >= Double.PRECISION) { // fp arithmetic gives too few correct bits + // Set the initial estimate to 2 << (rootLen - 1) + s.value[0] = 2; + s.intLen = 1; + } else { + // Discard wrong integer bits from the initial estimate + // The reduced radicand rad has Math.getExponent(rad)+1 integer bits, but only + // the first Double.PRECISION leftmost bits are correct + // We scale the corresponding wrong bits of approx in the fraction part. + int wrongBits = ((Math.getExponent(rad) + 1) - Double.PRECISION) / n; + // Since rad <= 2^(bitLength - sh), then + // wrongBits <= ((bitLength - sh + 1) - Double.PRECISION) / n, + // so wrongBits is less than ⌈(bitLength - sh) / n⌉, + // the bit length of the exact shifted root, + // hence wrongBits + rootSh < ⌈(bitLength - sh) / n⌉ + rootSh == rootLen + rootSh += wrongBits; + approx = Math.scalb(approx, -wrongBits); + + // now approx has at most ⌈Double.PRECISION / n⌉ + 1 ≤ 19 integer bits + s.value[0] = (int) approx + 1; + s.intLen = 1; + } + + /* The Newton's recurrence roughly doubles the correct bits at each iteration. + * Instead of shifting the approximate root into the original range right now, + * we only double its bit length and then refine it with Newton's recurrence, + * using a suitable shifted radicand, in order to avoid computing and + * carrying trash bits in the approximate root. + * The shifted radicand is determined by the same reasoning used to get the + * initial estimate. + */ + // Refine the estimate to avoid computing non-significant bits + // rootSh is always less than rootLen, so correctBits >= 1 + for (int correctBits = rootLen - rootSh; correctBits < rootSh; correctBits <<= 1) { + s.leftShift(correctBits); + rootSh -= correctBits; + // Remove useless bits from the radicand + MutableBigInteger x = new MutableBigInteger(this); + x.rightShift(rootSh * n); + + newtonRecurrenceNthRoot(x, s, n, s.toBigInteger().pow(n - 1)); + s.add(ONE); // round up to ensure s is an upper bound of the root + } + + // Shift the approximate root back into the original range. + s.leftShift(rootSh); // Here rootSh > 0 always + } + } + + // Do the 1st iteration outside the loop to ensure an overestimate + newtonRecurrenceNthRoot(this, s, n, s.toBigInteger().pow(n - 1)); + // Refine the estimate. + do { + BigInteger sBig = s.toBigInteger(); + BigInteger sToN1 = sBig.pow(n - 1); + MutableBigInteger rem = new MutableBigInteger(sToN1.multiply(sBig).mag); + if (rem.subtract(this) <= 0) + return new MutableBigInteger[] { s, rem }; + + newtonRecurrenceNthRoot(this, s, n, sToN1); + } while (true); + } + + private static double nthRootApprox(double x, int n) { + return Math.nextUp(n == 3 ? Math.cbrt(x) : Math.pow(x, Math.nextUp(1.0 / n))); + } + + /** + * Computes {@code ((n-1)*s + x/sToN1)/n} and places the result in {@code s}. + */ + private static void newtonRecurrenceNthRoot( + MutableBigInteger x, MutableBigInteger s, int n, BigInteger sToN1) { + MutableBigInteger dividend = new MutableBigInteger(); + s.mul(n - 1, dividend); + MutableBigInteger xDivSToN1 = new MutableBigInteger(); + x.divide(new MutableBigInteger(sToN1.mag), xDivSToN1, false); + dividend.add(xDivSToN1); + dividend.divideOneWord(n, s); + } + /** * Calculate the integer square root {@code floor(sqrt(this))} and the remainder * if needed, where {@code sqrt(.)} denotes the mathematical square root. diff --git a/test/jdk/java/math/BigInteger/BigIntegerTest.java b/test/jdk/java/math/BigInteger/BigIntegerTest.java index 7da3fdac618..9018db6bb3c 100644 --- a/test/jdk/java/math/BigInteger/BigIntegerTest.java +++ b/test/jdk/java/math/BigInteger/BigIntegerTest.java @@ -1,5 +1,5 @@ /* - * Copyright (c) 1998, 2024, Oracle and/or its affiliates. All rights reserved. + * Copyright (c) 1998, 2025, Oracle and/or its affiliates. All rights reserved. * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. * * This code is free software; you can redistribute it and/or modify it @@ -26,7 +26,7 @@ * @library /test/lib * @build jdk.test.lib.RandomFactory * @run main BigIntegerTest - * @bug 4181191 4161971 4227146 4194389 4823171 4624738 4812225 4837946 4026465 8074460 8078672 8032027 8229845 + * @bug 4181191 4161971 4227146 4194389 4823171 4624738 4812225 4837946 4026465 8074460 8078672 8032027 8229845 8077587 * @summary tests methods in BigInteger (use -Dseed=X to set PRNG seed) * @run main/timeout=400 BigIntegerTest * @author madbot @@ -231,6 +231,9 @@ public class BigIntegerTest { if (!y.equals(z)) failCount1++; + + failCount1 += checkResult(x.signum() < 0 && power % 2 == 0 ? x.negate() : x, + y.nthRoot(power), "BigInteger.pow() inconsistent with BigInteger.nthRoot()"); } report("pow for " + order + " bits", failCount1); } @@ -412,6 +415,154 @@ public class BigIntegerTest { BigInteger.valueOf(x)).collect(Collectors.summingInt(g))); } + private static void nthRootSmall() { + int failCount = 0; + + // A non-positive degree should cause an exception. + int n = 0; + BigInteger x = BigInteger.ONE; + BigInteger s; + try { + s = x.nthRoot(n); + // If nthRoot() does not throw an exception that is a failure. + failCount++; + printErr("nthRoot() of non-positive degree did not throw an exception"); + } catch (ArithmeticException expected) { + // Not a failure + } + + // A negative value with even degree should cause an exception. + n = 4; + x = BigInteger.valueOf(-1); + try { + s = x.nthRoot(n); + // If nthRoot() does not throw an exception that is a failure. + failCount++; + printErr("nthRoot() of negative number and even degree did not throw an exception"); + } catch (ArithmeticException expected) { + // Not a failure + } + + // A negative value with odd degree should return -nthRoot(-x, n) + n = 3; + x = BigInteger.valueOf(-8); + failCount += checkResult(x.negate().nthRoot(n).negate(), x.nthRoot(n), + "nthRoot(" + x + ", " + n + ") != -nthRoot(" + x.negate() + ", " + n + ")"); + + // A zero value should return BigInteger.ZERO. + failCount += checkResult(BigInteger.ZERO, BigInteger.ZERO.nthRoot(n), + "nthRoot(0, " + n + ") != 0"); + + // A one degree should return x. + x = BigInteger.TWO; + failCount += checkResult(x, x.nthRoot(1), "nthRoot(" + x + ", 1) != " + x); + + n = 8; + // 1 <= value < 2^n should return BigInteger.ONE. + int end = 1 << n; + for (int i = 1; i < end; i++) { + failCount += checkResult(BigInteger.ONE, + BigInteger.valueOf(i).nthRoot(n), "nthRoot(" + i + ", " + n + ") != 1"); + } + + report("nthRootSmall", failCount); + } + + public static void nthRoot() { + nthRootSmall(); + + ToIntFunction f = (x) -> { + int n = random.nextInt(x.bitLength()) + 2; + int failCount = 0; + + // nth root of x^n -> x + BigInteger xN = x.pow(n); + failCount += checkResult(x, xN.nthRoot(n), "nthRoot() x^n -> x"); + + // nth root of x^n + 1 -> x + BigInteger xNup = xN.add(BigInteger.ONE); + failCount += checkResult(x, xNup.nthRoot(n), "nthRoot() x^n + 1 -> x"); + + // nth root of (x + 1)^n - 1 -> x + BigInteger up = + x.add(BigInteger.ONE).pow(n).subtract(BigInteger.ONE); + failCount += checkResult(x, up.nthRoot(n), "nthRoot() (x + 1)^n - 1 -> x"); + + // nthRoot(x, n)^n <= x + BigInteger r = x.nthRoot(n); + if (r.pow(n).compareTo(x) > 0) { + failCount++; + printErr("nthRoot(x, n)^n > x for x = " + x + ", n = " + n); + } + + // (nthRoot(x, n) + 1)^n > x + if (r.add(BigInteger.ONE).pow(n).compareTo(x) <= 0) { + failCount++; + printErr("(nthRoot(x, n) + 1)^n <= x for x = " + x + ", n = " + n); + } + + return failCount; + }; + + Stream.Builder sb = Stream.builder(); + int maxExponent = 256; + for (int i = 1; i <= maxExponent; i++) { + BigInteger p2 = BigInteger.ONE.shiftLeft(i); + sb.add(p2.subtract(BigInteger.ONE)); + sb.add(p2); + sb.add(p2.add(BigInteger.ONE)); + } + sb.add((new BigDecimal(Double.MAX_VALUE)).toBigInteger()); + sb.add((new BigDecimal(Double.MAX_VALUE)).toBigInteger().add(BigInteger.ONE)); + report("nthRoot for 2^N, 2^N - 1 and 2^N + 1, 1 <= N <= " + maxExponent, + sb.build().collect(Collectors.summingInt(f))); + + IntStream ints = random.ints(SIZE, 2, Integer.MAX_VALUE); + report("nthRoot for int", ints.mapToObj(x -> + BigInteger.valueOf(x)).collect(Collectors.summingInt(f))); + + LongStream longs = random.longs(SIZE, Integer.MAX_VALUE + 1L, Long.MAX_VALUE); + report("nthRoot for long", longs.mapToObj(x -> + BigInteger.valueOf(x)).collect(Collectors.summingInt(f))); + + DoubleStream doubles = random.doubles(SIZE, 0x1p63, Math.scalb(1.0, maxExponent)); + report("nthRoot for double", doubles.mapToObj(x -> + BigDecimal.valueOf(x).toBigInteger()).collect(Collectors.summingInt(f))); + } + + public static void nthRootAndRemainder() { + ToIntFunction g = (x) -> { + int failCount = 0; + int n = random.nextInt(x.bitLength()) + 2; + BigInteger xN = x.pow(n); + + // nth root of x^n -> x + BigInteger[] actual = xN.nthRootAndRemainder(n); + failCount += checkResult(x, actual[0], "nthRootAndRemainder()[0]"); + failCount += checkResult(BigInteger.ZERO, actual[1], "nthRootAndRemainder()[1]"); + + // nth root of x^n + 1 -> x + BigInteger xNup = xN.add(BigInteger.ONE); + actual = xNup.nthRootAndRemainder(n); + failCount += checkResult(x, actual[0], "nthRootAndRemainder()[0]"); + failCount += checkResult(BigInteger.ONE, actual[1], "nthRootAndRemainder()[1]"); + + // nth root of (x + 1)^n - 1 -> x + BigInteger up = + x.add(BigInteger.ONE).pow(n).subtract(BigInteger.ONE); + actual = up.nthRootAndRemainder(n); + failCount += checkResult(x, actual[0], "nthRootAndRemainder()[0]"); + BigInteger r = up.subtract(xN); + failCount += checkResult(r, actual[1], "nthRootAndRemainder()[1]"); + + return failCount; + }; + + IntStream bits = random.ints(SIZE, 3, Short.MAX_VALUE); + report("nthRootAndRemainder", bits.mapToObj(x -> + BigInteger.valueOf(x)).collect(Collectors.summingInt(g))); + } + public static void arithmetic(int order) { int failCount = 0; @@ -1294,6 +1445,8 @@ public class BigIntegerTest { squareRoot(); squareRootAndRemainder(); + nthRoot(); + nthRootAndRemainder(); bitCount(); bitLength();