--- old/src/java.base/share/classes/java/math/BigDecimal.java 2016-05-11 10:37:28.126558223 -0700 +++ new/src/java.base/share/classes/java/math/BigDecimal.java 2016-05-11 10:37:27.730558229 -0700 @@ -346,6 +346,16 @@ 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 +2005,219 @@ return result; } + + /** + * Returns an approximation to the square root of {@code this} + * with rounding according to the context settings. + * + *

The preferred scale of the returned result is equal to + * {@code floor(this.scale()/2.0)}. 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. + * + *

Special case: + *

+ * + * @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 follow, 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. + + BigDecimal zeroWithFinalPreferredScale = valueOf(0L, this.scale()/2); + + // 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 (BigInteger.ONE.equals(stripped.unscaledValue()) && + strippedScale % 2 == 0) { + BigDecimal result = valueOf(1L, strippedScale/2); + // Adjust to requested precision and preferred scale as possible + return result.add(zeroWithFinalPreferredScale, mc); + } + + System.out.println("\tStarting: " + this); // DEBUG + System.out.println("\tStripped: " + stripped); // DEBUG + + + // 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 as so + // that the significand portion is to right of the decimal + // point, roughly: + // + // (unscaledValue * (10^-precision) * 10^(-scale)) * (10^precision) + // + // so that + // + // sqrt(unscaledValue * (10^-precision) * 10^(-scale) * (10^precision)) = + // + // sqrt(unscaledValue * (10^-precision) * 10^(-scale)) * 10^(precision/2) + // + // Therefore, this adjustment occurs for by 10^-precision is precision is even or + // (adjust as needed, +/-1 + // + // (A double value might have as much as 17 decimal digits + // of precision; it depends on the relative density of + // binary and decimal numbers at different points of the + // number line.) + + // 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; + stripped = stripped.divide(TEN); + } + + // At least document potential exponent overflow issues here??? + BigDecimal working = stripped.scaleByPowerOfTen(scaleAdjust); + System.out.println("\tWorking2: " + working + "\tscale: " + working.scale()); // DEBUG + + // Verify 1 <= working < 10 (verify range) + assert working.compareTo(ONE) >= 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. + BigDecimal guess = new BigDecimal(Math.sqrt(working.doubleValue())); + int guessPrecision = 15; + int targetPrecision = mc.getPrecision(); + RoundingMode rm = mc.getRoundingMode(); + + // When setting the precision to use inside the 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 { + System.out.println(approx); // DEBUG + 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 = approx.round(mc); + + // TODO: need code for Round.UNNECESSARY. Probably also + // need additinal checking for HALF_EVEN since we are only + // using two extra precision digits. + + // For now, skip the rounding up / down fix up + + // Read through "What every ..." for 2p + 2 discussion... + + // Handle 1/2 ulp, 1 ulp, and unnecessary rounding separately? + // UP, DOWN, + // CEILING (equiv to UP), FLOOR (equiv to DOWN), + // + // HALF_UP, HALF_DOWN, HALF_EVEN -- look up 2p + 2 derivation + // -- same for decimal as well as binary? + + // Could use direct definition to test adjacent values to + // target precision -- need to lookup if Newton's iteration + // converges from above or below to avoid checking three values... + + // UNNECESSARY + + // Unnecessary rounding -- in the worst case, only need about + // (1/2 + 1) digits to produce an exact output, sqrt(100) = + // 10, sqrt(10_000) = 100, etc.. How to handle cases that are + // powers of 10 with trailing zeros? + + // Normalize to a range [.1, 1] or [1, 10] -- presunmably this + // is open on one end to account for exact square roots... + // + + + // Convert to double and take double sqrt -- easy way to get + // (at least) 15 digits of a decimal approximation :-) + + // If did more work, could make sure came in from uniformly + // above / below the root? + + // Final answer: rescale fraction and then add a zero with the + // preferred scale with a rounding precision equal to the + // target precision. + + result = result.scaleByPowerOfTen(-scaleAdjust/2); + return result.add(zeroWithFinalPreferredScale, + new MathContext(mc.getPrecision(), RoundingMode.UNNECESSARY)); + } 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"); + } + } + } + /** * Returns a {@code BigDecimal} whose value is * (thisn), The power is computed exactly, to --- /dev/null 2016-04-30 12:24:07.320478178 -0700 +++ new/test/java/math/BigDecimal/SquareRootTests.java 2016-05-11 10:37:28.902558211 -0700 @@ -0,0 +1,192 @@ +/* + * 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(); + + // Add trickier rounding cases -- try to find examples of + + // When setting the precision to use inside the loop, take + // care to avoid the case where the precision of the input + // exceeds the requested precision and rounding the input + // value too soon. + + + 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 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; + } + +}