| /* | |
| * Copyright (c) 2016, 2024, 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 8233452 8341402 | |
| * @summary Tests of BigDecimal.sqrt(). | |
| */ | |
| import java.math.BigDecimal; | |
| import java.math.BigInteger; | |
| import java.math.MathContext; | |
| import java.math.RoundingMode; | |
| import java.util.List; | |
| import static java.math.BigDecimal.ONE; | |
| import static java.math.BigDecimal.TWO; | |
| import static java.math.BigDecimal.TEN; | |
| import static java.math.BigDecimal.ZERO; | |
| import static java.math.BigDecimal.valueOf; | |
| public class SquareRootTests { | |
| /** | |
| * The value 0.1, with a scale of 1. | |
| */ | |
| private static final BigDecimal ONE_TENTH = valueOf(1L, 1); | |
| public static void main(String... args) { | |
| int failures = 0; | |
| failures += negativeTests(); | |
| failures += zeroTests(); | |
| failures += oneDigitTests(); | |
| failures += twoDigitTests(); | |
| failures += evenPowersOfTenTests(); | |
| failures += squareRootTwoTests(); | |
| failures += lowPrecisionPerfectSquares(); | |
| failures += almostFourRoundingDown(); | |
| failures += almostFourRoundingUp(); | |
| failures += nearTen(); | |
| failures += nearOne(); | |
| failures += halfWay(); | |
| failures += exactResultTests(); | |
| failures += scaleOverflowTest(); | |
| 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; | |
| } | |
| /** | |
| * 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); | |
| for (int i = 1; i < 20; i++) { | |
| for (RoundingMode rm : modes) { | |
| for (BigDecimal bd : oneToNine) { | |
| MathContext mc = new MathContext(i, rm); | |
| failures += compareSqrtImplementations(bd, mc); | |
| bd = bd.multiply(ONE_TENTH); | |
| failures += compareSqrtImplementations(bd, mc); | |
| } | |
| } | |
| } | |
| return failures; | |
| } | |
| /** | |
| * 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); | |
| for (int i = 10; i < 100; i++) { | |
| BigDecimal bd0 = BigDecimal.valueOf(i); | |
| BigDecimal bd1 = bd0.multiply(ONE_TENTH); | |
| BigDecimal bd2 = bd1.multiply(ONE_TENTH); | |
| for (BigDecimal bd : List.of(bd0, bd1, bd2)) { | |
| for (int precision = 1; i < 20; i++) { | |
| for (RoundingMode rm : modes) { | |
| MathContext mc = new MathContext(precision, rm); | |
| failures += compareSqrtImplementations(bd, mc); | |
| } | |
| } | |
| } | |
| } | |
| return failures; | |
| } | |
| private static int compareSqrtImplementations(BigDecimal bd, MathContext mc) { | |
| return equalNumerically(BigSquareRoot.sqrt(bd, mc), | |
| bd.sqrt(mc), "sqrt(" + bd + ") under " + mc); | |
| } | |
| private static int exactResultTests() { | |
| int failures = 0; | |
| 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, | |
| }; | |
| for (BigDecimal input : errCases) { | |
| BigDecimal result; | |
| // mc.roundingMode == RoundingMode.UNNECESSARY | |
| try { | |
| result = input.sqrt(unnecessary); | |
| System.err.println("Unexpected sqrt with UNNECESSARY RoundingMode: (" + input + ").sqrt() = " + result); | |
| failures += 1; | |
| } catch (ArithmeticException e) { | |
| // Expected | |
| } | |
| // mc.roundingMode != RoundingMode.UNNECESSARY && mc.precision == 0 | |
| try { | |
| result = input.sqrt(arbitrary); | |
| System.err.println("Unexpected sqrt with mc.precision == 0: (" + input + ").sqrt() = " + result); | |
| failures += 1; | |
| } catch (ArithmeticException e) { | |
| // Expected | |
| } | |
| } | |
| // (strippedScale & 1) == 0 | |
| // !stripped.isPowerOfTen() && sqrtRem[1].signum == 0 && (mc.precision != 0 && result.precision() > mc.precision) | |
| try { | |
| BigDecimal input = BigDecimal.valueOf(121); | |
| BigDecimal result = input.sqrt(unnecessary); | |
| System.err.println("Unexpected sqrt with result.precision() > mc.precision: (" | |
| + input + ").sqrt() = " + result); | |
| failures += 1; | |
| } catch (ArithmeticException e) { | |
| // Expected | |
| } | |
| BigDecimal four = BigDecimal.valueOf(4); | |
| 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.precision != 0 && result.precision() <= mc.precision) | |
| { four, unnecessary, BigDecimal.TWO }, | |
| }; | |
| for (Object[] testCase : cases) { | |
| BigDecimal expected = (BigDecimal) testCase[2]; | |
| BigDecimal result = ((BigDecimal) testCase[0]).sqrt((MathContext) testCase[1]); | |
| failures += compare(expected, result, true, "Exact results"); | |
| } | |
| return failures; | |
| } | |
| private static int scaleOverflowTest() { | |
| int failures = 0; | |
| try { | |
| BigDecimal.valueOf(1, -1).sqrt(new MathContext((1 << 30) + 1, RoundingMode.UP)); | |
| System.err.println("ArithmeticException expected: possible overflow undetected " | |
| + "or the range of supported values for the algorithm has extended."); | |
| failures += 1; | |
| } catch (ArithmeticException e) { | |
| // Expected | |
| } | |
| 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; | |
| // 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; | |
| } | |
| 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; | |
| } | |
| /** | |
| * Test around 3.9999 that the sqrt doesn't improperly round-up to | |
| * a numerical value of 2. | |
| */ | |
| private static int almostFourRoundingDown() { | |
| int failures = 0; | |
| 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 ; | |
| } | |
| return failures; | |
| } | |
| /** | |
| * Test around 4.000...1 that the sqrt doesn't improperly | |
| * round-down to a numerical value of 2. | |
| */ | |
| private static int almostFourRoundingUp() { | |
| int failures = 0; | |
| 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 ; | |
| } | |
| 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 half even"); | |
| } | |
| 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), | |
| mc.toString()); | |
| } | |
| } | |
| 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" + 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 + | |
| "(" + 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; | |
| } | |
| /** | |
| * 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.5, with a scale of 1. | |
| */ | |
| private static final BigDecimal ONE_HALF = valueOf(5L, 1); | |
| public static boolean isPowerOfTen(BigDecimal bd) { | |
| return BigInteger.ONE.equals(bd.unscaledValue()); | |
| } | |
| public static BigDecimal square(BigDecimal bd) { | |
| return bd.multiply(bd); | |
| } | |
| 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(2 * Math.max(targetPrecision, workingPrecision) + 2, | |
| 34); // Force at least two Netwon | |
| // iterations on the Math.sqrt | |
| // result. | |
| do { | |
| 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(square(result)).compareTo(ZERO) != 0) { | |
| throw new ArithmeticException("Computed square root not exact."); | |
| } | |
| } else { | |
| result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mc); | |
| } | |
| assert squareRootResultAssertions(bd, result, 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"); | |
| } | |
| } | |
| } | |
| /** | |
| * 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 static boolean squareRootResultAssertions(BigDecimal input, BigDecimal result, MathContext mc) { | |
| if (result.signum() == 0) { | |
| return squareRootZeroResultAssertions(input, 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 (isPowerOfTen(result)) { | |
| ulp = ulp.divide(TEN); | |
| } | |
| BigDecimal neighborDown = result.subtract(ulp); | |
| // Both the starting value and result should be nonzero and positive. | |
| if (result.signum() != 1 || | |
| input.signum() != 1) { | |
| return false; | |
| } | |
| switch (rm) { | |
| case DOWN: | |
| case FLOOR: | |
| assert | |
| square(result).compareTo(input) <= 0 && | |
| square(neighborUp).compareTo(input) > 0: | |
| "Square of result out for bounds rounding " + rm; | |
| return true; | |
| case UP: | |
| case CEILING: | |
| assert | |
| square(result).compareTo(input) >= 0 : | |
| "Square of result too small rounding " + rm; | |
| assert | |
| square(neighborDown).compareTo(input) < 0 : | |
| "Square of down neighbor too large rounding " + rm + "\n" + | |
| "\t input: " + input + "\t neighborDown: " + neighborDown +"\t sqrt: " + result + | |
| "\t" + mc; | |
| return true; | |
| case HALF_DOWN: | |
| case HALF_EVEN: | |
| case HALF_UP: | |
| BigDecimal err = square(result).subtract(input).abs(); | |
| BigDecimal errUp = square(neighborUp).subtract(input); | |
| BigDecimal errDown = input.subtract(square(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); | |
| assert | |
| errUp.signum() == 1 && | |
| errDown.signum() == 1 : | |
| "Errors of neighbors squared don't have correct signs"; | |
| // At least one of these must be true, but not both | |
| // assert | |
| // err_comp_errUp <= 0 : "Upper neighbor is closer than result: " + rm + | |
| // "\t" + input + "\t result" + result; | |
| // assert | |
| // err_comp_errDown <= 0 : "Lower neighbor is closer than result: " + rm + | |
| // "\t" + input + "\t result " + result + "\t lower neighbor: " + neighborDown; | |
| 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 static boolean squareRootZeroResultAssertions(BigDecimal input, | |
| BigDecimal result, | |
| MathContext mc) { | |
| return input.compareTo(ZERO) == 0; | |
| } | |
| } | |
| } | |