8077587: BigInteger Roots

Reviewed-by: rgiulietti
This commit is contained in:
Fabio Romano 2025-09-08 16:10:22 +00:00 committed by Raffaello Giulietti
parent 6765a9d775
commit ab12fbfda2
3 changed files with 433 additions and 3 deletions

View File

@ -2742,6 +2742,85 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
return new BigInteger[] { sqrtRem[0].toBigInteger(), sqrtRem[1].toBigInteger() };
}
/**
* Returns the integer {@code n}<sup>th</sup> root of this BigInteger. The integer
* {@code n}<sup>th</sup> root {@code r} of the corresponding mathematical integer {@code x}
* is defined as follows:
* <ul>
* <li>if {@code x} &ge; 0, then {@code r} &ge; 0 is the largest integer such that
* {@code r}<sup>{@code n}</sup> &le; {@code x};
* <li>if {@code x} &lt; 0, then {@code r} &le; 0 is the smallest integer such that
* {@code r}<sup>{@code n}</sup> &ge; {@code x}.
* </ul>
* If the root is defined, it is equal to the value of
* {@code x.signum()}&sdot; &lfloor;{@code |nthRoot(x, n)|}&rfloor;,
* where {@code nthRoot(x, n)} denotes the real {@code n}<sup>th</sup> root of {@code x}
* treated as a real.
* Otherwise, the method throws an {@code ArithmeticException}.
*
* <p>Note that the magnitude of the integer {@code n}<sup>th</sup> root will be less than
* the magnitude of the real {@code n}<sup>th</sup> root if the latter is not representable
* as an integral value.
*
* @param n the root degree
* @return the integer {@code n}<sup>th</sup> 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}<sup>th</sup> root
* {@code r} of {@code this} and its remainder {@code this - r}<sup>{@code n}</sup>,
* respectively.
*
* @param n the root degree
* @return an array of two BigIntegers with the integer {@code n}<sup>th</sup> 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

View File

@ -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 <em>not</em> 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, <a href="https://maths-people.anu.edu.au/~brent/pd/mca-cup-0.5.9.pdf">
* Modern Computer Arithmetic</a>, 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.

View File

@ -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<BigInteger> 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<BigInteger> 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<BigInteger> 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();