--- old/test/jdk/java/math/BigDecimal/SquareRootTests.java 2019-11-27 06:27:40.616000000 -0800 +++ new/test/jdk/java/math/BigDecimal/SquareRootTests.java 2019-11-27 06:27:40.352000000 -0800 @@ -1,5 +1,5 @@ /* - * Copyright (c) 2016, Oracle and/or its affiliates. All rights reserved. + * 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 @@ -23,23 +23,40 @@ /* * @test - * @bug 4851777 + * @bug 4851777 8233452 * @summary Tests of BigDecimal.sqrt(). */ import java.math.*; import java.util.*; +import static java.math.BigDecimal.ONE; +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); + + /** + * 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 += evenPowersOfTenTests(); failures += squareRootTwoTests(); failures += lowPrecisionPerfectSquares(); + failures += almostFourRoundingDown(); + failures += almostFourRoundingUp(); + failures += nearTen(); + failures += nearOne(); + failures += halfWay(); if (failures > 0 ) { throw new RuntimeException("Incurred " + failures + " failures" + @@ -84,6 +101,40 @@ } /** + * Probe inputs with one digit of precision, 1 .... 9 and those + * values scaled by 10^-1. + */ + private static int oneDigitTests() { + int failures = 0; + + List oneToNine = + List.of(ONE, valueOf(2), valueOf(3), + valueOf(4), valueOf(5), valueOf(6), + valueOf(7), valueOf(8), valueOf(9)); + + List 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 += equalNumerically(BigSquareRoot.sqrt(bd, mc), + bd.sqrt(mc), "sqrt(" + bd + ") under " + mc); + bd = bd.multiply(ONE_TENTH); + failures += equalNumerically(BigSquareRoot.sqrt(bd, mc), + bd.sqrt(mc), "sqrt(" + bd + ") under " + mc); + } + } + } + + return failures; + } + + /** * sqrt(10^2N) is 10^N * Both numerical value and representation should be verified */ @@ -97,7 +148,6 @@ BigDecimal result; - failures += equalNumerically(expectedNumericalResult, result = testValue.sqrt(MathContext.DECIMAL64), "Even powers of 10, DECIMAL64"); @@ -111,9 +161,7 @@ System.err.println("Excess precision for " + result); } - // If rounding to more than one digit, do precision / scale checking... - } return failures; @@ -121,30 +169,31 @@ 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 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); - equalNumerically(expected, computed, "sqrt(2)"); + failures += equalNumerically(expected, computed, "sqrt(2)"); + failures += equalNumerically(computed, altComputed, "computed & altComputed"); } } @@ -181,8 +230,8 @@ int computedScale = computedRoot.scale(); if (precision >= expectedScale + 1 && computedScale != expectedScale) { - System.err.printf("%s\tprecision=%d\trm=%s%n", - computedRoot.toString(), precision, rm); + 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); @@ -195,6 +244,135 @@ 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.... + + 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; @@ -224,4 +402,305 @@ 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 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); + } + + 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 + result.multiply(result).compareTo(input) <= 0 && + neighborUp.multiply(neighborUp).compareTo(input) > 0: + "Square of result out for bounds rounding " + rm; + return true; + + case UP: + case CEILING: + assert + result.multiply(result).compareTo(input) >= 0 : + "Square of result too small rounding " + rm; + + assert + neighborDown.multiply(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 = result.multiply(result).subtract(input).abs(); + BigDecimal errUp = neighborUp.multiply(neighborUp).subtract(input); + BigDecimal errDown = input.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); + + assert + errUp.signum() == 1 && + errDown.signum() == 1 : + "Errors of neighbors squared don't have correct signs"; + + 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; + } + } } +