Path: blob/master/test/jdk/java/math/BigDecimal/SquareRootTests.java
41149 views
/*1* Copyright (c) 2016, 2020, Oracle and/or its affiliates. All rights reserved.2* DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.3*4* This code is free software; you can redistribute it and/or modify it5* under the terms of the GNU General Public License version 2 only, as6* published by the Free Software Foundation.7*8* This code is distributed in the hope that it will be useful, but WITHOUT9* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or10* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License11* version 2 for more details (a copy is included in the LICENSE file that12* accompanied this code).13*14* You should have received a copy of the GNU General Public License version15* 2 along with this work; if not, write to the Free Software Foundation,16* Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.17*18* Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA19* or visit www.oracle.com if you need additional information or have any20* questions.21*/2223/*24* @test25* @bug 4851777 823345226* @summary Tests of BigDecimal.sqrt().27*/2829import java.math.*;30import java.util.*;3132import static java.math.BigDecimal.ONE;33import static java.math.BigDecimal.TEN;34import static java.math.BigDecimal.ZERO;35import static java.math.BigDecimal.valueOf;3637public class SquareRootTests {38private static BigDecimal TWO = new BigDecimal(2);3940/**41* The value 0.1, with a scale of 1.42*/43private static final BigDecimal ONE_TENTH = valueOf(1L, 1);4445public static void main(String... args) {46int failures = 0;4748failures += negativeTests();49failures += zeroTests();50failures += oneDigitTests();51failures += twoDigitTests();52failures += evenPowersOfTenTests();53failures += squareRootTwoTests();54failures += lowPrecisionPerfectSquares();55failures += almostFourRoundingDown();56failures += almostFourRoundingUp();57failures += nearTen();58failures += nearOne();59failures += halfWay();6061if (failures > 0 ) {62throw new RuntimeException("Incurred " + failures + " failures" +63" testing BigDecimal.sqrt().");64}65}6667private static int negativeTests() {68int failures = 0;6970for (long i = -10; i < 0; i++) {71for (int j = -5; j < 5; j++) {72try {73BigDecimal input = BigDecimal.valueOf(i, j);74BigDecimal result = input.sqrt(MathContext.DECIMAL64);75System.err.println("Unexpected sqrt of negative: (" +76input + ").sqrt() = " + result );77failures += 1;78} catch (ArithmeticException e) {79; // Expected80}81}82}8384return failures;85}8687private static int zeroTests() {88int failures = 0;8990for (int i = -100; i < 100; i++) {91BigDecimal expected = BigDecimal.valueOf(0L, i/2);92// These results are independent of rounding mode93failures += compare(BigDecimal.valueOf(0L, i).sqrt(MathContext.UNLIMITED),94expected, true, "zeros");9596failures += compare(BigDecimal.valueOf(0L, i).sqrt(MathContext.DECIMAL64),97expected, true, "zeros");98}99100return failures;101}102103/**104* Probe inputs with one digit of precision, 1 ... 9 and those105* values scaled by 10^-1, 0.1, ... 0.9.106*/107private static int oneDigitTests() {108int failures = 0;109110List<BigDecimal> oneToNine =111List.of(ONE, TWO, valueOf(3),112valueOf(4), valueOf(5), valueOf(6),113valueOf(7), valueOf(8), valueOf(9));114115List<RoundingMode> modes =116List.of(RoundingMode.UP, RoundingMode.DOWN,117RoundingMode.CEILING, RoundingMode.FLOOR,118RoundingMode.HALF_UP, RoundingMode.HALF_DOWN, RoundingMode.HALF_EVEN);119120for (int i = 1; i < 20; i++) {121for (RoundingMode rm : modes) {122for (BigDecimal bd : oneToNine) {123MathContext mc = new MathContext(i, rm);124125failures += compareSqrtImplementations(bd, mc);126bd = bd.multiply(ONE_TENTH);127failures += compareSqrtImplementations(bd, mc);128}129}130}131132return failures;133}134135/**136* Probe inputs with two digits of precision, (10 ... 99) and137* those values scaled by 10^-1 (1, ... 9.9) and scaled by 10^-2138* (0.1 ... 0.99).139*/140private static int twoDigitTests() {141int failures = 0;142143List<RoundingMode> modes =144List.of(RoundingMode.UP, RoundingMode.DOWN,145RoundingMode.CEILING, RoundingMode.FLOOR,146RoundingMode.HALF_UP, RoundingMode.HALF_DOWN, RoundingMode.HALF_EVEN);147148for (int i = 10; i < 100; i++) {149BigDecimal bd0 = BigDecimal.valueOf(i);150BigDecimal bd1 = bd0.multiply(ONE_TENTH);151BigDecimal bd2 = bd1.multiply(ONE_TENTH);152153for (BigDecimal bd : List.of(bd0, bd1, bd2)) {154for (int precision = 1; i < 20; i++) {155for (RoundingMode rm : modes) {156MathContext mc = new MathContext(precision, rm);157failures += compareSqrtImplementations(bd, mc);158}159}160}161}162163return failures;164}165166private static int compareSqrtImplementations(BigDecimal bd, MathContext mc) {167return equalNumerically(BigSquareRoot.sqrt(bd, mc),168bd.sqrt(mc), "sqrt(" + bd + ") under " + mc);169}170171/**172* sqrt(10^2N) is 10^N173* Both numerical value and representation should be verified174*/175private static int evenPowersOfTenTests() {176int failures = 0;177MathContext oneDigitExactly = new MathContext(1, RoundingMode.UNNECESSARY);178179for (int scale = -100; scale <= 100; scale++) {180BigDecimal testValue = BigDecimal.valueOf(1, 2*scale);181BigDecimal expectedNumericalResult = BigDecimal.valueOf(1, scale);182183BigDecimal result;184185failures += equalNumerically(expectedNumericalResult,186result = testValue.sqrt(MathContext.DECIMAL64),187"Even powers of 10, DECIMAL64");188189// Can round to one digit of precision exactly190failures += equalNumerically(expectedNumericalResult,191result = testValue.sqrt(oneDigitExactly),192"even powers of 10, 1 digit");193194if (result.precision() > 1) {195failures += 1;196System.err.println("Excess precision for " + result);197}198199// If rounding to more than one digit, do precision / scale checking...200}201202return failures;203}204205private static int squareRootTwoTests() {206int failures = 0;207208// Square root of 2 truncated to 65 digits209BigDecimal highPrecisionRoot2 =210new BigDecimal("1.41421356237309504880168872420969807856967187537694807317667973799");211212RoundingMode[] modes = {213RoundingMode.UP, RoundingMode.DOWN,214RoundingMode.CEILING, RoundingMode.FLOOR,215RoundingMode.HALF_UP, RoundingMode.HALF_DOWN, RoundingMode.HALF_EVEN216};217218219// For each interesting rounding mode, for precisions 1 to, say,220// 63 numerically compare TWO.sqrt(mc) to221// highPrecisionRoot2.round(mc) and the alternative internal high-precision222// implementation of square root.223for (RoundingMode mode : modes) {224for (int precision = 1; precision < 63; precision++) {225MathContext mc = new MathContext(precision, mode);226BigDecimal expected = highPrecisionRoot2.round(mc);227BigDecimal computed = TWO.sqrt(mc);228BigDecimal altComputed = BigSquareRoot.sqrt(TWO, mc);229230failures += equalNumerically(expected, computed, "sqrt(2)");231failures += equalNumerically(computed, altComputed, "computed & altComputed");232}233}234235return failures;236}237238private static int lowPrecisionPerfectSquares() {239int failures = 0;240241// For 5^2 through 9^2, if the input is rounded to one digit242// first before the root is computed, the wrong answer will243// result. Verify results and scale for different rounding244// modes and precisions.245long[][] squaresWithOneDigitRoot = {{ 4, 2},246{ 9, 3},247{25, 5},248{36, 6},249{49, 7},250{64, 8},251{81, 9}};252253for (long[] squareAndRoot : squaresWithOneDigitRoot) {254BigDecimal square = new BigDecimal(squareAndRoot[0]);255BigDecimal expected = new BigDecimal(squareAndRoot[1]);256257for (int scale = 0; scale <= 4; scale++) {258BigDecimal scaledSquare = square.setScale(scale, RoundingMode.UNNECESSARY);259int expectedScale = scale/2;260for (int precision = 0; precision <= 5; precision++) {261for (RoundingMode rm : RoundingMode.values()) {262MathContext mc = new MathContext(precision, rm);263BigDecimal computedRoot = scaledSquare.sqrt(mc);264failures += equalNumerically(expected, computedRoot, "simple squares");265int computedScale = computedRoot.scale();266if (precision >= expectedScale + 1 &&267computedScale != expectedScale) {268System.err.printf("%s\tprecision=%d\trm=%s%n",269computedRoot.toString(), precision, rm);270failures++;271System.err.printf("\t%s does not have expected scale of %d%n.",272computedRoot, expectedScale);273}274}275}276}277}278279return failures;280}281282/**283* Test around 3.9999 that the sqrt doesn't improperly round-up to284* a numerical value of 2.285*/286private static int almostFourRoundingDown() {287int failures = 0;288BigDecimal nearFour = new BigDecimal("3.999999999999999999999999999999");289290// Sqrt is 1.9999...291292for (int i = 1; i < 64; i++) {293MathContext mc = new MathContext(i, RoundingMode.FLOOR);294BigDecimal result = nearFour.sqrt(mc);295BigDecimal expected = BigSquareRoot.sqrt(nearFour, mc);296failures += equalNumerically(expected, result, "near four rounding down");297failures += (result.compareTo(TWO) < 0) ? 0 : 1 ;298}299300return failures;301}302303/**304* Test around 4.000...1 that the sqrt doesn't improperly305* round-down to a numerical value of 2.306*/307private static int almostFourRoundingUp() {308int failures = 0;309BigDecimal nearFour = new BigDecimal("4.000000000000000000000000000001");310311// Sqrt is 2.0000....<non-zero digits>312313for (int i = 1; i < 64; i++) {314MathContext mc = new MathContext(i, RoundingMode.CEILING);315BigDecimal result = nearFour.sqrt(mc);316BigDecimal expected = BigSquareRoot.sqrt(nearFour, mc);317failures += equalNumerically(expected, result, "near four rounding up");318failures += (result.compareTo(TWO) > 0) ? 0 : 1 ;319}320321return failures;322}323324private static int nearTen() {325int failures = 0;326327BigDecimal near10 = new BigDecimal("9.99999999999999999999");328329BigDecimal near10sq = near10.multiply(near10);330331BigDecimal near10sq_ulp = near10sq.add(near10sq.ulp());332333for (int i = 10; i < 23; i++) {334MathContext mc = new MathContext(i, RoundingMode.HALF_EVEN);335336failures += equalNumerically(BigSquareRoot.sqrt(near10sq_ulp, mc),337near10sq_ulp.sqrt(mc),338"near 10 rounding half even");339}340341return failures;342}343344345/*346* Probe for rounding failures near a power of ten, 1 = 10^0,347* where an ulp has a different size above and below the value.348*/349private static int nearOne() {350int failures = 0;351352BigDecimal near1 = new BigDecimal(".999999999999999999999");353BigDecimal near1sq = near1.multiply(near1);354BigDecimal near1sq_ulp = near1sq.add(near1sq.ulp());355356for (int i = 10; i < 23; i++) {357for (RoundingMode rm : List.of(RoundingMode.HALF_EVEN,358RoundingMode.UP,359RoundingMode.DOWN )) {360MathContext mc = new MathContext(i, rm);361failures += equalNumerically(BigSquareRoot.sqrt(near1sq_ulp, mc),362near1sq_ulp.sqrt(mc),363mc.toString());364}365}366367return failures;368}369370371private static int halfWay() {372int failures = 0;373374/*375* Use enough digits that the exact result cannot be computed376* from the sqrt of a double.377*/378BigDecimal[] halfWayCases = {379// Odd next digit, truncate on HALF_EVEN380new BigDecimal("123456789123456789.5"),381382// Even next digit, round up on HALF_EVEN383new BigDecimal("123456789123456788.5"),384};385386for (BigDecimal halfWayCase : halfWayCases) {387// Round result to next-to-last place388int precision = halfWayCase.precision() - 1;389BigDecimal square = halfWayCase.multiply(halfWayCase);390391for (RoundingMode rm : List.of(RoundingMode.HALF_EVEN,392RoundingMode.HALF_UP,393RoundingMode.HALF_DOWN)) {394MathContext mc = new MathContext(precision, rm);395396System.out.println("\nRounding mode " + rm);397System.out.println("\t" + halfWayCase.round(mc) + "\t" + halfWayCase);398System.out.println("\t" + BigSquareRoot.sqrt(square, mc));399400failures += equalNumerically(/*square.sqrt(mc),*/401BigSquareRoot.sqrt(square, mc),402halfWayCase.round(mc),403"Rounding halway " + rm);404}405}406407return failures;408}409410private static int compare(BigDecimal a, BigDecimal b, boolean expected, String prefix) {411boolean result = a.equals(b);412int failed = (result==expected) ? 0 : 1;413if (failed == 1) {414System.err.println("Testing " + prefix +415"(" + a + ").compareTo(" + b + ") => " + result +416"\n\tExpected " + expected);417}418return failed;419}420421private static int equalNumerically(BigDecimal a, BigDecimal b,422String prefix) {423return compareNumerically(a, b, 0, prefix);424}425426427private static int compareNumerically(BigDecimal a, BigDecimal b,428int expected, String prefix) {429int result = a.compareTo(b);430int failed = (result==expected) ? 0 : 1;431if (failed == 1) {432System.err.println("Testing " + prefix +433"(" + a + ").compareTo(" + b + ") => " + result +434"\n\tExpected " + expected);435}436return failed;437}438439/**440* Alternative implementation of BigDecimal square root which uses441* higher-precision for a simpler set of termination conditions442* for the Newton iteration.443*/444private static class BigSquareRoot {445446/**447* The value 0.5, with a scale of 1.448*/449private static final BigDecimal ONE_HALF = valueOf(5L, 1);450451public static boolean isPowerOfTen(BigDecimal bd) {452return BigInteger.ONE.equals(bd.unscaledValue());453}454455public static BigDecimal square(BigDecimal bd) {456return bd.multiply(bd);457}458459public static BigDecimal sqrt(BigDecimal bd, MathContext mc) {460int signum = bd.signum();461if (signum == 1) {462/*463* The following code draws on the algorithm presented in464* "Properly Rounded Variable Precision Square Root," Hull and465* Abrham, ACM Transactions on Mathematical Software, Vol 11,466* No. 3, September 1985, Pages 229-237.467*468* The BigDecimal computational model differs from the one469* presented in the paper in several ways: first BigDecimal470* numbers aren't necessarily normalized, second many more471* rounding modes are supported, including UNNECESSARY, and472* exact results can be requested.473*474* The main steps of the algorithm below are as follows,475* first argument reduce the value to the numerical range476* [1, 10) using the following relations:477*478* x = y * 10 ^ exp479* sqrt(x) = sqrt(y) * 10^(exp / 2) if exp is even480* sqrt(x) = sqrt(y/10) * 10 ^((exp+1)/2) is exp is odd481*482* Then use Newton's iteration on the reduced value to compute483* the numerical digits of the desired result.484*485* Finally, scale back to the desired exponent range and486* perform any adjustment to get the preferred scale in the487* representation.488*/489490// The code below favors relative simplicity over checking491// for special cases that could run faster.492493int preferredScale = bd.scale()/2;494BigDecimal zeroWithFinalPreferredScale =495BigDecimal.valueOf(0L, preferredScale);496497// First phase of numerical normalization, strip trailing498// zeros and check for even powers of 10.499BigDecimal stripped = bd.stripTrailingZeros();500int strippedScale = stripped.scale();501502// Numerically sqrt(10^2N) = 10^N503if (isPowerOfTen(stripped) &&504strippedScale % 2 == 0) {505BigDecimal result = BigDecimal.valueOf(1L, strippedScale/2);506if (result.scale() != preferredScale) {507// Adjust to requested precision and preferred508// scale as appropriate.509result = result.add(zeroWithFinalPreferredScale, mc);510}511return result;512}513514// After stripTrailingZeros, the representation is normalized as515//516// unscaledValue * 10^(-scale)517//518// where unscaledValue is an integer with the mimimum519// precision for the cohort of the numerical value. To520// allow binary floating-point hardware to be used to get521// approximately a 15 digit approximation to the square522// root, it is helpful to instead normalize this so that523// the significand portion is to right of the decimal524// point by roughly (scale() - precision() + 1).525526// Now the precision / scale adjustment527int scaleAdjust = 0;528int scale = stripped.scale() - stripped.precision() + 1;529if (scale % 2 == 0) {530scaleAdjust = scale;531} else {532scaleAdjust = scale - 1;533}534535BigDecimal working = stripped.scaleByPowerOfTen(scaleAdjust);536537assert // Verify 0.1 <= working < 10538ONE_TENTH.compareTo(working) <= 0 && working.compareTo(TEN) < 0;539540// Use good ole' Math.sqrt to get the initial guess for541// the Newton iteration, good to at least 15 decimal542// digits. This approach does incur the cost of a543//544// BigDecimal -> double -> BigDecimal545//546// conversion cycle, but it avoids the need for several547// Newton iterations in BigDecimal arithmetic to get the548// working answer to 15 digits of precision. If many fewer549// than 15 digits were needed, it might be faster to do550// the loop entirely in BigDecimal arithmetic.551//552// (A double value might have as much many as 17 decimal553// digits of precision; it depends on the relative density554// of binary and decimal numbers at different regions of555// the number line.)556//557// (It would be possible to check for certain special558// cases to avoid doing any Newton iterations. For559// example, if the BigDecimal -> double conversion was560// known to be exact and the rounding mode had a561// low-enough precision, the post-Newton rounding logic562// could be applied directly.)563564BigDecimal guess = new BigDecimal(Math.sqrt(working.doubleValue()));565int guessPrecision = 15;566int originalPrecision = mc.getPrecision();567int targetPrecision;568569// If an exact value is requested, it must only need570// about half of the input digits to represent since571// multiplying an N digit number by itself yield a (2N572// - 1) digit or 2N digit result.573if (originalPrecision == 0) {574targetPrecision = stripped.precision()/2 + 1;575} else {576targetPrecision = originalPrecision;577}578579// When setting the precision to use inside the Newton580// iteration loop, take care to avoid the case where the581// precision of the input exceeds the requested precision582// and rounding the input value too soon.583BigDecimal approx = guess;584int workingPrecision = working.precision();585// Use "2p + 2" property to guarantee enough586// intermediate precision so that a double-rounding587// error does not occur when rounded to the final588// destination precision.589int loopPrecision =590Math.max(2 * Math.max(targetPrecision, workingPrecision) + 2,59134); // Force at least two Netwon592// iterations on the Math.sqrt593// result.594do {595MathContext mcTmp = new MathContext(loopPrecision, RoundingMode.HALF_EVEN);596// approx = 0.5 * (approx + fraction / approx)597approx = ONE_HALF.multiply(approx.add(working.divide(approx, mcTmp), mcTmp));598guessPrecision *= 2;599} while (guessPrecision < loopPrecision);600601BigDecimal result;602RoundingMode targetRm = mc.getRoundingMode();603if (targetRm == RoundingMode.UNNECESSARY || originalPrecision == 0) {604RoundingMode tmpRm =605(targetRm == RoundingMode.UNNECESSARY) ? RoundingMode.DOWN : targetRm;606MathContext mcTmp = new MathContext(targetPrecision, tmpRm);607result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mcTmp);608609// If result*result != this numerically, the square610// root isn't exact611if (bd.subtract(square(result)).compareTo(ZERO) != 0) {612throw new ArithmeticException("Computed square root not exact.");613}614} else {615result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mc);616}617618assert squareRootResultAssertions(bd, result, mc);619if (result.scale() != preferredScale) {620// The preferred scale of an add is621// max(addend.scale(), augend.scale()). Therefore, if622// the scale of the result is first minimized using623// stripTrailingZeros(), adding a zero of the624// preferred scale rounding the correct precision will625// perform the proper scale vs precision tradeoffs.626result = result.stripTrailingZeros().627add(zeroWithFinalPreferredScale,628new MathContext(originalPrecision, RoundingMode.UNNECESSARY));629}630return result;631} else {632switch (signum) {633case -1:634throw new ArithmeticException("Attempted square root " +635"of negative BigDecimal");636case 0:637return valueOf(0L, bd.scale()/2);638639default:640throw new AssertionError("Bad value from signum");641}642}643}644645/**646* For nonzero values, check numerical correctness properties of647* the computed result for the chosen rounding mode.648*649* For the directed roundings, for DOWN and FLOOR, result^2 must650* be {@code <=} the input and (result+ulp)^2 must be {@code >} the651* input. Conversely, for UP and CEIL, result^2 must be {@code >=} the652* input and (result-ulp)^2 must be {@code <} the input.653*/654private static boolean squareRootResultAssertions(BigDecimal input, BigDecimal result, MathContext mc) {655if (result.signum() == 0) {656return squareRootZeroResultAssertions(input, result, mc);657} else {658RoundingMode rm = mc.getRoundingMode();659BigDecimal ulp = result.ulp();660BigDecimal neighborUp = result.add(ulp);661// Make neighbor down accurate even for powers of ten662if (isPowerOfTen(result)) {663ulp = ulp.divide(TEN);664}665BigDecimal neighborDown = result.subtract(ulp);666667// Both the starting value and result should be nonzero and positive.668if (result.signum() != 1 ||669input.signum() != 1) {670return false;671}672673switch (rm) {674case DOWN:675case FLOOR:676assert677square(result).compareTo(input) <= 0 &&678square(neighborUp).compareTo(input) > 0:679"Square of result out for bounds rounding " + rm;680return true;681682case UP:683case CEILING:684assert685square(result).compareTo(input) >= 0 :686"Square of result too small rounding " + rm;687688assert689square(neighborDown).compareTo(input) < 0 :690"Square of down neighbor too large rounding " + rm + "\n" +691"\t input: " + input + "\t neighborDown: " + neighborDown +"\t sqrt: " + result +692"\t" + mc;693return true;694695696case HALF_DOWN:697case HALF_EVEN:698case HALF_UP:699BigDecimal err = square(result).subtract(input).abs();700BigDecimal errUp = square(neighborUp).subtract(input);701BigDecimal errDown = input.subtract(square(neighborDown));702// All error values should be positive so don't need to703// compare absolute values.704705int err_comp_errUp = err.compareTo(errUp);706int err_comp_errDown = err.compareTo(errDown);707708assert709errUp.signum() == 1 &&710errDown.signum() == 1 :711"Errors of neighbors squared don't have correct signs";712713// At least one of these must be true, but not both714// assert715// err_comp_errUp <= 0 : "Upper neighbor is closer than result: " + rm +716// "\t" + input + "\t result" + result;717// assert718// err_comp_errDown <= 0 : "Lower neighbor is closer than result: " + rm +719// "\t" + input + "\t result " + result + "\t lower neighbor: " + neighborDown;720721assert722((err_comp_errUp == 0 ) ? err_comp_errDown < 0 : true) &&723((err_comp_errDown == 0 ) ? err_comp_errUp < 0 : true) :724"Incorrect error relationships";725// && could check for digit conditions for ties too726return true;727728default: // Definition of UNNECESSARY already verified.729return true;730}731}732}733734private static boolean squareRootZeroResultAssertions(BigDecimal input,735BigDecimal result,736MathContext mc) {737return input.compareTo(ZERO) == 0;738}739}740}741742743744