8366478: BigDecimal roots

Reviewed-by: rgiulietti
This commit is contained in:
Fabio Romano 2026-05-07 23:59:07 +00:00 committed by Joe Darcy
parent f5aae68315
commit 8ae19e1ea4
2 changed files with 478 additions and 296 deletions

View File

@ -1,5 +1,5 @@
/*
* Copyright (c) 1996, 2025, Oracle and/or its affiliates. All rights reserved.
* Copyright (c) 1996, 2026, 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
@ -155,6 +155,7 @@ import jdk.internal.vm.annotation.Stable;
* <tr><th scope="row">Multiply</th><td>multiplier.scale() + multiplicand.scale()</td>
* <tr><th scope="row">Divide</th><td>dividend.scale() - divisor.scale()</td>
* <tr><th scope="row">Square root</th><td>ceil(radicand.scale()/2.0)</td>
* <tr><th scope="row">n<sup>th</sup> root</th><td>ceil((double) radicand.scale()/n)</td>
* </tbody>
* </table>
*
@ -2142,149 +2143,224 @@ public class BigDecimal extends Number implements Comparable<BigDecimal> {
* @since 9
*/
public BigDecimal sqrt(MathContext mc) {
return rootn(2, mc);
}
/**
* Returns an approximation to the {@code n}<sup>th</sup> root of {@code this}
* with rounding according to the context settings.
*
* <p>The preferred scale of the returned result is equal to
* {@code Math.ceilDiv(this.scale(), n)}. The value of the returned result is
* always within one ulp of the exact decimal value for the
* precision in question. If the rounding mode is {@link
* RoundingMode#HALF_UP HALF_UP}, {@link RoundingMode#HALF_DOWN
* HALF_DOWN}, or {@link RoundingMode#HALF_EVEN HALF_EVEN}, the
* result is within one half an ulp of the exact decimal value.
*
* <p>Special case:
* <ul>
* <li> The {@code n}<sup>th</sup> root of a number numerically equal to {@code
* ZERO} is numerically equal to {@code ZERO} with a preferred
* scale according to the general rule above. In particular, for
* {@code ZERO}, {@code ZERO.rootn(n, mc).equals(ZERO)} is true with
* any {@code MathContext} as an argument.
* </ul>
*
* @param n the root degree
* @param mc the context to use.
* @return the {@code n}<sup>th</sup> root of {@code this}.
* @throws ArithmeticException if {@code n == 0 || n == Integer.MIN_VALUE}.
* @throws ArithmeticException if {@code n} is even and {@code this} is negative.
* @throws ArithmeticException if {@code n} is negative and {@code this} is zero.
* @throws ArithmeticException if an exact result is requested
* ({@code mc.getPrecision() == 0}) and there is no finite decimal
* expansion of the exact result
* @throws ArithmeticException if
* {@code (mc.getRoundingMode() == RoundingMode.UNNECESSARY}) and
* the exact result cannot fit in {@code mc.getPrecision()} digits.
* @see #sqrt(MathContext)
* @see BigInteger#rootn(int)
* @since 27
* @apiNote Note that calling {@code rootn(2, mc)} is equivalent to calling {@code sqrt(mc)}.
*/
public BigDecimal rootn(int n, MathContext mc) {
// Special cases
if (n == 0)
throw new ArithmeticException("Zero root degree");
final int signum = signum();
if (signum != 1) {
switch (signum) {
case -1 -> throw new ArithmeticException("Attempted square root of negative BigDecimal");
case 0 -> {
BigDecimal result = valueOf(0L, scale/2);
assert squareRootResultAssertions(result, mc);
return result;
}
default -> throw new AssertionError("Bad value from signum");
}
if (signum < 0 && (n & 1) == 0)
throw new ArithmeticException("Negative radicand with even root degree");
final int preferredScale = saturateLong(Math.ceilDiv((long) this.scale, n));
if (signum == 0) {
if (n < 0)
throw new ArithmeticException("Zero radicand with negative root degree");
return zeroValueOf(preferredScale);
}
/*
* The main steps of the algorithm below are as follows,
* first argument reduce the value to an integer
* using the following relations:
* using the following relations, assuming n > 0:
*
* x = y * 10 ^ exp
* sqrt(x) = sqrt(y) * 10^(exp / 2) if exp is even
* sqrt(x) = sqrt(y*10) * 10^((exp-1)/2) is exp is odd
* rootn(x, n) = rootn(y, n) * 10^(exp / n) if exp mod n == 0
* rootn(x, n) = rootn(y*10^(exp mod n), n) * 10^((exp - (exp mod n))/n) otherwise
*
* Then use BigInteger.sqrt() on the reduced value to compute
* where exp mod n == Math.floorMod(exp, n).
*
* Then use BigInteger.rootn() on the reduced value to compute
* the numerical digits of the desired result.
*
* Finally, scale back to the desired exponent range and
* perform any adjustment to get the preferred scale in the
* representation.
*/
// The code below favors relative simplicity over checking
// for special cases that could run faster.
final int preferredScale = Math.ceilDiv(this.scale, 2);
final int nAbs = Math.absExact(n);
BigDecimal result;
if (mc.roundingMode == RoundingMode.UNNECESSARY || mc.precision == 0) { // Exact result requested
// To avoid trailing zeros in the result, strip trailing zeros.
final BigDecimal stripped = this.stripTrailingZeros();
final int strippedScale = stripped.scale;
if ((strippedScale & 1) != 0) // 10*stripped.unscaledValue() can't be an exact square
throw new ArithmeticException("Computed square root not exact.");
// Check for even powers of 10. Numerically sqrt(10^2N) = 10^N
if (stripped.isPowerOfTen()) {
result = valueOf(1L, strippedScale >> 1);
// Adjust to requested precision and preferred
// scale as appropriate.
return result.adjustToPreferredScale(preferredScale, mc.precision);
}
// if stripped.scale is not a multiple of n,
// 10^((-stripped.scale) mod n)*stripped.unscaledValue() can't be an exact power
if (stripped.scale % n != 0)
throw new ArithmeticException("Computed root not exact.");
// After stripTrailingZeros, the representation is normalized as
//
// unscaledValue * 10^(-scale)
//
// where unscaledValue is an integer with the minimum
// precision for the cohort of the numerical value and the scale is even.
BigInteger[] sqrtRem = stripped.unscaledValue().sqrtAndRemainder();
result = new BigDecimal(sqrtRem[0], strippedScale >> 1);
// precision for the cohort of the numerical value and the scale is a multiple of n.
BigInteger[] rootRem = stripped.unscaledValue().rootnAndRemainder(nAbs);
result = new BigDecimal(rootRem[0], stripped.scale / nAbs);
// If result^nAbs != this numerically, the root isn't exact
if (rootRem[1].signum != 0)
throw new ArithmeticException("Computed root not exact.");
// If result*result != this numerically or requires too high precision,
// the square root isn't exact
if (sqrtRem[1].signum != 0 || mc.precision != 0 && result.precision() > mc.precision)
throw new ArithmeticException("Computed square root not exact.");
// Test numerical properties at full precision before any
// scale adjustments.
assert squareRootResultAssertions(result, mc);
// Adjust to requested precision and preferred
// scale as appropriate.
if (n > 0) {
// If result requires too high precision, the root isn't exact
if (mc.precision != 0 && result.precision() > mc.precision)
throw new ArithmeticException("Computed root not exact.");
} else {
try {
result = ONE.divide(result, mc);
} catch (ArithmeticException e) {
// The exact result requires too high precision,
// including non-terminating decimal expansions
throw new ArithmeticException("Computed root not exact.");
}
}
// Test numerical properties at full precision before any scale adjustments.
assert rootnResultAssertions(result, mc, n);
// Adjust to requested precision and preferred scale as appropriate.
return result.adjustToPreferredScale(preferredScale, mc.precision);
}
// To allow BigInteger.sqrt() to be used to get the square
// root, it is necessary to normalize the input so that
// its integer part is sufficient to get the square root
// Handle negative radicands
BigDecimal x = this;
if (signum < 0) {
x = x.negate();
if (mc.roundingMode == RoundingMode.FLOOR) {
mc = new MathContext(mc.precision, RoundingMode.UP);
} else if (mc.roundingMode == RoundingMode.CEILING) {
mc = new MathContext(mc.precision, RoundingMode.DOWN);
}
}
// To allow BigInteger.rootn() to be used to get the root,
// it is necessary to normalize the input so that
// its integer part is sufficient to get the root
// with the desired precision.
final boolean halfWay = isHalfWay(mc.roundingMode);
// To obtain a square root with N digits,
// the radicand must have at least 2*(N-1)+1 == 2*N-1 digits.
final long minWorkingPrec = ((mc.precision + (halfWay ? 1L : 0L)) << 1) - 1L;
// normScale is the number of digits to take from the fraction of the input
long normScale = minWorkingPrec - this.precision() + this.scale;
normScale += normScale & 1L; // the scale for normalizing must be even
final long rootDigits = mc.precision + (halfWay ? 1L : 0L);
// To obtain an n-th root with k digits,
// the radicand must have at least n*(k-1)+1 digits.
final long minWorkingPrec = nAbs * (rootDigits - 1L) + 1L;
final long workingScale = this.scale - normScale;
if (workingScale != (int) workingScale)
throw new ArithmeticException("Overflow");
long normScale; // the number of digits to take from the fraction of the input
BigDecimal working = null, xInv = null;
BigInteger workingInt;
if (n > 0) {
normScale = minWorkingPrec - x.precision() + x.scale;
int mod = Math.floorMod(normScale, n);
if (mod != 0) // the scale for normalizing must be a multiple of n
normScale += n - mod;
BigDecimal working = new BigDecimal(this.intVal, this.intCompact, (int) workingScale, this.precision);
BigInteger workingInt = working.toBigInteger();
working = new BigDecimal(x.intVal, x.intCompact, checkScaleNonZero(x.scale - normScale), x.precision);
workingInt = working.toBigInteger();
} else { // Handle negative degrees
/* Computing the n-th root of x is equivalent
* to computing the (-n)-th root of 1/x.
*/
// Compute the scale for xInv, in order to ensure
// that xInv's precision is at least minWorkingPrec
final int fracZeros = x.precision() - 1 - (x.isPowerOfTen() ? 1 : 0);
normScale = minWorkingPrec + fracZeros - x.scale;
int mod = Math.floorMod(normScale, nAbs);
if (mod != 0)
normScale += nAbs - mod;
BigInteger sqrt;
long resultScale = normScale >> 1;
// Round sqrt with the specified settings
xInv = ONE.divide(x, checkScaleNonZero(normScale), RoundingMode.DOWN);
workingInt = xInv.unscaledValue();
}
// Compute and round the root with the specified settings
BigInteger root;
long resultScale = normScale / nAbs;
boolean increment = false;
if (halfWay) { // half-way rounding
BigInteger workingSqrt = workingInt.sqrt();
BigInteger[] rootRem = workingInt.rootnAndRemainder(nAbs);
// remove the one-tenth digit
BigInteger[] quotRem10 = workingSqrt.divideAndRemainder(BigInteger.TEN);
sqrt = quotRem10[0];
BigInteger[] quotRem10 = rootRem[0].divideAndRemainder(BigInteger.TEN);
root = quotRem10[0];
resultScale--;
boolean increment = false;
int digit = quotRem10[1].intValue();
if (digit > 5) {
increment = true;
} else if (digit == 5) {
if (mc.roundingMode == RoundingMode.HALF_UP
|| mc.roundingMode == RoundingMode.HALF_EVEN && sqrt.testBit(0)
|| mc.roundingMode == RoundingMode.HALF_EVEN && root.testBit(0)
// Check if remainder is non-zero
|| !workingInt.equals(workingSqrt.multiply(workingSqrt))
|| !working.isInteger()) {
|| rootRem[1].signum != 0
|| (n > 0 ? !working.isInteger() : xInv.multiply(x).compareMagnitude(ONE) != 0)) {
increment = true;
}
}
if (increment)
sqrt = sqrt.add(1L);
} else {
switch (mc.roundingMode) {
case DOWN, FLOOR -> sqrt = workingInt.sqrt(); // No need to round
case DOWN, FLOOR -> root = workingInt.rootn(nAbs); // No need to round
case UP, CEILING -> {
BigInteger[] sqrtRem = workingInt.sqrtAndRemainder();
sqrt = sqrtRem[0];
BigInteger[] rootRem = workingInt.rootnAndRemainder(nAbs);
root = rootRem[0];
// Check if remainder is non-zero
if (sqrtRem[1].signum != 0 || !working.isInteger())
sqrt = sqrt.add(1L);
if (rootRem[1].signum != 0
|| (n > 0 ? !working.isInteger() : xInv.multiply(x).compareMagnitude(ONE) != 0))
increment = true;
}
default -> throw new AssertionError("Unexpected value for RoundingMode: " + mc.roundingMode);
}
}
if (increment) {
root = root.add(1L);
}
result = new BigDecimal(sqrt, checkScale(sqrt, resultScale), mc); // mc ensures no increase of precision
// Test numerical properties at full precision before any
// scale adjustments.
assert squareRootResultAssertions(result, mc);
// Adjust to requested precision and preferred
// scale as appropriate.
if (result.scale > preferredScale) // else can't increase the result's precision to fit the preferred scale
result = new BigDecimal(root, checkScale(root, resultScale), mc); // mc ensures no increase of precision
// Test numerical properties at full precision before any scale adjustments.
assert rootnResultAssertions(result, mc, n);
// Adjust to requested precision and preferred scale as appropriate.
if (result.scale > preferredScale) // else can't increase result's precision to fit the preferred scale
result = stripZerosToMatchScale(result.intVal, result.intCompact, result.scale, preferredScale);
return result;
return signum > 0 ? result : result.negate();
}
/**
@ -2315,12 +2391,8 @@ public class BigDecimal extends Number implements Comparable<BigDecimal> {
};
}
private BigDecimal square() {
return this.multiply(this);
}
private boolean isPowerOfTen() {
return BigInteger.ONE.equals(this.unscaledValue());
return this.stripTrailingZeros().unscaledValue().equals(BigInteger.ONE);
}
/**
@ -2331,92 +2403,102 @@ public class BigDecimal extends Number implements Comparable<BigDecimal> {
*
* <ul>
*
* <li> For DOWN and FLOOR, result^2 must be {@code <=} the input
* and (result+ulp)^2 must be {@code >} the input.
* <li> For DOWN and FLOOR if input > 0 and CEIL if input < 0,
* |result|^n must be {@code <=} |input|
* and (|result|+ulp)^n must be {@code >} |input|.
*
* <li>Conversely, for UP and CEIL, result^2 must be {@code >=}
* the input and (result-ulp)^2 must be {@code <} the input.
* <li>Conversely, for UP and FLOOR if input < 0 and CEIL if input > 0,
* |result|^n must be {@code >=} |input|
* and (|result|-ulp)^n must be {@code <} |input|.
* </ul>
*/
private boolean squareRootResultAssertions(BigDecimal result, MathContext mc) {
if (result.signum() == 0) {
return squareRootZeroResultAssertions(result, mc);
} else {
RoundingMode rm = mc.getRoundingMode();
BigDecimal ulp = result.ulp();
BigDecimal neighborUp = result.add(ulp);
// Make neighbor down accurate even for powers of ten
if (result.isPowerOfTen()) {
ulp = ulp.divide(TEN);
}
BigDecimal neighborDown = result.subtract(ulp);
// Both the starting value and result should be nonzero and positive.
assert (result.signum() == 1 &&
this.signum() == 1) :
"Bad signum of this and/or its sqrt.";
switch (rm) {
case DOWN:
case FLOOR:
assert
result.square().compareTo(this) <= 0 &&
neighborUp.square().compareTo(this) > 0:
"Square of result out for bounds rounding " + rm;
return true;
case UP:
case CEILING:
assert
result.square().compareTo(this) >= 0 &&
neighborDown.square().compareTo(this) < 0:
"Square of result out for bounds rounding " + rm;
return true;
case HALF_DOWN:
case HALF_EVEN:
case HALF_UP:
BigDecimal err = result.square().subtract(this).abs();
BigDecimal errUp = neighborUp.square().subtract(this);
BigDecimal errDown = this.subtract(neighborDown.square());
// All error values should be positive so don't need to
// compare absolute values.
int err_comp_errUp = err.compareTo(errUp);
int err_comp_errDown = err.compareTo(errDown);
assert
errUp.signum() == 1 &&
errDown.signum() == 1 :
"Errors of neighbors squared don't have correct signs";
// For breaking a half-way tie, the return value may
// have a larger error than one of the neighbors. For
// example, the square root of 2.25 to a precision of
// 1 digit is either 1 or 2 depending on how the exact
// value of 1.5 is rounded. If 2 is returned, it will
// have a larger rounding error than its neighbor 1.
assert
err_comp_errUp <= 0 ||
err_comp_errDown <= 0 :
"Computed square root has larger error than neighbors for " + rm;
assert
((err_comp_errUp == 0 ) ? err_comp_errDown < 0 : true) &&
((err_comp_errDown == 0 ) ? err_comp_errUp < 0 : true) :
"Incorrect error relationships";
// && could check for digit conditions for ties too
return true;
default: // Definition of UNNECESSARY already verified.
return true;
private boolean rootnResultAssertions(BigDecimal result, MathContext mc, int n) {
BigDecimal rad = this.abs(), resAbs = result.abs();
RoundingMode rm = mc.roundingMode;
if (this.signum() < 0) {
if (rm == RoundingMode.FLOOR) {
rm = RoundingMode.UP;
} else if (rm == RoundingMode.CEILING) {
rm = RoundingMode.DOWN;
}
}
}
private boolean squareRootZeroResultAssertions(BigDecimal result, MathContext mc) {
return this.compareTo(ZERO) == 0;
int nAbs = Math.abs(n);
BigDecimal ulp = resAbs.ulp();
BigDecimal neighborUp = resAbs.add(ulp);
// Make neighbor down accurate even for powers of ten
if (resAbs.isPowerOfTen()) {
ulp = ulp.scaleByPowerOfTen(-1);
}
BigDecimal neighborDown = resAbs.subtract(ulp);
switch (rm) {
case DOWN:
case FLOOR:
assert
(n > 0 ? resAbs.pow(nAbs).compareTo(rad) <= 0 &&
neighborUp.pow(nAbs).compareTo(rad) > 0
: resAbs.pow(nAbs).multiply(rad).compareTo(ONE) <= 0 &&
neighborUp.pow(nAbs).multiply(rad).compareTo(ONE) > 0)
: "Power of result out for bounds rounding " + rm;
return true;
case UP:
case CEILING:
assert
(n > 0 ? resAbs.pow(nAbs).compareTo(rad) >= 0 &&
neighborDown.pow(nAbs).compareTo(rad) < 0
: resAbs.pow(nAbs).multiply(rad).compareTo(ONE) >= 0 &&
neighborDown.pow(nAbs).multiply(rad).compareTo(ONE) < 0)
: "Power of result out for bounds rounding " + rm;
return true;
case HALF_DOWN:
case HALF_EVEN:
case HALF_UP:
BigDecimal err, errUp, errDown;
if (n > 0) {
err = resAbs.pow(nAbs).subtract(rad).abs();
errUp = neighborUp.pow(nAbs).subtract(rad);
errDown = rad.subtract(neighborDown.pow(nAbs));
} else {
err = resAbs.pow(nAbs).multiply(rad).subtract(ONE).abs();
errUp = neighborUp.pow(nAbs).multiply(rad).subtract(ONE);
errDown = ONE.subtract(neighborDown.pow(nAbs).multiply(rad));
}
// All error values should be positive
// so don't need to compare absolute values.
int err_comp_errUp = err.compareTo(errUp);
int err_comp_errDown = err.compareTo(errDown);
assert
errUp.signum() == 1 &&
errDown.signum() == 1
: "Errors of neighbors powered don't have correct signs";
// For breaking a half-way tie, the return value may
// have a larger error than one of the neighbors. For
// example, the square root of 2.25 to a precision of
// 1 digit is either 1 or 2 depending on how the exact
// value of 1.5 is rounded. If 2 is returned, it will
// have a larger rounding error than its neighbor 1.
assert
err_comp_errUp <= 0 ||
err_comp_errDown <= 0 :
"Computed root has larger error than neighbors for " + rm;
assert
((err_comp_errUp == 0 ) ? err_comp_errDown < 0 : true) &&
((err_comp_errDown == 0 ) ? err_comp_errUp < 0 : true) :
"Incorrect error relationships";
// && could check for digit conditions for ties too
return true;
default: // Definition of UNNECESSARY already verified.
return true;
}
}
/**

View File

@ -1,5 +1,5 @@
/*
* Copyright (c) 2016, 2025, Oracle and/or its affiliates. All rights reserved.
* Copyright (c) 2016, 2026, 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
@ -31,7 +31,8 @@ import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.List;
import java.util.*;
import java.util.stream.IntStream;
import static java.math.BigDecimal.ONE;
import static java.math.BigDecimal.TWO;
@ -48,13 +49,14 @@ public class SquareRootTests {
public static void main(String... args) {
int failures = 0;
failures += negativeTests();
failures += zerothRootTests();
failures += negativeWithEvenDegreeTests();
failures += zeroTests();
failures += oneDigitTests();
failures += twoDigitTests();
failures += evenPowersOfTenTests();
failures += squareRootTwoTests();
failures += lowPrecisionPerfectSquares();
failures += rootTwoTests();
failures += lowPrecisionPerfectPowers();
failures += almostFourRoundingDown();
failures += almostFourRoundingUp();
failures += nearTen();
@ -69,19 +71,41 @@ public class SquareRootTests {
}
}
private static int negativeTests() {
private static int zerothRootTests() {
int failures = 0;
for (long i = -5; i < 5; i++) {
for (int j = -5; j < 5; j++) {
try {
BigDecimal input = BigDecimal.valueOf(i, j);
BigDecimal result = input.rootn(0, MathContext.DECIMAL64);
System.err.println("Unexpected 0th root: (" +
input + ").rootn(0) = " + result );
failures += 1;
} catch (ArithmeticException e) {
; // Expected
}
}
}
return failures;
}
private static int negativeWithEvenDegreeTests() {
int failures = 0;
for (long i = -10; i < 0; i++) {
for (int j = -5; j < 5; j++) {
try {
BigDecimal input = BigDecimal.valueOf(i, j);
BigDecimal result = input.sqrt(MathContext.DECIMAL64);
System.err.println("Unexpected sqrt of negative: (" +
input + ").sqrt() = " + result );
failures += 1;
} catch (ArithmeticException e) {
; // Expected
BigDecimal input = BigDecimal.valueOf(i, j);
for (int n = -4; n <= 4; n += 2) {
try {
BigDecimal result = input.rootn(n, MathContext.DECIMAL64);
System.err.println("Unexpected nth root of negative: (" +
input + ").rootn(" + n + ") = " + result );
failures += 1;
} catch (ArithmeticException e) {
; // Expected
}
}
}
}
@ -93,43 +117,72 @@ public class SquareRootTests {
int failures = 0;
for (int i = -100; i < 100; i++) {
BigDecimal expected = BigDecimal.valueOf(0L, i/2);
// These results are independent of rounding mode
failures += compare(BigDecimal.valueOf(0L, i).sqrt(MathContext.UNLIMITED),
expected, true, "zeros");
BigDecimal input = BigDecimal.valueOf(0L, i);
for (int n = -10; n <= 0; n++) {
try {
BigDecimal result = input.rootn(n, MathContext.DECIMAL64);
System.err.println("Unexpected nth root of zero: (" +
input + ").rootn(" + n + ") = " + result );
failures += 1;
} catch (ArithmeticException e) {
; // Expected
}
}
failures += compare(BigDecimal.valueOf(0L, i).sqrt(MathContext.DECIMAL64),
expected, true, "zeros");
for (int n = 1; n < 10; n++) {
BigDecimal expected = BigDecimal.valueOf(0L, Math.ceilDiv(i, n));
// These results are independent of rounding mode
failures += compare(input.rootn(n, MathContext.UNLIMITED),
expected, true, "zeros");
failures += compare(input.rootn(n, MathContext.DECIMAL64),
expected, true, "zeros");
}
}
return failures;
}
private static RoundingMode positiveRoundingMode(RoundingMode rm) {
return rm == RoundingMode.FLOOR ? RoundingMode.UP :
rm == RoundingMode.CEILING ? RoundingMode.DOWN : rm;
}
private static List<RoundingMode> roundingModes() {
List<RoundingMode> modes = new ArrayList<>(Arrays.asList(RoundingMode.values()));
modes.remove(RoundingMode.UNNECESSARY);
return modes;
}
/**
* Probe inputs with one digit of precision, 1 ... 9 and those
* values scaled by 10^-1, 0.1, ... 0.9.
* Probe inputs with one digit of precision, ±1 ... ±9 and those
* values scaled by 10^-1, ±0.1 ... ±0.9.
*/
private static int oneDigitTests() {
int failures = 0;
List<BigDecimal> oneToNine =
List.of(ONE, TWO, valueOf(3),
valueOf(4), valueOf(5), valueOf(6),
valueOf(7), valueOf(8), valueOf(9));
List<RoundingMode> modes =
List.of(RoundingMode.UP, RoundingMode.DOWN,
RoundingMode.CEILING, RoundingMode.FLOOR,
RoundingMode.HALF_UP, RoundingMode.HALF_DOWN, RoundingMode.HALF_EVEN);
List<BigDecimal> oneToNine = IntStream.rangeClosed(1, 9).mapToObj(BigDecimal::valueOf).toList();
List<RoundingMode> modes = roundingModes();
for (int i = 1; i < 20; i++) {
for (RoundingMode rm : modes) {
MathContext mc = new MathContext(i, rm);
MathContext positiveMC = new MathContext(i, positiveRoundingMode(rm));
for (BigDecimal bd : oneToNine) {
MathContext mc = new MathContext(i, rm);
failures += compareSqrtImplementations(bd, mc);
BigDecimal minus_bd = bd.negate();
for (int n = 1; n < 11; n += 2) {
failures += compare(minus_bd.rootn( n, mc), bd.rootn( n, positiveMC).negate(), true, "one digit");
failures += compare(minus_bd.rootn(-n, mc), bd.rootn(-n, positiveMC).negate(), true, "one digit");
}
bd = bd.scaleByPowerOfTen(-1);
failures += compareSqrtImplementations(bd, mc);
bd = bd.multiply(ONE_TENTH);
failures += compareSqrtImplementations(bd, mc);
minus_bd = bd.negate();
for (int n = 1; n < 11; n += 2) {
failures += compare(minus_bd.rootn( n, mc), bd.rootn( n, positiveMC).negate(), true, "one digit");
failures += compare(minus_bd.rootn(-n, mc), bd.rootn(-n, positiveMC).negate(), true, "one digit");
}
}
}
}
@ -138,28 +191,32 @@ public class SquareRootTests {
}
/**
* Probe inputs with two digits of precision, (10 ... 99) and
* those values scaled by 10^-1 (1, ... 9.9) and scaled by 10^-2
* (0.1 ... 0.99).
* Probe inputs with two digits of precision, (±10 ... ±99) and
* those values scaled by 10^-1 (±1 ... ±9.9) and scaled by 10^-2
* (±0.1 ... ±0.99).
*/
private static int twoDigitTests() {
int failures = 0;
List<RoundingMode> modes =
List.of(RoundingMode.UP, RoundingMode.DOWN,
RoundingMode.CEILING, RoundingMode.FLOOR,
RoundingMode.HALF_UP, RoundingMode.HALF_DOWN, RoundingMode.HALF_EVEN);
List<RoundingMode> modes = roundingModes();
for (int i = 10; i < 100; i++) {
BigDecimal bd0 = BigDecimal.valueOf(i);
BigDecimal bd1 = bd0.multiply(ONE_TENTH);
BigDecimal bd2 = bd1.multiply(ONE_TENTH);
BigDecimal bd1 = bd0.scaleByPowerOfTen(-1);
BigDecimal bd2 = bd1.scaleByPowerOfTen(-1);
for (BigDecimal bd : List.of(bd0, bd1, bd2)) {
for (int precision = 1; i < 20; i++) {
for (int prec = 1; prec < 20; prec++) {
for (RoundingMode rm : modes) {
MathContext mc = new MathContext(precision, rm);
MathContext mc = new MathContext(prec, rm);
MathContext positiveMC = new MathContext(prec, positiveRoundingMode(rm));
failures += compareSqrtImplementations(bd, mc);
BigDecimal minus_bd = bd.negate();
for (int n = 1; n < 11; n += 2) {
failures += compare(minus_bd.rootn( n, mc), bd.rootn( n, positiveMC).negate(), true, "two digits");
failures += compare(minus_bd.rootn(-n, mc), bd.rootn(-n, positiveMC).negate(), true, "two digits");
}
}
}
}
@ -178,19 +235,24 @@ public class SquareRootTests {
MathContext unnecessary = new MathContext(1, RoundingMode.UNNECESSARY);
MathContext arbitrary = new MathContext(0, RoundingMode.CEILING);
BigDecimal[] errCases = {
// (strippedScale & 1) != 0
BigDecimal.TEN,
// (strippedScale & 1) == 0 && !stripped.isPowerOfTen() && sqrtRem[1].signum != 0
BigDecimal.TWO,
Object[][] errCases = {
// strippedScale % n != 0
{ BigDecimal.TEN, 2 },
// strippedScale % n == 0 && sqrtRem[1].signum != 0
{ BigDecimal.TWO, 2 },
// sqrtRem[1].signum == 0 && n < 0
{ BigDecimal.valueOf(9L), -2 },
};
for (BigDecimal input : errCases) {
for (Object[] errCase : errCases) {
BigDecimal input = (BigDecimal) errCase[0];
int n = (int) errCase[1];
BigDecimal result;
// mc.roundingMode == RoundingMode.UNNECESSARY
try {
result = input.sqrt(unnecessary);
System.err.println("Unexpected sqrt with UNNECESSARY RoundingMode: (" + input + ").sqrt() = " + result);
result = input.rootn(n, unnecessary);
System.err.println("Unexpected nth root with UNNECESSARY RoundingMode: ("
+ input + ").rootn(" + n + ") = " + result);
failures += 1;
} catch (ArithmeticException e) {
// Expected
@ -198,17 +260,18 @@ public class SquareRootTests {
// mc.roundingMode != RoundingMode.UNNECESSARY && mc.precision == 0
try {
result = input.sqrt(arbitrary);
System.err.println("Unexpected sqrt with mc.precision == 0: (" + input + ").sqrt() = " + result);
result = input.rootn(n, arbitrary);
System.err.println("Unexpected nth root with mc.precision == 0: ("
+ input + ").rootn(" + n + ") = " + result);
failures += 1;
} catch (ArithmeticException e) {
// Expected
}
}
// (strippedScale & 1) == 0
// strippedScale % n == 0
// !stripped.isPowerOfTen() && sqrtRem[1].signum == 0 && (mc.precision != 0 && result.precision() > mc.precision)
// sqrtRem[1].signum == 0 && n > 0 && (mc.precision != 0 && result.precision() > mc.precision)
try {
BigDecimal input = BigDecimal.valueOf(121);
BigDecimal result = input.sqrt(unnecessary);
@ -219,26 +282,27 @@ public class SquareRootTests {
// Expected
}
BigDecimal four = BigDecimal.valueOf(4);
BigDecimal four = BigDecimal.valueOf(4), oneHalf = BigDecimal.valueOf(5, 1);
Object[][] cases = {
// stripped.isPowerOfTen() && mc.roundingMode == RoundingMode.UNNECESSARY
{ BigDecimal.ONE, unnecessary, BigDecimal.ONE },
// stripped.isPowerOfTen() && mc.roundingMode != RoundingMode.UNNECESSARY && mc.precision == 0
{ BigDecimal.ONE, arbitrary, BigDecimal.ONE },
// !stripped.isPowerOfTen() && mc.roundingMode == RoundingMode.UNNECESSARY
// && sqrtRem[1].signum == 0 && mc.precision == 0
{ four, new MathContext(0, RoundingMode.UNNECESSARY), BigDecimal.TWO },
// !stripped.isPowerOfTen() && mc.roundingMode != RoundingMode.UNNECESSARY
// && sqrtRem[1].signum == 0 && mc.precision == 0
{ four, arbitrary, BigDecimal.TWO },
// !stripped.isPowerOfTen() && sqrtRem[1].signum == 0
// mc.roundingMode == RoundingMode.UNNECESSARY && sqrtRem[1].signum == 0 && n > 0
// && mc.precision == 0
{ four, 2, new MathContext(0, RoundingMode.UNNECESSARY), BigDecimal.TWO },
// mc.roundingMode != RoundingMode.UNNECESSARY && sqrtRem[1].signum == 0 && n > 0
// && mc.precision == 0
{ four, 2, arbitrary, BigDecimal.TWO },
// sqrtRem[1].signum == 0 && n > 0
// && (mc.precision != 0 && result.precision() <= mc.precision)
{ four, unnecessary, BigDecimal.TWO },
{ four, 2, unnecessary, BigDecimal.TWO },
// mc.roundingMode == RoundingMode.UNNECESSARY && sqrtRem[1].signum == 0 && n < 0
{ four, -2, unnecessary, oneHalf },
// mc.roundingMode != RoundingMode.UNNECESSARY && sqrtRem[1].signum == 0 && n < 0
// && mc.precision == 0
{ four, -2, arbitrary, oneHalf },
};
for (Object[] testCase : cases) {
BigDecimal expected = (BigDecimal) testCase[2];
BigDecimal result = ((BigDecimal) testCase[0]).sqrt((MathContext) testCase[1]);
BigDecimal expected = (BigDecimal) testCase[3];
BigDecimal result = ((BigDecimal) testCase[0]).rootn((int) testCase[1], (MathContext) testCase[2]);
failures += compare(expected, result, true, "Exact results");
}
@ -294,74 +358,69 @@ public class SquareRootTests {
return failures;
}
private static int squareRootTwoTests() {
private static int rootTwoTests() {
int failures = 0;
// Square root of 2 truncated to 65 digits
BigDecimal highPrecisionRoot2 =
BigDecimal highPrecisionSqrt2 =
new BigDecimal("1.41421356237309504880168872420969807856967187537694807317667973799");
// Cube root of 2 truncated to 65 digits
BigDecimal highPrecisionCbrt2 =
new BigDecimal("1.25992104989487316476721060727822835057025146470150798008197511215");
RoundingMode[] modes = {
RoundingMode.UP, RoundingMode.DOWN,
RoundingMode.CEILING, RoundingMode.FLOOR,
RoundingMode.HALF_UP, RoundingMode.HALF_DOWN, RoundingMode.HALF_EVEN
};
List<RoundingMode> modes = roundingModes();
// For each interesting rounding mode, for precisions 1 to, say,
// 63 numerically compare TWO.sqrt(mc) to
// highPrecisionRoot2.round(mc) and the alternative internal high-precision
// implementation of square root.
// highPrecisionSqrt2.round(mc) and the alternative internal high-precision
// implementation of square root, and compare TWO.rootn(3, mc) to
// highPrecisionCbrt2.round(mc).
for (RoundingMode mode : modes) {
for (int precision = 1; precision < 63; precision++) {
MathContext mc = new MathContext(precision, mode);
BigDecimal expected = highPrecisionRoot2.round(mc);
BigDecimal computed = TWO.sqrt(mc);
BigDecimal altComputed = BigSquareRoot.sqrt(TWO, mc);
BigDecimal expectedSqrt = highPrecisionSqrt2.round(mc);
BigDecimal computedSqrt = TWO.sqrt(mc);
BigDecimal altComputedSqrt = BigSquareRoot.sqrt(TWO, mc);
failures += equalNumerically(expected, computed, "sqrt(2)");
failures += equalNumerically(computed, altComputed, "computed & altComputed");
failures += equalNumerically(expectedSqrt, computedSqrt, "sqrt(2)");
failures += equalNumerically(computedSqrt, altComputedSqrt, "computed & altComputed");
BigDecimal expectedCbrt = highPrecisionCbrt2.round(mc);
BigDecimal computedCbrt = TWO.rootn(3, mc);
failures += equalNumerically(expectedCbrt, computedCbrt, "rootn(2, 3)");
}
}
return failures;
}
private static int lowPrecisionPerfectSquares() {
private static int lowPrecisionPerfectPowers() {
int failures = 0;
// For 5^2 through 9^2, if the input is rounded to one digit
// first before the root is computed, the wrong answer will
// result. Verify results and scale for different rounding
// modes and precisions.
long[][] squaresWithOneDigitRoot = {{ 4, 2},
{ 9, 3},
{25, 5},
{36, 6},
{49, 7},
{64, 8},
{81, 9}};
for (int i = -9; i <= 9; i++) {
for (int n = 1; n < 10; n++) {
BigDecimal expected = BigDecimal.valueOf(n % 2 != 0 ? i : Math.abs(i));
BigDecimal pow = BigDecimal.valueOf(Math.powExact(i, n));
for (long[] squareAndRoot : squaresWithOneDigitRoot) {
BigDecimal square = new BigDecimal(squareAndRoot[0]);
BigDecimal expected = new BigDecimal(squareAndRoot[1]);
for (int scale = 0; scale <= 4; scale++) {
BigDecimal scaledSquare = square.setScale(scale, RoundingMode.UNNECESSARY);
int expectedScale = Math.ceilDiv(scale, 2);
for (int precision = 0; precision <= 5; precision++) {
for (RoundingMode rm : RoundingMode.values()) {
MathContext mc = new MathContext(precision, rm);
BigDecimal computedRoot = scaledSquare.sqrt(mc);
failures += equalNumerically(expected, computedRoot, "simple squares");
int computedScale = computedRoot.scale();
if (precision >= expectedScale + 1 &&
computedScale != expectedScale) {
System.err.printf("%s\tprecision=%d\trm=%s%n",
computedRoot.toString(), precision, rm);
failures++;
System.err.printf("\t%s does not have expected scale of %d%n.",
computedRoot, expectedScale);
for (int scale = 0; scale <= 4; scale++) {
BigDecimal scaledPow = pow.setScale(scale);
int expectedScale = Math.ceilDiv(scale, n);
for (int precision = 0; precision <= 5; precision++) {
for (RoundingMode rm : RoundingMode.values()) {
MathContext mc = new MathContext(precision, rm);
BigDecimal computedRoot = scaledPow.rootn(n, mc);
failures += equalNumerically(expected, computedRoot, "simple powers");
int computedScale = computedRoot.scale();
if (precision >= expectedScale + 1 &&
computedScale != expectedScale) {
System.err.printf("%s\tprecision=%d\trm=%s%n",
computedRoot.toString(), precision, rm);
failures++;
System.err.printf("\t%s does not have expected scale of %d%n.",
computedRoot, expectedScale);
}
}
}
}
@ -372,7 +431,7 @@ public class SquareRootTests {
}
/**
* Test around 3.9999 that the sqrt doesn't improperly round-up to
* Test around 2^n-0.000...1 that the root doesn't improperly round-up to
* a numerical value of 2.
*/
private static int almostFourRoundingDown() {
@ -380,20 +439,30 @@ public class SquareRootTests {
BigDecimal nearFour = new BigDecimal("3.999999999999999999999999999999");
// Sqrt is 1.9999...
for (int i = 1; i < 64; i++) {
MathContext mc = new MathContext(i, RoundingMode.FLOOR);
BigDecimal result = nearFour.sqrt(mc);
BigDecimal expected = BigSquareRoot.sqrt(nearFour, mc);
failures += equalNumerically(expected, result, "near four rounding down");
failures += (result.compareTo(TWO) < 0) ? 0 : 1 ;
failures += (result.compareTo(TWO) < 0) ? 0 : 1;
}
BigDecimal ulp = nearFour.ulp();
for (int n = 1; n < 10; n++) {
BigDecimal near2ToN = BigDecimal.valueOf(Math.powExact(2, n)).subtract(ulp);
// Root is 1.9999...
for (int i = 1; i < 64; i++) {
MathContext mc = new MathContext(i, RoundingMode.FLOOR);
BigDecimal result = near2ToN.rootn(n, mc);
failures += (result.compareTo(TWO) < 0) ? 0 : 1;
}
}
return failures;
}
/**
* Test around 4.000...1 that the sqrt doesn't improperly
* Test around 2^n+0.000...1 that the root doesn't improperly
* round-down to a numerical value of 2.
*/
private static int almostFourRoundingUp() {
@ -401,13 +470,23 @@ public class SquareRootTests {
BigDecimal nearFour = new BigDecimal("4.000000000000000000000000000001");
// Sqrt is 2.0000....<non-zero digits>
for (int i = 1; i < 64; i++) {
MathContext mc = new MathContext(i, RoundingMode.CEILING);
BigDecimal result = nearFour.sqrt(mc);
BigDecimal expected = BigSquareRoot.sqrt(nearFour, mc);
failures += equalNumerically(expected, result, "near four rounding up");
failures += (result.compareTo(TWO) > 0) ? 0 : 1 ;
failures += (result.compareTo(TWO) > 0) ? 0 : 1;
}
BigDecimal ulp = nearFour.ulp();
for (int n = 1; n < 10; n++) {
BigDecimal near2ToN = BigDecimal.valueOf(Math.powExact(2, n)).add(ulp);
// Root is 2.0000....<non-zero digits>
for (int i = 1; i < 64; i++) {
MathContext mc = new MathContext(i, RoundingMode.CEILING);
BigDecimal result = near2ToN.rootn(n, mc);
failures += (result.compareTo(TWO) > 0) ? 0 : 1;
}
}
return failures;
@ -416,11 +495,11 @@ public class SquareRootTests {
private static int nearTen() {
int failures = 0;
BigDecimal near10 = new BigDecimal("9.99999999999999999999");
BigDecimal near10 = new BigDecimal("9.99999999999999999999");
BigDecimal near10sq = near10.multiply(near10);
BigDecimal near10sq = near10.multiply(near10);
BigDecimal near10sq_ulp = near10sq.add(near10sq.ulp());
BigDecimal near10sq_ulp = near10sq.add(near10sq.ulp());
for (int i = 10; i < 23; i++) {
MathContext mc = new MathContext(i, RoundingMode.HALF_EVEN);
@ -496,6 +575,27 @@ public class SquareRootTests {
}
}
for (int n = 1; n < 10; n++) {
for (BigDecimal halfWayCase : halfWayCases) {
// Round result to next-to-last place
int precision = halfWayCase.precision() - 1;
BigDecimal pow = halfWayCase.pow(n);
for (RoundingMode rm : List.of(RoundingMode.HALF_EVEN,
RoundingMode.HALF_UP,
RoundingMode.HALF_DOWN)) {
MathContext mc = new MathContext(precision, rm);
System.out.println("\nRounding mode " + rm);
System.out.println("\t" + halfWayCase.round(mc) + "\t" + halfWayCase);
failures += equalNumerically(pow.rootn(n, mc),
halfWayCase.round(mc),
"Rounding halway " + rm);
}
}
}
return failures;
}