< prev index next >
test/jdk/java/math/BigDecimal/SquareRootTests.java
Print this page
*** 1,7 ****
/*
! * 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.
--- 1,7 ----
/*
! * Copyright (c) 2016, 2019, 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.
*** 21,47 ****
* 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().");
}
--- 21,57 ----
* questions.
*/
/*
* @test
! * @bug 4851777 8233452
* @summary Tests of BigDecimal.sqrt().
*/
import java.math.*;
import java.util.*;
+ import static java.math.BigDecimal.TEN;
+ import static java.math.BigDecimal.ZERO;
+ import static java.math.BigDecimal.valueOf;
+
public class SquareRootTests {
+ private static BigDecimal TWO = new BigDecimal(2);
public static void main(String... args) {
int failures = 0;
failures += negativeTests();
failures += zeroTests();
failures += evenPowersOfTenTests();
failures += squareRootTwoTests();
failures += lowPrecisionPerfectSquares();
+ failures += almostFourRoundingDown();
+ failures += almostFourRoundingUp();
+ failures += nearTen();
+ failures += nearOne();
+ failures += halfWay();
if (failures > 0 ) {
throw new RuntimeException("Incurred " + failures + " failures" +
" testing BigDecimal.sqrt().");
}
*** 95,105 ****
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
--- 105,114 ----
*** 109,152 ****
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;
}
--- 118,160 ----
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;
// 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 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.
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);
! failures += equalNumerically(expected, computed, "sqrt(2)");
! failures += equalNumerically(computed, altComputed, "computed & altComputed");
}
}
return failures;
}
*** 193,202 ****
--- 201,339 ----
}
return failures;
}
+ /**
+ * Test around 3.9999 that the result doesn't not improperly
+ * round-up to a numerical value of 2.
+ */
+ private static int almostFourRoundingDown() {
+ int failures = 0;
+ BigDecimal nearFour = new BigDecimal("3.999999999999999999999999999999");
+
+ // Sqrt root 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 ;
+ }
+
+ return failures;
+ }
+
+ /**
+ * Test around 4.000...1 that the result doesn't not improperly
+ * round-down to a numerical value of 2.
+ */
+ private static int almostFourRoundingUp() {
+ int failures = 0;
+ BigDecimal nearFour = new BigDecimal("4.000000000000000000000000000001");
+
+ // Sqrt root 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 down");
+ failures += (result.compareTo(TWO) > 0) ? 0 : 1 ;
+ }
+
+ return failures;
+ }
+
+ private static int nearTen() {
+ int failures = 0;
+
+ BigDecimal near10 = new BigDecimal("9.99999999999999999999");
+
+ BigDecimal near10sq = near10.multiply(near10);
+
+ BigDecimal near10sq_ulp = near10sq.add(near10sq.ulp());
+
+ for (int i = 10; i < 23; i++) {
+ MathContext mc = new MathContext(i, RoundingMode.HALF_EVEN);
+
+ failures += equalNumerically(BigSquareRoot.sqrt(near10sq_ulp, mc),
+ near10sq_ulp.sqrt(mc),
+ "near 10 rounding down");
+ }
+
+ return failures;
+ }
+
+
+ /*
+ * Probe for rounding failures near a power of ten, 1 = 10^0,
+ * where an ulp has a different size above and below the value.
+ */
+ private static int nearOne() {
+ int failures = 0;
+
+ BigDecimal near1 = new BigDecimal(".999999999999999999999");
+ BigDecimal near1sq = near1.multiply(near1);
+ BigDecimal near1sq_ulp = near1sq.add(near1sq.ulp());
+
+ for (int i = 10; i < 23; i++) {
+ for (RoundingMode rm : List.of(RoundingMode.HALF_EVEN,
+ RoundingMode.UP,
+ RoundingMode.DOWN )) {
+ MathContext mc = new MathContext(i, rm);
+ failures += equalNumerically(BigSquareRoot.sqrt(near1sq_ulp, mc),
+ near1sq_ulp.sqrt(mc),
+ "near 1 half even");
+ }
+ }
+
+ return failures;
+ }
+
+
+ private static int halfWay() {
+ int failures = 0;
+
+ /*
+ * Use enough digits that the exact result cannot be computed
+ * from the sqrt of a double.
+ */
+ BigDecimal[] halfWayCases = {
+ // Odd next digit, truncate on HALF_EVEN
+ new BigDecimal("123456789123456789.5"),
+
+ // Even next digit, round up on HALF_EVEN
+ new BigDecimal("123456789123456788.5"),
+ };
+
+ for (BigDecimal halfWayCase : halfWayCases) {
+ // Round result to next-to-last place
+ int precision = halfWayCase.precision() - 1;
+ BigDecimal square = halfWayCase.multiply(halfWayCase);
+
+ 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);
+ /*System.out.println("\t" + square.sqrt(mc));*/
+ System.out.println("\t" + BigSquareRoot.sqrt(square, mc));
+
+ failures += equalNumerically(/*square.sqrt(mc),*/
+ BigSquareRoot.sqrt(square, mc),
+ halfWayCase.round(mc),
+ "Rounding halway " + rm);
+ }
+ }
+
+ 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 +
*** 222,227 ****
--- 359,572 ----
"\n\tExpected " + expected);
}
return failed;
}
+ /**
+ * Alternative implementation of BigDecimal square root which uses
+ * higher-precision for a simpler set of termination conditions
+ * for the Newton iteration.
+ */
+ private static class BigSquareRoot {
+ /**
+ * 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);
+
+ private static boolean isPowerOfTen(BigDecimal bd) {
+ return BigInteger.ONE.equals(bd.unscaledValue());
+ }
+
+ public static BigDecimal sqrt(BigDecimal bd, MathContext mc) {
+ int signum = bd.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 = bd.scale()/2;
+ BigDecimal zeroWithFinalPreferredScale =
+ BigDecimal.valueOf(0L, preferredScale);
+
+ // First phase of numerical normalization, strip trailing
+ // zeros and check for even powers of 10.
+ BigDecimal stripped = bd.stripTrailingZeros();
+ int strippedScale = stripped.scale();
+
+ // Numerically sqrt(10^2N) = 10^N
+ if (isPowerOfTen(stripped) &&
+ strippedScale % 2 == 0) {
+ BigDecimal result = BigDecimal.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();
+ // Use "2p + 2" property to guarantee enough
+ // intermediate precision so that a double-rounding
+ // error does not occur when rounded to the final
+ // destination precision.
+ int loopPrecision = Math.max(Math.max(2 * targetPrecision + 2,
+ workingPrecision),
+ 34); // Force at least
+ // two Netwon
+ // iterations on the
+ // Math.sqrt result.
+ do {
+ int tmpPrecision = Math.max(Math.max(guessPrecision, targetPrecision + 2),
+ workingPrecision);
+ MathContext mcTmp = new MathContext(loopPrecision, RoundingMode.HALF_EVEN);
+ // approx = 0.5 * (approx + fraction / approx)
+ approx = ONE_HALF.multiply(approx.add(working.divide(approx, mcTmp), mcTmp));
+ guessPrecision *= 2;
+ } while (guessPrecision < loopPrecision);
+
+ 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 (bd.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));
+ }
+ return result;
+ } else {
+ switch (signum) {
+ case -1:
+ throw new ArithmeticException("Attempted square root " +
+ "of negative BigDecimal");
+ case 0:
+ return valueOf(0L, bd.scale()/2);
+
+ default:
+ throw new AssertionError("Bad value from signum");
+ }
+ }
+ }
+ }
}
< prev index next >