test/jdk/java/math/BigDecimal/SquareRootTests.java - master - githubmirror/panama-vector

test/jdk/java/math/BigDecimal/SquareRootTests.java
933 строки39 KB
Fabio Romano
8366478: BigDecimal roots
08 май 2026, 02:59
08 май 2026, 02:598ae19e1
/*
 * Copyright (c) 2016, 2026, 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.*;
import java.util.stream.IntStream;

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 += zerothRootTests();
        failures += negativeWithEvenDegreeTests();
        failures += zeroTests();
        failures += oneDigitTests();
        failures += twoDigitTests();
        failures += evenPowersOfTenTests();
        failures += rootTwoTests();
        failures += lowPrecisionPerfectPowers();
        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 zerothRootTests() {
        int failures = 0;

        for (long i = -5; i < 5; i++) {
            for (int j = -5; j < 5; j++) {
                try {
                    BigDecimal input = BigDecimal.valueOf(i, j);
                    BigDecimal result = input.rootn(0, MathContext.DECIMAL64);
                    System.err.println("Unexpected 0th root: (" +
                                       input + ").rootn(0)  = " + result );
                    failures += 1;
                } catch (ArithmeticException e) {
                    ; // Expected
                }
            }
        }

        return failures;
    }

    private static int negativeWithEvenDegreeTests() {
        int failures = 0;

        for (long i = -10; i < 0; i++) {
            for (int j = -5; j < 5; j++) {
                BigDecimal input = BigDecimal.valueOf(i, j);
                for (int n = -4; n <= 4; n += 2) {
                    try {
                        BigDecimal result = input.rootn(n, MathContext.DECIMAL64);
                        System.err.println("Unexpected nth root of negative: (" +
                                           input + ").rootn(" + n + ")  = " + result );
                        failures += 1;
                    } catch (ArithmeticException e) {
                        ; // Expected
                    }
                }
            }
        }

        return failures;
    }

    private static int zeroTests() {
        int failures = 0;

        for (int i = -100; i < 100; i++) {
            BigDecimal input = BigDecimal.valueOf(0L, i);
            for (int n = -10; n <= 0; n++) {
                try {
                    BigDecimal result = input.rootn(n, MathContext.DECIMAL64);
                    System.err.println("Unexpected nth root of zero: (" +
                                       input + ").rootn(" + n + ")  = " + result );
                    failures += 1;
                } catch (ArithmeticException e) {
                    ; // Expected
                }
            }

            for (int n = 1; n < 10; n++) {
                BigDecimal expected = BigDecimal.valueOf(0L, Math.ceilDiv(i, n));
                // These results are independent of rounding mode
                failures += compare(input.rootn(n, MathContext.UNLIMITED),
                                    expected, true, "zeros");

                failures += compare(input.rootn(n, MathContext.DECIMAL64),
                                    expected, true, "zeros");
            }
        }

        return failures;
    }

    private static RoundingMode positiveRoundingMode(RoundingMode rm) {
        return rm == RoundingMode.FLOOR   ? RoundingMode.UP :
               rm == RoundingMode.CEILING ? RoundingMode.DOWN : rm;
    }

    private static List<RoundingMode> roundingModes() {
        List<RoundingMode> modes = new ArrayList<>(Arrays.asList(RoundingMode.values()));
        modes.remove(RoundingMode.UNNECESSARY);
        return modes;
    }

    /**
     * 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 = IntStream.rangeClosed(1, 9).mapToObj(BigDecimal::valueOf).toList();
        List<RoundingMode> modes = roundingModes();

        for (int i = 1; i < 20; i++) {
            for (RoundingMode rm : modes) {
                MathContext mc = new MathContext(i, rm);
                MathContext positiveMC = new MathContext(i, positiveRoundingMode(rm));
                for (BigDecimal bd  : oneToNine) {
                    failures += compareSqrtImplementations(bd, mc);
                    BigDecimal minus_bd = bd.negate();
                    for (int n = 1; n < 11; n += 2) {
                        failures += compare(minus_bd.rootn( n, mc), bd.rootn( n, positiveMC).negate(), true, "one digit");
                        failures += compare(minus_bd.rootn(-n, mc), bd.rootn(-n, positiveMC).negate(), true, "one digit");
                    }

                    bd = bd.scaleByPowerOfTen(-1);
                    failures += compareSqrtImplementations(bd, mc);
                    minus_bd = bd.negate();
                    for (int n = 1; n < 11; n += 2) {
                        failures += compare(minus_bd.rootn( n, mc), bd.rootn( n, positiveMC).negate(), true, "one digit");
                        failures += compare(minus_bd.rootn(-n, mc), bd.rootn(-n, positiveMC).negate(), true, "one digit");
                    }
                }
            }
        }

        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 = roundingModes();

        for (int i = 10; i < 100; i++) {
            BigDecimal bd0 = BigDecimal.valueOf(i);
            BigDecimal bd1 = bd0.scaleByPowerOfTen(-1);
            BigDecimal bd2 = bd1.scaleByPowerOfTen(-1);

            for (BigDecimal bd : List.of(bd0, bd1, bd2)) {
                for (int prec = 1; prec < 20; prec++) {
                    for (RoundingMode rm : modes) {
                        MathContext mc = new MathContext(prec, rm);
                        MathContext positiveMC = new MathContext(prec, positiveRoundingMode(rm));
                        failures += compareSqrtImplementations(bd, mc);

                        BigDecimal minus_bd = bd.negate();
                        for (int n = 1; n < 11; n += 2) {
                            failures += compare(minus_bd.rootn( n, mc), bd.rootn( n, positiveMC).negate(), true, "two digits");
                            failures += compare(minus_bd.rootn(-n, mc), bd.rootn(-n, positiveMC).negate(), true, "two digits");
                        }
                    }
                }
            }
        }

        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);

        Object[][] errCases = {
                // strippedScale % n != 0
                { BigDecimal.TEN, 2 },
                // strippedScale % n == 0 && sqrtRem[1].signum != 0
                { BigDecimal.TWO, 2 },
                // sqrtRem[1].signum == 0 && n < 0
                { BigDecimal.valueOf(9L), -2 },
        };

        for (Object[] errCase : errCases) {
            BigDecimal input = (BigDecimal) errCase[0];
            int n = (int) errCase[1];
            BigDecimal result;
            // mc.roundingMode == RoundingMode.UNNECESSARY
            try {
                result = input.rootn(n, unnecessary);
                System.err.println("Unexpected nth root with UNNECESSARY RoundingMode: ("
                + input + ").rootn(" + n + ") = " + result);
                failures += 1;
            } catch (ArithmeticException e) {
                // Expected
            }

            // mc.roundingMode != RoundingMode.UNNECESSARY && mc.precision == 0
            try {
                result = input.rootn(n, arbitrary);
                System.err.println("Unexpected nth root with mc.precision == 0: ("
                + input + ").rootn(" + n + ") = " + result);
                failures += 1;
            } catch (ArithmeticException e) {
                // Expected
            }
        }

        // strippedScale % n == 0

        // sqrtRem[1].signum == 0 && n > 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), oneHalf = BigDecimal.valueOf(5, 1);
        Object[][] cases = {
                // mc.roundingMode == RoundingMode.UNNECESSARY && sqrtRem[1].signum == 0 && n > 0
                // && mc.precision == 0
                { four, 2, new MathContext(0, RoundingMode.UNNECESSARY), BigDecimal.TWO },
                // mc.roundingMode != RoundingMode.UNNECESSARY && sqrtRem[1].signum == 0 && n > 0
                // && mc.precision == 0
                { four, 2, arbitrary, BigDecimal.TWO },
                // sqrtRem[1].signum == 0 && n > 0
                // && (mc.precision != 0 && result.precision() <= mc.precision)
                { four, 2, unnecessary, BigDecimal.TWO },
                // mc.roundingMode == RoundingMode.UNNECESSARY && sqrtRem[1].signum == 0 && n < 0
                { four, -2, unnecessary, oneHalf },
                // mc.roundingMode != RoundingMode.UNNECESSARY && sqrtRem[1].signum == 0 && n < 0
                // && mc.precision == 0
                { four, -2, arbitrary, oneHalf },
        };

        for (Object[] testCase : cases) {
            BigDecimal expected = (BigDecimal) testCase[3];
            BigDecimal result = ((BigDecimal) testCase[0]).rootn((int) testCase[1], (MathContext) testCase[2]);
            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 rootTwoTests() {
        int failures = 0;

        // Square root of 2 truncated to 65 digits
        BigDecimal highPrecisionSqrt2 =
            new BigDecimal("1.41421356237309504880168872420969807856967187537694807317667973799");
        // Cube root of 2 truncated to 65 digits
        BigDecimal highPrecisionCbrt2 =
            new BigDecimal("1.25992104989487316476721060727822835057025146470150798008197511215");

        List<RoundingMode> modes = roundingModes();


        // For each interesting rounding mode, for precisions 1 to, say,
        // 63 numerically compare TWO.sqrt(mc) to
        // highPrecisionSqrt2.round(mc) and the alternative internal high-precision
        // implementation of square root, and compare TWO.rootn(3, mc) to
        // highPrecisionCbrt2.round(mc).
        for (RoundingMode mode : modes) {
            for (int precision = 1; precision < 63; precision++) {
                MathContext mc = new MathContext(precision, mode);
                BigDecimal expectedSqrt = highPrecisionSqrt2.round(mc);
                BigDecimal computedSqrt = TWO.sqrt(mc);
                BigDecimal altComputedSqrt = BigSquareRoot.sqrt(TWO, mc);

                failures += equalNumerically(expectedSqrt, computedSqrt, "sqrt(2)");
                failures += equalNumerically(computedSqrt, altComputedSqrt, "computed & altComputed");

                BigDecimal expectedCbrt = highPrecisionCbrt2.round(mc);
                BigDecimal computedCbrt = TWO.rootn(3, mc);

                failures += equalNumerically(expectedCbrt, computedCbrt, "rootn(2, 3)");
            }
        }

        return failures;
    }

    private static int lowPrecisionPerfectPowers() {
        int failures = 0;

        for (int i = -9; i <= 9; i++) {
            for (int n = 1; n < 10; n++) {
                BigDecimal expected = BigDecimal.valueOf(n % 2 != 0 ? i : Math.abs(i));
                BigDecimal pow = BigDecimal.valueOf(Math.powExact(i, n));

                for (int scale = 0; scale <= 4; scale++) {
                    BigDecimal scaledPow = pow.setScale(scale);
                    int expectedScale = Math.ceilDiv(scale, n);
                    for (int precision = 0; precision <= 5; precision++) {
                        for (RoundingMode rm : RoundingMode.values()) {
                            MathContext mc = new MathContext(precision, rm);
                            BigDecimal computedRoot = scaledPow.rootn(n, mc);
                            failures += equalNumerically(expected, computedRoot, "simple powers");
                            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 2^n-0.000...1 that the root 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;
        }

        BigDecimal ulp = nearFour.ulp();
        for (int n = 1; n < 10; n++) {
            BigDecimal near2ToN = BigDecimal.valueOf(Math.powExact(2, n)).subtract(ulp);
            // Root is 1.9999...
            for (int i = 1; i < 64; i++) {
                MathContext mc = new MathContext(i, RoundingMode.FLOOR);
                BigDecimal result = near2ToN.rootn(n, mc);
                failures += (result.compareTo(TWO) < 0) ? 0 : 1;
            }
        }

        return failures;
    }

    /**
     * Test around 2^n+0.000...1 that the root 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;
        }

        BigDecimal ulp = nearFour.ulp();
        for (int n = 1; n < 10; n++) {
            BigDecimal near2ToN = BigDecimal.valueOf(Math.powExact(2, n)).add(ulp);
            // Root is 2.0000....<non-zero digits>
            for (int i = 1; i < 64; i++) {
                MathContext mc = new MathContext(i, RoundingMode.CEILING);
                BigDecimal result = near2ToN.rootn(n, mc);
                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);
            }
        }

        for (int n = 1; n < 10; n++) {
            for (BigDecimal halfWayCase : halfWayCases) {
                // Round result to next-to-last place
                int precision = halfWayCase.precision() - 1;
                BigDecimal pow = halfWayCase.pow(n);

                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);

                    failures += equalNumerically(pow.rootn(n, 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 = Math.ceilDiv(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;
        }
    }
}