stan/jdk-24
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

4851777: Add BigDecimal sqrt method

Browse Source 
Reviewed-by: bpb
This commit is contained in:
Joe Darcy 2016-05-20 15:34:37 -07:00
parent 8c58aff49d
commit 4045a8be07
2 changed files with 527 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

300
jdk/src/java.base/share/classes/java/math/BigDecimal.java
View File

@ -128,6 +128,7 @@ import java.util.Arrays;
* <tr><td>Subtract</td><td>max(minuend.scale(), subtrahend.scale())</td>
* <tr><td>Multiply</td><td>multiplier.scale() + multiplicand.scale()</td>
* <tr><td>Divide</td><td>dividend.scale() - divisor.scale()</td>
* <tr><td>Square root</td><td>radicand.scale()/2</td>
* </table>
*
* These scales are the ones used by the methods which return exact
@ -346,6 +347,16 @@ public class BigDecimal extends Number implements Comparable<BigDecimal> {
public static final BigDecimal TEN =
ZERO_THROUGH_TEN[10];
/**
* The value 0.1, with a scale of 1.
*/
private static final BigDecimal ONE_TENTH = valueOf(1L, 1);
/**
* The value 0.5, with a scale of 1.
*/
private static final BigDecimal ONE_HALF = valueOf(5L, 1);
// Constructors
/**
@ -1995,6 +2006,295 @@ public class BigDecimal extends Number implements Comparable<BigDecimal> {
return result;
}
/**
* Returns an approximation to the square root of {@code this}
* with rounding according to the context settings.
*
* <p>The preferred scale of the returned result is equal to
* {@code this.scale()/2}. 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 square 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.sqrt(mc).equals(ZERO)} is true with
* any {@code MathContext} as an argument.
* </ul>
*
* @param mc the context to use.
* @return the square root of {@code this}.
* @throws ArithmeticException if {@code this} is less than 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.
* @since 9
*/
public BigDecimal sqrt(MathContext mc) {
int signum = signum();
if (signum == 1) {
/*
* The following code draws on the algorithm presented in
* "Properly Rounded Variable Precision Square Root," Hull and
* Abrham, ACM Transactions on Mathematical Software, Vol 11,
* No. 3, September 1985, Pages 229-237.
*
* The BigDecimal computational model differs from the one
* presented in the paper in several ways: first BigDecimal
* numbers aren't necessarily normalized, second many more
* rounding modes are supported, including UNNECESSARY, and
* exact results can be requested.
*
* The main steps of the algorithm below are as follows,
* first argument reduce the value to the numerical range
* [1, 10) using the following relations:
*
* 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
*
* Then use Newton's iteration 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.
int preferredScale = this.scale()/2;
BigDecimal zeroWithFinalPreferredScale = valueOf(0L, preferredScale);
// First phase of numerical normalization, strip trailing
// zeros and check for even powers of 10.
BigDecimal stripped = this.stripTrailingZeros();
int strippedScale = stripped.scale();
// Numerically sqrt(10^2N) = 10^N
if (stripped.isPowerOfTen() &&
strippedScale % 2 == 0) {
BigDecimal result = valueOf(1L, strippedScale/2);
if (result.scale() != preferredScale) {
// Adjust to requested precision and preferred
// scale as appropriate.
result = result.add(zeroWithFinalPreferredScale, mc);
}
return result;
}
// After stripTrailingZeros, the representation is normalized as
//
// unscaledValue * 10^(-scale)
//
// where unscaledValue is an integer with the mimimum
// precision for the cohort of the numerical value. To
// allow binary floating-point hardware to be used to get
// approximately a 15 digit approximation to the square
// root, it is helpful to instead normalize this so that
// the significand portion is to right of the decimal
// point by roughly (scale() - precision() +1).
// Now the precision / scale adjustment
int scaleAdjust = 0;
int scale = stripped.scale() - stripped.precision() + 1;
if (scale % 2 == 0) {
scaleAdjust = scale;
} else {
scaleAdjust = scale - 1;
}
BigDecimal working = stripped.scaleByPowerOfTen(scaleAdjust);
assert // Verify 0.1 <= working < 10
ONE_TENTH.compareTo(working) <= 0 && working.compareTo(TEN) < 0;
// Use good ole' Math.sqrt to get the initial guess for
// the Newton iteration, good to at least 15 decimal
// digits. This approach does incur the cost of a
//
// BigDecimal -> double -> BigDecimal
//
// conversion cycle, but it avoids the need for several
// Newton iterations in BigDecimal arithmetic to get the
// working answer to 15 digits of precision. If many fewer
// than 15 digits were needed, it might be faster to do
// the loop entirely in BigDecimal arithmetic.
//
// (A double value might have as much many as 17 decimal
// digits of precision; it depends on the relative density
// of binary and decimal numbers at different regions of
// the number line.)
//
// (It would be possible to check for certain special
// cases to avoid doing any Newton iterations. For
// example, if the BigDecimal -> double conversion was
// known to be exact and the rounding mode had a
// low-enough precision, the post-Newton rounding logic
// could be applied directly.)
BigDecimal guess = new BigDecimal(Math.sqrt(working.doubleValue()));
int guessPrecision = 15;
int originalPrecision = mc.getPrecision();
int targetPrecision;
// If an exact value is requested, it must only need about
// half of the input digits to represent since multiplying
// an N digit number by itself yield a 2N-1 digit or 2N
// digit result.
if (originalPrecision == 0) {
targetPrecision = stripped.precision()/2 + 1;
} else {
targetPrecision = originalPrecision;
}
// When setting the precision to use inside the Newton
// iteration loop, take care to avoid the case where the
// precision of the input exceeds the requested precision
// and rounding the input value too soon.
BigDecimal approx = guess;
int workingPrecision = working.precision();
do {
int tmpPrecision = Math.max(Math.max(guessPrecision, targetPrecision + 2),
workingPrecision);
MathContext mcTmp = new MathContext(tmpPrecision, RoundingMode.HALF_EVEN);
// approx = 0.5 * (approx + fraction / approx)
approx = ONE_HALF.multiply(approx.add(working.divide(approx, mcTmp), mcTmp));
guessPrecision *= 2;
} while (guessPrecision < targetPrecision + 2);
BigDecimal result;
RoundingMode targetRm = mc.getRoundingMode();
if (targetRm == RoundingMode.UNNECESSARY || originalPrecision == 0) {
RoundingMode tmpRm =
(targetRm == RoundingMode.UNNECESSARY) ? RoundingMode.DOWN : targetRm;
MathContext mcTmp = new MathContext(targetPrecision, tmpRm);
result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mcTmp);
// If result*result != this numerically, the square
// root isn't exact
if (this.subtract(result.multiply(result)).compareTo(ZERO) != 0) {
throw new ArithmeticException("Computed square root not exact.");
}
} else {
result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mc);
}
if (result.scale() != preferredScale) {
// The preferred scale of an add is
// max(addend.scale(), augend.scale()). Therefore, if
// the scale of the result is first minimized using
// stripTrailingZeros(), adding a zero of the
// preferred scale rounding the correct precision will
// perform the proper scale vs precision tradeoffs.
result = result.stripTrailingZeros().
add(zeroWithFinalPreferredScale,
new MathContext(originalPrecision, RoundingMode.UNNECESSARY));
}
assert squareRootResultAssertions(result, mc);
return result;
} else {
switch (signum) {
case -1:
throw new ArithmeticException("Attempted square root " +
"of negative BigDecimal");
case 0:
return valueOf(0L, scale()/2);
default:
throw new AssertionError("Bad value from signum");
}
}
}
private boolean isPowerOfTen() {
return BigInteger.ONE.equals(this.unscaledValue());
}
/**
* For nonzero values, check numerical correctness properties of
* the computed result for the chosen rounding mode.
*
* For the directed roundings, for DOWN and FLOOR, result^2 must
* be {@code <=} the input and (result+ulp)^2 must be {@code >} the
* input. Conversely, for UP and CEIL, result^2 must be {@code >=} the
* input and (result-ulp)^2 must be {@code <} the input.
*/
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 (this.isPowerOfTen()) {
ulp = ulp.divide(TEN);
}
BigDecimal neighborDown = result.subtract(ulp);
// Both the starting value and result should be nonzero and positive.
if (result.signum() != 1 ||
this.signum() != 1) {
return false;
}
switch (rm) {
case DOWN:
case FLOOR:
return
result.multiply(result).compareTo(this) <= 0 &&
neighborUp.multiply(neighborUp).compareTo(this) > 0;
case UP:
case CEILING:
return
result.multiply(result).compareTo(this) >= 0 &&
neighborDown.multiply(neighborDown).compareTo(this) < 0;
case HALF_DOWN:
case HALF_EVEN:
case HALF_UP:
BigDecimal err = result.multiply(result).subtract(this).abs();
BigDecimal errUp = neighborUp.multiply(neighborUp).subtract(this);
BigDecimal errDown = this.subtract(neighborDown.multiply(neighborDown));
// 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);
return
errUp.signum() == 1 &&
errDown.signum() == 1 &&
err_comp_errUp <= 0 &&
err_comp_errDown <= 0 &&
((err_comp_errUp == 0 ) ? err_comp_errDown < 0 : true) &&
((err_comp_errDown == 0 ) ? err_comp_errUp < 0 : true);
// && could check for digit conditions for ties too
default: // Definition of UNNECESSARY already verified.
return true;
}
}
}
private boolean squareRootZeroResultAssertions(BigDecimal result, MathContext mc) {
return this.compareTo(ZERO) == 0;
}
/**
* Returns a {@code BigDecimal} whose value is
* <code>(this<sup>n</sup>)</code>, The power is computed exactly, to

227
jdk/test/java/math/BigDecimal/SquareRootTests.java Normal file
View File

@ -0,0 +1,227 @@
/*
* Copyright (c) 2016, 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
* under the terms of the GNU General Public License version 2 only, as
* published by the Free Software Foundation.
*
* This code is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
* version 2 for more details (a copy is included in the LICENSE file that
* accompanied this code).
*
* You should have received a copy of the GNU General Public License version
* 2 along with this work; if not, write to the Free Software Foundation,
* Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
*
* Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
* or visit www.oracle.com if you need additional information or have any
* questions.
*/
/*
* @test
* @bug 4851777
* @summary Tests of BigDecimal.sqrt().
*/
import java.math.*;
import java.util.*;
public class SquareRootTests {
public static void main(String... args) {
int failures = 0;
failures += negativeTests();
failures += zeroTests();
failures += evenPowersOfTenTests();
failures += squareRootTwoTests();
failures += lowPrecisionPerfectSquares();
if (failures > 0 ) {
throw new RuntimeException("Incurred " + failures + " failures" +
" testing BigDecimal.sqrt().");
}
}
private static int negativeTests() {
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
}
}
}
return failures;
}
private static int zeroTests() {
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");
failures += compare(BigDecimal.valueOf(0L, i).sqrt(MathContext.DECIMAL64),
expected, true, "zeros");
}
return failures;
}
/**
* sqrt(10^2N) is 10^N
* Both numerical value and representation should be verified
*/
private static int evenPowersOfTenTests() {
int failures = 0;
MathContext oneDigitExactly = new MathContext(1, RoundingMode.UNNECESSARY);
for (int scale = -100; scale <= 100; scale++) {
BigDecimal testValue = BigDecimal.valueOf(1, 2*scale);
BigDecimal expectedNumericalResult = BigDecimal.valueOf(1, scale);
BigDecimal result;
failures += equalNumerically(expectedNumericalResult,
result = testValue.sqrt(MathContext.DECIMAL64),
"Even powers of 10, DECIMAL64");
// Can round to one digit of precision exactly
failures += equalNumerically(expectedNumericalResult,
result = testValue.sqrt(oneDigitExactly),
"even powers of 10, 1 digit");
if (result.precision() > 1) {
failures += 1;
System.err.println("Excess precision for " + result);
}
// If rounding to more than one digit, do precision / scale checking...
}
return failures;
}
private static int squareRootTwoTests() {
int failures = 0;
BigDecimal TWO = new BigDecimal(2);
// Square root of 2 truncated to 65 digits
BigDecimal highPrecisionRoot2 =
new BigDecimal("1.41421356237309504880168872420969807856967187537694807317667973799");
RoundingMode[] modes = {
RoundingMode.UP, RoundingMode.DOWN,
RoundingMode.CEILING, RoundingMode.FLOOR,
RoundingMode.HALF_UP, RoundingMode.HALF_DOWN, RoundingMode.HALF_EVEN
};
// For each iteresting rounding mode, for precisions 1 to, say
// 63 numerically compare TWO.sqrt(mc) to
// highPrecisionRoot2.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);
equalNumerically(expected, computed, "sqrt(2)");
}
}
return failures;
}
private static int lowPrecisionPerfectSquares() {
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 (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 = 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);
}
}
}
}
}
return failures;
}
private static int compare(BigDecimal a, BigDecimal b, boolean expected, String prefix) {
boolean result = a.equals(b);
int failed = (result==expected) ? 0 : 1;
if (failed == 1) {
System.err.println("Testing " + prefix +
"(" + a + ").compareTo(" + b + ") => " + result +
"\n\tExpected " + expected);
}
return failed;
}
private static int equalNumerically(BigDecimal a, BigDecimal b,
String prefix) {
return compareNumerically(a, b, 0, prefix);
}
private static int compareNumerically(BigDecimal a, BigDecimal b,
int expected, String prefix) {
int result = a.compareTo(b);
int failed = (result==expected) ? 0 : 1;
if (failed == 1) {
System.err.println("Testing " + prefix +
"(" + a + ").compareTo(" + b + ") => " + result +
"\n\tExpected " + expected);
}
return failed;
}
}