1 /*
2  * Copyright (c) 2016, 2020, Oracle and/or its affiliates. All rights reserved.
3  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4  *
5  * This code is free software; you can redistribute it and/or modify it
6  * under the terms of the GNU General Public License version 2 only, as
7  * published by the Free Software Foundation.
8  *
9  * This code is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
12  * version 2 for more details (a copy is included in the LICENSE file that
13  * accompanied this code).
14  *
15  * You should have received a copy of the GNU General Public License version
16  * 2 along with this work; if not, write to the Free Software Foundation,
17  * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
18  *
19  * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
20  * or visit www.oracle.com if you need additional information or have any
21  * questions.
22  */
23 
24 /*
25  * @test
26  * @bug 4851777 8233452
27  * @summary Tests of BigDecimal.sqrt().
28  */
29 
30 import java.math.*;
31 import java.util.*;
32 
33 import static java.math.BigDecimal.ONE;
34 import static java.math.BigDecimal.TEN;
35 import static java.math.BigDecimal.ZERO;
36 import static java.math.BigDecimal.valueOf;
37 
38 public class SquareRootTests {
39     private static BigDecimal TWO = new BigDecimal(2);
40 
41     /**
42      * The value 0.1, with a scale of 1.
43      */
44     private static final BigDecimal ONE_TENTH = valueOf(1L, 1);
45 
main(String... args)46     public static void main(String... args) {
47         int failures = 0;
48 
49         failures += negativeTests();
50         failures += zeroTests();
51         failures += oneDigitTests();
52         failures += twoDigitTests();
53         failures += evenPowersOfTenTests();
54         failures += squareRootTwoTests();
55         failures += lowPrecisionPerfectSquares();
56         failures += almostFourRoundingDown();
57         failures += almostFourRoundingUp();
58         failures += nearTen();
59         failures += nearOne();
60         failures += halfWay();
61 
62         if (failures > 0 ) {
63             throw new RuntimeException("Incurred " + failures + " failures" +
64                                        " testing BigDecimal.sqrt().");
65         }
66     }
67 
negativeTests()68     private static int negativeTests() {
69         int failures = 0;
70 
71         for (long i = -10; i < 0; i++) {
72             for (int j = -5; j < 5; j++) {
73                 try {
74                     BigDecimal input = BigDecimal.valueOf(i, j);
75                     BigDecimal result = input.sqrt(MathContext.DECIMAL64);
76                     System.err.println("Unexpected sqrt of negative: (" +
77                                        input + ").sqrt()  = " + result );
78                     failures += 1;
79                 } catch (ArithmeticException e) {
80                     ; // Expected
81                 }
82             }
83         }
84 
85         return failures;
86     }
87 
zeroTests()88     private static int zeroTests() {
89         int failures = 0;
90 
91         for (int i = -100; i < 100; i++) {
92             BigDecimal expected = BigDecimal.valueOf(0L, i/2);
93             // These results are independent of rounding mode
94             failures += compare(BigDecimal.valueOf(0L, i).sqrt(MathContext.UNLIMITED),
95                                 expected, true, "zeros");
96 
97             failures += compare(BigDecimal.valueOf(0L, i).sqrt(MathContext.DECIMAL64),
98                                 expected, true, "zeros");
99         }
100 
101         return failures;
102     }
103 
104     /**
105      * Probe inputs with one digit of precision, 1 ... 9 and those
106      * values scaled by 10^-1, 0.1, ... 0.9.
107      */
oneDigitTests()108     private static int oneDigitTests() {
109         int failures = 0;
110 
111         List<BigDecimal> oneToNine =
112             List.of(ONE,        TWO,        valueOf(3),
113                     valueOf(4), valueOf(5), valueOf(6),
114                     valueOf(7), valueOf(8), valueOf(9));
115 
116         List<RoundingMode> modes =
117             List.of(RoundingMode.UP,      RoundingMode.DOWN,
118                     RoundingMode.CEILING, RoundingMode.FLOOR,
119                     RoundingMode.HALF_UP, RoundingMode.HALF_DOWN, RoundingMode.HALF_EVEN);
120 
121         for (int i = 1; i < 20; i++) {
122             for (RoundingMode rm : modes) {
123                 for (BigDecimal bd  : oneToNine) {
124                     MathContext mc = new MathContext(i, rm);
125 
126                     failures += compareSqrtImplementations(bd, mc);
127                     bd = bd.multiply(ONE_TENTH);
128                     failures += compareSqrtImplementations(bd, mc);
129                 }
130             }
131         }
132 
133         return failures;
134     }
135 
136     /**
137      * Probe inputs with two digits of precision, (10 ... 99) and
138      * those values scaled by 10^-1 (1, ... 9.9) and scaled by 10^-2
139      * (0.1 ... 0.99).
140      */
twoDigitTests()141     private static int twoDigitTests() {
142         int failures = 0;
143 
144         List<RoundingMode> modes =
145             List.of(RoundingMode.UP,      RoundingMode.DOWN,
146                     RoundingMode.CEILING, RoundingMode.FLOOR,
147                     RoundingMode.HALF_UP, RoundingMode.HALF_DOWN, RoundingMode.HALF_EVEN);
148 
149         for (int i = 10; i < 100; i++) {
150             BigDecimal bd0 = BigDecimal.valueOf(i);
151             BigDecimal bd1 = bd0.multiply(ONE_TENTH);
152             BigDecimal bd2 = bd1.multiply(ONE_TENTH);
153 
154             for (BigDecimal bd : List.of(bd0, bd1, bd2)) {
155                 for (int precision = 1; i < 20; i++) {
156                     for (RoundingMode rm : modes) {
157                         MathContext mc = new MathContext(precision, rm);
158                         failures += compareSqrtImplementations(bd, mc);
159                     }
160                 }
161             }
162         }
163 
164         return failures;
165     }
166 
compareSqrtImplementations(BigDecimal bd, MathContext mc)167     private static int compareSqrtImplementations(BigDecimal bd, MathContext mc) {
168         return equalNumerically(BigSquareRoot.sqrt(bd, mc),
169                                 bd.sqrt(mc), "sqrt(" + bd + ") under " + mc);
170     }
171 
172     /**
173      * sqrt(10^2N) is 10^N
174      * Both numerical value and representation should be verified
175      */
evenPowersOfTenTests()176     private static int evenPowersOfTenTests() {
177         int failures = 0;
178         MathContext oneDigitExactly = new MathContext(1, RoundingMode.UNNECESSARY);
179 
180         for (int scale = -100; scale <= 100; scale++) {
181             BigDecimal testValue               = BigDecimal.valueOf(1, 2*scale);
182             BigDecimal expectedNumericalResult = BigDecimal.valueOf(1,   scale);
183 
184             BigDecimal result;
185 
186             failures += equalNumerically(expectedNumericalResult,
187                                          result = testValue.sqrt(MathContext.DECIMAL64),
188                                          "Even powers of 10, DECIMAL64");
189 
190             // Can round to one digit of precision exactly
191             failures += equalNumerically(expectedNumericalResult,
192                                          result = testValue.sqrt(oneDigitExactly),
193                                          "even powers of 10, 1 digit");
194 
195             if (result.precision() > 1) {
196                 failures += 1;
197                 System.err.println("Excess precision for " + result);
198             }
199 
200             // If rounding to more than one digit, do precision / scale checking...
201         }
202 
203         return failures;
204     }
205 
squareRootTwoTests()206     private static int squareRootTwoTests() {
207         int failures = 0;
208 
209         // Square root of 2 truncated to 65 digits
210         BigDecimal highPrecisionRoot2 =
211             new BigDecimal("1.41421356237309504880168872420969807856967187537694807317667973799");
212 
213         RoundingMode[] modes = {
214             RoundingMode.UP,       RoundingMode.DOWN,
215             RoundingMode.CEILING, RoundingMode.FLOOR,
216             RoundingMode.HALF_UP, RoundingMode.HALF_DOWN, RoundingMode.HALF_EVEN
217         };
218 
219 
220         // For each interesting rounding mode, for precisions 1 to, say,
221         // 63 numerically compare TWO.sqrt(mc) to
222         // highPrecisionRoot2.round(mc) and the alternative internal high-precision
223         // implementation of square root.
224         for (RoundingMode mode : modes) {
225             for (int precision = 1; precision < 63; precision++) {
226                 MathContext mc = new MathContext(precision, mode);
227                 BigDecimal expected = highPrecisionRoot2.round(mc);
228                 BigDecimal computed = TWO.sqrt(mc);
229                 BigDecimal altComputed = BigSquareRoot.sqrt(TWO, mc);
230 
231                 failures += equalNumerically(expected, computed, "sqrt(2)");
232                 failures += equalNumerically(computed, altComputed, "computed & altComputed");
233             }
234         }
235 
236         return failures;
237     }
238 
lowPrecisionPerfectSquares()239     private static int lowPrecisionPerfectSquares() {
240         int failures = 0;
241 
242         // For 5^2 through 9^2, if the input is rounded to one digit
243         // first before the root is computed, the wrong answer will
244         // result. Verify results and scale for different rounding
245         // modes and precisions.
246         long[][] squaresWithOneDigitRoot = {{ 4, 2},
247                                             { 9, 3},
248                                             {25, 5},
249                                             {36, 6},
250                                             {49, 7},
251                                             {64, 8},
252                                             {81, 9}};
253 
254         for (long[] squareAndRoot : squaresWithOneDigitRoot) {
255             BigDecimal square     = new BigDecimal(squareAndRoot[0]);
256             BigDecimal expected   = new BigDecimal(squareAndRoot[1]);
257 
258             for (int scale = 0; scale <= 4; scale++) {
259                 BigDecimal scaledSquare = square.setScale(scale, RoundingMode.UNNECESSARY);
260                 int expectedScale = scale/2;
261                 for (int precision = 0; precision <= 5; precision++) {
262                     for (RoundingMode rm : RoundingMode.values()) {
263                         MathContext mc = new MathContext(precision, rm);
264                         BigDecimal computedRoot = scaledSquare.sqrt(mc);
265                         failures += equalNumerically(expected, computedRoot, "simple squares");
266                         int computedScale = computedRoot.scale();
267                         if (precision >=  expectedScale + 1 &&
268                             computedScale != expectedScale) {
269                             System.err.printf("%s\tprecision=%d\trm=%s%n",
270                                               computedRoot.toString(), precision, rm);
271                             failures++;
272                             System.err.printf("\t%s does not have expected scale of %d%n.",
273                                               computedRoot, expectedScale);
274                         }
275                     }
276                 }
277             }
278         }
279 
280         return failures;
281     }
282 
283     /**
284      * Test around 3.9999 that the sqrt doesn't improperly round-up to
285      * a numerical value of 2.
286      */
almostFourRoundingDown()287     private static int almostFourRoundingDown() {
288         int failures = 0;
289         BigDecimal nearFour = new BigDecimal("3.999999999999999999999999999999");
290 
291         // Sqrt is 1.9999...
292 
293         for (int i = 1; i < 64; i++) {
294             MathContext mc = new MathContext(i, RoundingMode.FLOOR);
295             BigDecimal result = nearFour.sqrt(mc);
296             BigDecimal expected = BigSquareRoot.sqrt(nearFour, mc);
297             failures += equalNumerically(expected, result, "near four rounding down");
298             failures += (result.compareTo(TWO) < 0) ? 0  : 1 ;
299         }
300 
301         return failures;
302     }
303 
304     /**
305      * Test around 4.000...1 that the sqrt doesn't improperly
306      * round-down to a numerical value of 2.
307      */
almostFourRoundingUp()308     private static int almostFourRoundingUp() {
309         int failures = 0;
310         BigDecimal nearFour = new BigDecimal("4.000000000000000000000000000001");
311 
312         // Sqrt is 2.0000....<non-zero digits>
313 
314         for (int i = 1; i < 64; i++) {
315             MathContext mc = new MathContext(i, RoundingMode.CEILING);
316             BigDecimal result = nearFour.sqrt(mc);
317             BigDecimal expected = BigSquareRoot.sqrt(nearFour, mc);
318             failures += equalNumerically(expected, result, "near four rounding up");
319             failures += (result.compareTo(TWO) > 0) ? 0  : 1 ;
320         }
321 
322         return failures;
323     }
324 
nearTen()325     private static int nearTen() {
326         int failures = 0;
327 
328          BigDecimal near10 = new BigDecimal("9.99999999999999999999");
329 
330          BigDecimal near10sq = near10.multiply(near10);
331 
332          BigDecimal near10sq_ulp = near10sq.add(near10sq.ulp());
333 
334         for (int i = 10; i < 23; i++) {
335             MathContext mc = new MathContext(i, RoundingMode.HALF_EVEN);
336 
337             failures += equalNumerically(BigSquareRoot.sqrt(near10sq_ulp, mc),
338                                          near10sq_ulp.sqrt(mc),
339                                          "near 10 rounding half even");
340         }
341 
342         return failures;
343     }
344 
345 
346     /*
347      * Probe for rounding failures near a power of ten, 1 = 10^0,
348      * where an ulp has a different size above and below the value.
349      */
nearOne()350     private static int nearOne() {
351         int failures = 0;
352 
353          BigDecimal near1 = new BigDecimal(".999999999999999999999");
354          BigDecimal near1sq = near1.multiply(near1);
355          BigDecimal near1sq_ulp = near1sq.add(near1sq.ulp());
356 
357          for (int i = 10; i < 23; i++) {
358              for (RoundingMode rm : List.of(RoundingMode.HALF_EVEN,
359                                             RoundingMode.UP,
360                                             RoundingMode.DOWN )) {
361                  MathContext mc = new MathContext(i, rm);
362                  failures += equalNumerically(BigSquareRoot.sqrt(near1sq_ulp, mc),
363                                               near1sq_ulp.sqrt(mc),
364                                               mc.toString());
365              }
366          }
367 
368          return failures;
369     }
370 
371 
halfWay()372     private static int halfWay() {
373         int failures = 0;
374 
375         /*
376          * Use enough digits that the exact result cannot be computed
377          * from the sqrt of a double.
378          */
379         BigDecimal[] halfWayCases = {
380             // Odd next digit, truncate on HALF_EVEN
381             new BigDecimal("123456789123456789.5"),
382 
383              // Even next digit, round up on HALF_EVEN
384             new BigDecimal("123456789123456788.5"),
385         };
386 
387         for (BigDecimal halfWayCase : halfWayCases) {
388             // Round result to next-to-last place
389             int precision = halfWayCase.precision() - 1;
390             BigDecimal square = halfWayCase.multiply(halfWayCase);
391 
392             for (RoundingMode rm : List.of(RoundingMode.HALF_EVEN,
393                                            RoundingMode.HALF_UP,
394                                            RoundingMode.HALF_DOWN)) {
395                 MathContext mc = new MathContext(precision, rm);
396 
397                 System.out.println("\nRounding mode " + rm);
398                 System.out.println("\t" + halfWayCase.round(mc) + "\t" + halfWayCase);
399                 System.out.println("\t" + BigSquareRoot.sqrt(square, mc));
400 
401                 failures += equalNumerically(/*square.sqrt(mc),*/
402                                              BigSquareRoot.sqrt(square, mc),
403                                              halfWayCase.round(mc),
404                                              "Rounding halway " + rm);
405             }
406         }
407 
408         return failures;
409     }
410 
compare(BigDecimal a, BigDecimal b, boolean expected, String prefix)411     private static int compare(BigDecimal a, BigDecimal b, boolean expected, String prefix) {
412         boolean result = a.equals(b);
413         int failed = (result==expected) ? 0 : 1;
414         if (failed == 1) {
415             System.err.println("Testing " + prefix +
416                                "(" + a + ").compareTo(" + b + ") => " + result +
417                                "\n\tExpected " + expected);
418         }
419         return failed;
420     }
421 
equalNumerically(BigDecimal a, BigDecimal b, String prefix)422     private static int equalNumerically(BigDecimal a, BigDecimal b,
423                                         String prefix) {
424         return compareNumerically(a, b, 0, prefix);
425     }
426 
427 
compareNumerically(BigDecimal a, BigDecimal b, int expected, String prefix)428     private static int compareNumerically(BigDecimal a, BigDecimal b,
429                                           int expected, String prefix) {
430         int result = a.compareTo(b);
431         int failed = (result==expected) ? 0 : 1;
432         if (failed == 1) {
433             System.err.println("Testing " + prefix +
434                                "(" + a + ").compareTo(" + b + ") => " + result +
435                                "\n\tExpected " + expected);
436         }
437         return failed;
438     }
439 
440     /**
441      * Alternative implementation of BigDecimal square root which uses
442      * higher-precision for a simpler set of termination conditions
443      * for the Newton iteration.
444      */
445     private static class BigSquareRoot {
446 
447         /**
448          * The value 0.5, with a scale of 1.
449          */
450         private static final BigDecimal ONE_HALF = valueOf(5L, 1);
451 
isPowerOfTen(BigDecimal bd)452         public static boolean isPowerOfTen(BigDecimal bd) {
453             return BigInteger.ONE.equals(bd.unscaledValue());
454         }
455 
square(BigDecimal bd)456         public static BigDecimal square(BigDecimal bd) {
457             return bd.multiply(bd);
458         }
459 
sqrt(BigDecimal bd, MathContext mc)460         public static BigDecimal sqrt(BigDecimal bd, MathContext mc) {
461             int signum = bd.signum();
462             if (signum == 1) {
463                 /*
464                  * The following code draws on the algorithm presented in
465                  * "Properly Rounded Variable Precision Square Root," Hull and
466                  * Abrham, ACM Transactions on Mathematical Software, Vol 11,
467                  * No. 3, September 1985, Pages 229-237.
468                  *
469                  * The BigDecimal computational model differs from the one
470                  * presented in the paper in several ways: first BigDecimal
471                  * numbers aren't necessarily normalized, second many more
472                  * rounding modes are supported, including UNNECESSARY, and
473                  * exact results can be requested.
474                  *
475                  * The main steps of the algorithm below are as follows,
476                  * first argument reduce the value to the numerical range
477                  * [1, 10) using the following relations:
478                  *
479                  * x = y * 10 ^ exp
480                  * sqrt(x) = sqrt(y) * 10^(exp / 2) if exp is even
481                  * sqrt(x) = sqrt(y/10) * 10 ^((exp+1)/2) is exp is odd
482                  *
483                  * Then use Newton's iteration on the reduced value to compute
484                  * the numerical digits of the desired result.
485                  *
486                  * Finally, scale back to the desired exponent range and
487                  * perform any adjustment to get the preferred scale in the
488                  * representation.
489                  */
490 
491                 // The code below favors relative simplicity over checking
492                 // for special cases that could run faster.
493 
494                 int preferredScale = bd.scale()/2;
495                 BigDecimal zeroWithFinalPreferredScale =
496                     BigDecimal.valueOf(0L, preferredScale);
497 
498                 // First phase of numerical normalization, strip trailing
499                 // zeros and check for even powers of 10.
500                 BigDecimal stripped = bd.stripTrailingZeros();
501                 int strippedScale = stripped.scale();
502 
503                 // Numerically sqrt(10^2N) = 10^N
504                 if (isPowerOfTen(stripped) &&
505                     strippedScale % 2 == 0) {
506                     BigDecimal result = BigDecimal.valueOf(1L, strippedScale/2);
507                     if (result.scale() != preferredScale) {
508                         // Adjust to requested precision and preferred
509                         // scale as appropriate.
510                         result = result.add(zeroWithFinalPreferredScale, mc);
511                     }
512                     return result;
513                 }
514 
515                 // After stripTrailingZeros, the representation is normalized as
516                 //
517                 // unscaledValue * 10^(-scale)
518                 //
519                 // where unscaledValue is an integer with the mimimum
520                 // precision for the cohort of the numerical value. To
521                 // allow binary floating-point hardware to be used to get
522                 // approximately a 15 digit approximation to the square
523                 // root, it is helpful to instead normalize this so that
524                 // the significand portion is to right of the decimal
525                 // point by roughly (scale() - precision() + 1).
526 
527                 // Now the precision / scale adjustment
528                 int scaleAdjust = 0;
529                 int scale = stripped.scale() - stripped.precision() + 1;
530                 if (scale % 2 == 0) {
531                     scaleAdjust = scale;
532                 } else {
533                     scaleAdjust = scale - 1;
534                 }
535 
536                 BigDecimal working = stripped.scaleByPowerOfTen(scaleAdjust);
537 
538                 assert  // Verify 0.1 <= working < 10
539                     ONE_TENTH.compareTo(working) <= 0 && working.compareTo(TEN) < 0;
540 
541                 // Use good ole' Math.sqrt to get the initial guess for
542                 // the Newton iteration, good to at least 15 decimal
543                 // digits. This approach does incur the cost of a
544                 //
545                 // BigDecimal -> double -> BigDecimal
546                 //
547                 // conversion cycle, but it avoids the need for several
548                 // Newton iterations in BigDecimal arithmetic to get the
549                 // working answer to 15 digits of precision. If many fewer
550                 // than 15 digits were needed, it might be faster to do
551                 // the loop entirely in BigDecimal arithmetic.
552                 //
553                 // (A double value might have as much many as 17 decimal
554                 // digits of precision; it depends on the relative density
555                 // of binary and decimal numbers at different regions of
556                 // the number line.)
557                 //
558                 // (It would be possible to check for certain special
559                 // cases to avoid doing any Newton iterations. For
560                 // example, if the BigDecimal -> double conversion was
561                 // known to be exact and the rounding mode had a
562                 // low-enough precision, the post-Newton rounding logic
563                 // could be applied directly.)
564 
565                 BigDecimal guess = new BigDecimal(Math.sqrt(working.doubleValue()));
566                 int guessPrecision = 15;
567                 int originalPrecision = mc.getPrecision();
568                 int targetPrecision;
569 
570                 // If an exact value is requested, it must only need
571                 // about half of the input digits to represent since
572                 // multiplying an N digit number by itself yield a (2N
573                 // - 1) digit or 2N digit result.
574                 if (originalPrecision == 0) {
575                     targetPrecision = stripped.precision()/2 + 1;
576                 } else {
577                     targetPrecision = originalPrecision;
578                 }
579 
580                 // When setting the precision to use inside the Newton
581                 // iteration loop, take care to avoid the case where the
582                 // precision of the input exceeds the requested precision
583                 // and rounding the input value too soon.
584                 BigDecimal approx = guess;
585                 int workingPrecision = working.precision();
586                 // Use "2p + 2" property to guarantee enough
587                 // intermediate precision so that a double-rounding
588                 // error does not occur when rounded to the final
589                 // destination precision.
590                 int loopPrecision =
591                     Math.max(2 * Math.max(targetPrecision, workingPrecision) + 2,
592                              34); // Force at least two Netwon
593                                   // iterations on the Math.sqrt
594                                   // result.
595                 do {
596                     MathContext mcTmp = new MathContext(loopPrecision, RoundingMode.HALF_EVEN);
597                     // approx = 0.5 * (approx + fraction / approx)
598                     approx = ONE_HALF.multiply(approx.add(working.divide(approx, mcTmp), mcTmp));
599                     guessPrecision *= 2;
600                 } while (guessPrecision < loopPrecision);
601 
602                 BigDecimal result;
603                 RoundingMode targetRm = mc.getRoundingMode();
604                 if (targetRm == RoundingMode.UNNECESSARY || originalPrecision == 0) {
605                     RoundingMode tmpRm =
606                         (targetRm == RoundingMode.UNNECESSARY) ? RoundingMode.DOWN : targetRm;
607                     MathContext mcTmp = new MathContext(targetPrecision, tmpRm);
608                     result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mcTmp);
609 
610                     // If result*result != this numerically, the square
611                     // root isn't exact
612                     if (bd.subtract(square(result)).compareTo(ZERO) != 0) {
613                         throw new ArithmeticException("Computed square root not exact.");
614                     }
615                 } else {
616                     result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mc);
617                 }
618 
619                 assert squareRootResultAssertions(bd, result, mc);
620                 if (result.scale() != preferredScale) {
621                     // The preferred scale of an add is
622                     // max(addend.scale(), augend.scale()). Therefore, if
623                     // the scale of the result is first minimized using
624                     // stripTrailingZeros(), adding a zero of the
625                     // preferred scale rounding the correct precision will
626                     // perform the proper scale vs precision tradeoffs.
627                     result = result.stripTrailingZeros().
628                         add(zeroWithFinalPreferredScale,
629                             new MathContext(originalPrecision, RoundingMode.UNNECESSARY));
630                 }
631                 return result;
632             } else {
633                 switch (signum) {
634                 case -1:
635                     throw new ArithmeticException("Attempted square root " +
636                                                   "of negative BigDecimal");
637                 case 0:
638                     return valueOf(0L, bd.scale()/2);
639 
640                 default:
641                     throw new AssertionError("Bad value from signum");
642                 }
643             }
644         }
645 
646         /**
647          * For nonzero values, check numerical correctness properties of
648          * the computed result for the chosen rounding mode.
649          *
650          * For the directed roundings, for DOWN and FLOOR, result^2 must
651          * be {@code <=} the input and (result+ulp)^2 must be {@code >} the
652          * input. Conversely, for UP and CEIL, result^2 must be {@code >=} the
653          * input and (result-ulp)^2 must be {@code <} the input.
654          */
655         private static boolean squareRootResultAssertions(BigDecimal input, BigDecimal result, MathContext mc) {
656             if (result.signum() == 0) {
657                 return squareRootZeroResultAssertions(input, result, mc);
658             } else {
659                 RoundingMode rm = mc.getRoundingMode();
660                 BigDecimal ulp = result.ulp();
661                 BigDecimal neighborUp   = result.add(ulp);
662                 // Make neighbor down accurate even for powers of ten
663                 if (isPowerOfTen(result)) {
664                     ulp = ulp.divide(TEN);
665                 }
666                 BigDecimal neighborDown = result.subtract(ulp);
667 
668                 // Both the starting value and result should be nonzero and positive.
669                 if (result.signum() != 1 ||
670                     input.signum() != 1) {
671                     return false;
672                 }
673 
674                 switch (rm) {
675                 case DOWN:
676                 case FLOOR:
677                     assert
678                         square(result).compareTo(input)    <= 0 &&
679                         square(neighborUp).compareTo(input) > 0:
680                     "Square of result out for bounds rounding " + rm;
681                     return true;
682 
683                 case UP:
684                 case CEILING:
685                     assert
686                         square(result).compareTo(input) >= 0 :
687                     "Square of result too small rounding " + rm;
688 
689                     assert
690                         square(neighborDown).compareTo(input) < 0 :
691                     "Square of down neighbor too large rounding  " + rm + "\n" +
692                         "\t input: " + input + "\t neighborDown: " +  neighborDown +"\t sqrt: " + result +
693                         "\t" + mc;
694                     return true;
695 
696 
697                 case HALF_DOWN:
698                 case HALF_EVEN:
699                 case HALF_UP:
700                     BigDecimal err = square(result).subtract(input).abs();
701                     BigDecimal errUp = square(neighborUp).subtract(input);
702                     BigDecimal errDown =  input.subtract(square(neighborDown));
703                     // All error values should be positive so don't need to
704                     // compare absolute values.
705 
706                     int err_comp_errUp = err.compareTo(errUp);
707                     int err_comp_errDown = err.compareTo(errDown);
708 
709                     assert
710                         errUp.signum()   == 1 &&
711                         errDown.signum() == 1 :
712                     "Errors of neighbors squared don't have correct signs";
713 
714                     // At least one of these must be true, but not both
715 //                     assert
716 //                         err_comp_errUp   <= 0 : "Upper neighbor is closer than result: " + rm +
717 //                         "\t" + input + "\t result" + result;
718 //                     assert
719 //                         err_comp_errDown <= 0 : "Lower neighbor is closer than result: " + rm +
720 //                         "\t" + input + "\t result " + result + "\t lower neighbor: " + neighborDown;
721 
722                     assert
723                         ((err_comp_errUp   == 0 ) ? err_comp_errDown < 0 : true) &&
724                         ((err_comp_errDown == 0 ) ? err_comp_errUp   < 0 : true) :
725                             "Incorrect error relationships";
726                         // && could check for digit conditions for ties too
727                         return true;
728 
729                 default: // Definition of UNNECESSARY already verified.
730                     return true;
731                 }
732             }
733         }
734 
735         private static boolean squareRootZeroResultAssertions(BigDecimal input,
736                                                               BigDecimal result,
737                                                               MathContext mc) {
738             return input.compareTo(ZERO) == 0;
739         }
740     }
741 }
742 
743