1 /*
2 * Copyright (c) 2016, 2019, 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.TEN;
34 import static java.math.BigDecimal.ZERO;
35 import static java.math.BigDecimal.valueOf;
36
37 public class SquareRootTests {
38 private static BigDecimal TWO = new BigDecimal(2);
39
40 public static void main(String... args) {
41 int failures = 0;
42
43 failures += negativeTests();
44 failures += zeroTests();
45 failures += evenPowersOfTenTests();
46 failures += squareRootTwoTests();
47 failures += lowPrecisionPerfectSquares();
48 failures += almostFourRoundingDown();
49 failures += almostFourRoundingUp();
50 failures += nearTen();
51 failures += nearOne();
52 failures += halfWay();
53
54 if (failures > 0 ) {
55 throw new RuntimeException("Incurred " + failures + " failures" +
56 " testing BigDecimal.sqrt().");
57 }
58 }
59
60 private static int negativeTests() {
61 int failures = 0;
62
63 for (long i = -10; i < 0; i++) {
64 for (int j = -5; j < 5; j++) {
65 try {
66 BigDecimal input = BigDecimal.valueOf(i, j);
67 BigDecimal result = input.sqrt(MathContext.DECIMAL64);
68 System.err.println("Unexpected sqrt of negative: (" +
69 input + ").sqrt() = " + result );
70 failures += 1;
71 } catch (ArithmeticException e) {
72 ; // Expected
73 }
74 }
75 }
76
77 return failures;
78 }
79
80 private static int zeroTests() {
81 int failures = 0;
82
83 for (int i = -100; i < 100; i++) {
84 BigDecimal expected = BigDecimal.valueOf(0L, i/2);
85 // These results are independent of rounding mode
86 failures += compare(BigDecimal.valueOf(0L, i).sqrt(MathContext.UNLIMITED),
87 expected, true, "zeros");
88
89 failures += compare(BigDecimal.valueOf(0L, i).sqrt(MathContext.DECIMAL64),
90 expected, true, "zeros");
91 }
92
93 return failures;
94 }
95
96 /**
97 * sqrt(10^2N) is 10^N
98 * Both numerical value and representation should be verified
99 */
100 private static int evenPowersOfTenTests() {
101 int failures = 0;
102 MathContext oneDigitExactly = new MathContext(1, RoundingMode.UNNECESSARY);
103
104 for (int scale = -100; scale <= 100; scale++) {
105 BigDecimal testValue = BigDecimal.valueOf(1, 2*scale);
106 BigDecimal expectedNumericalResult = BigDecimal.valueOf(1, scale);
107
108 BigDecimal result;
109
110 failures += equalNumerically(expectedNumericalResult,
111 result = testValue.sqrt(MathContext.DECIMAL64),
112 "Even powers of 10, DECIMAL64");
113
114 // Can round to one digit of precision exactly
115 failures += equalNumerically(expectedNumericalResult,
116 result = testValue.sqrt(oneDigitExactly),
117 "even powers of 10, 1 digit");
118 if (result.precision() > 1) {
119 failures += 1;
120 System.err.println("Excess precision for " + result);
121 }
122
123 // If rounding to more than one digit, do precision / scale checking...
124 }
125
126 return failures;
127 }
128
129 private static int squareRootTwoTests() {
130 int failures = 0;
131
132 // Square root of 2 truncated to 65 digits
133 BigDecimal highPrecisionRoot2 =
134 new BigDecimal("1.41421356237309504880168872420969807856967187537694807317667973799");
135
136 RoundingMode[] modes = {
137 RoundingMode.UP, RoundingMode.DOWN,
138 RoundingMode.CEILING, RoundingMode.FLOOR,
139 RoundingMode.HALF_UP, RoundingMode.HALF_DOWN, RoundingMode.HALF_EVEN
140 };
141
142
143 // For each interesting rounding mode, for precisions 1 to, say,
144 // 63 numerically compare TWO.sqrt(mc) to
145 // highPrecisionRoot2.round(mc) and the alternative internal high-precision
146 // implementation of square root.
147 for (RoundingMode mode : modes) {
148 for (int precision = 1; precision < 63; precision++) {
149 MathContext mc = new MathContext(precision, mode);
150 BigDecimal expected = highPrecisionRoot2.round(mc);
151 BigDecimal computed = TWO.sqrt(mc);
152 BigDecimal altComputed = BigSquareRoot.sqrt(TWO, mc);
153
154 failures += equalNumerically(expected, computed, "sqrt(2)");
155 failures += equalNumerically(computed, altComputed, "computed & altComputed");
156 }
157 }
158
159 return failures;
160 }
161
162 private static int lowPrecisionPerfectSquares() {
163 int failures = 0;
164
165 // For 5^2 through 9^2, if the input is rounded to one digit
166 // first before the root is computed, the wrong answer will
167 // result. Verify results and scale for different rounding
168 // modes and precisions.
169 long[][] squaresWithOneDigitRoot = {{ 4, 2},
170 { 9, 3},
171 {25, 5},
172 {36, 6},
173 {49, 7},
174 {64, 8},
175 {81, 9}};
176
177 for (long[] squareAndRoot : squaresWithOneDigitRoot) {
178 BigDecimal square = new BigDecimal(squareAndRoot[0]);
179 BigDecimal expected = new BigDecimal(squareAndRoot[1]);
180
181 for (int scale = 0; scale <= 4; scale++) {
182 BigDecimal scaledSquare = square.setScale(scale, RoundingMode.UNNECESSARY);
183 int expectedScale = scale/2;
184 for (int precision = 0; precision <= 5; precision++) {
185 for (RoundingMode rm : RoundingMode.values()) {
186 MathContext mc = new MathContext(precision, rm);
187 BigDecimal computedRoot = scaledSquare.sqrt(mc);
188 failures += equalNumerically(expected, computedRoot, "simple squares");
189 int computedScale = computedRoot.scale();
190 if (precision >= expectedScale + 1 &&
191 computedScale != expectedScale) {
192 System.err.printf("%s\tprecision=%d\trm=%s%n",
193 computedRoot.toString(), precision, rm);
194 failures++;
195 System.err.printf("\t%s does not have expected scale of %d%n.",
196 computedRoot, expectedScale);
197 }
198 }
199 }
200 }
201 }
202
203 return failures;
204 }
205
206 /**
207 * Test around 3.9999 that the result doesn't not improperly
208 * round-up to a numerical value of 2.
209 */
210 private static int almostFourRoundingDown() {
211 int failures = 0;
212 BigDecimal nearFour = new BigDecimal("3.999999999999999999999999999999");
213
214 // Sqrt root is 1.9999...
215
216 for (int i = 1; i < 64; i++) {
217 MathContext mc = new MathContext(i, RoundingMode.FLOOR);
218 BigDecimal result = nearFour.sqrt(mc);
219 BigDecimal expected = BigSquareRoot.sqrt(nearFour, mc);
220 failures += equalNumerically(expected, result, "near four rounding down");
221 failures += (result.compareTo(TWO) < 0) ? 0 : 1 ;
222 }
223
224 return failures;
225 }
226
227 /**
228 * Test around 4.000...1 that the result doesn't not improperly
229 * round-down to a numerical value of 2.
230 */
231 private static int almostFourRoundingUp() {
232 int failures = 0;
233 BigDecimal nearFour = new BigDecimal("4.000000000000000000000000000001");
234
235 // Sqrt root is 2.0000....<non-zero digits>
236
237 for (int i = 1; i < 64; i++) {
238 MathContext mc = new MathContext(i, RoundingMode.CEILING);
239 BigDecimal result = nearFour.sqrt(mc);
240 BigDecimal expected = BigSquareRoot.sqrt(nearFour, mc);
241 failures += equalNumerically(expected, result, "near four rounding down");
242 failures += (result.compareTo(TWO) > 0) ? 0 : 1 ;
243 }
244
245 return failures;
246 }
247
248 private static int nearTen() {
249 int failures = 0;
250
251 BigDecimal near10 = new BigDecimal("9.99999999999999999999");
252
253 BigDecimal near10sq = near10.multiply(near10);
254
255 BigDecimal near10sq_ulp = near10sq.add(near10sq.ulp());
256
257 for (int i = 10; i < 23; i++) {
258 MathContext mc = new MathContext(i, RoundingMode.HALF_EVEN);
259
260 failures += equalNumerically(BigSquareRoot.sqrt(near10sq_ulp, mc),
261 near10sq_ulp.sqrt(mc),
262 "near 10 rounding down");
263 }
264
265 return failures;
266 }
267
268
269 /*
270 * Probe for rounding failures near a power of ten, 1 = 10^0,
271 * where an ulp has a different size above and below the value.
272 */
273 private static int nearOne() {
274 int failures = 0;
275
276 BigDecimal near1 = new BigDecimal(".999999999999999999999");
277 BigDecimal near1sq = near1.multiply(near1);
278 BigDecimal near1sq_ulp = near1sq.add(near1sq.ulp());
279
280 for (int i = 10; i < 23; i++) {
281 for (RoundingMode rm : List.of(RoundingMode.HALF_EVEN,
282 RoundingMode.UP,
283 RoundingMode.DOWN )) {
284 MathContext mc = new MathContext(i, rm);
285 failures += equalNumerically(BigSquareRoot.sqrt(near1sq_ulp, mc),
286 near1sq_ulp.sqrt(mc),
287 "near 1 half even");
288 }
289 }
290
291 return failures;
292 }
293
294
295 private static int halfWay() {
296 int failures = 0;
297
298 /*
299 * Use enough digits that the exact result cannot be computed
300 * from the sqrt of a double.
301 */
302 BigDecimal[] halfWayCases = {
303 // Odd next digit, truncate on HALF_EVEN
304 new BigDecimal("123456789123456789.5"),
305
306 // Even next digit, round up on HALF_EVEN
307 new BigDecimal("123456789123456788.5"),
308 };
309
310 for (BigDecimal halfWayCase : halfWayCases) {
311 // Round result to next-to-last place
312 int precision = halfWayCase.precision() - 1;
313 BigDecimal square = halfWayCase.multiply(halfWayCase);
314
315 for (RoundingMode rm : List.of(RoundingMode.HALF_EVEN,
316 RoundingMode.HALF_UP,
317 RoundingMode.HALF_DOWN)) {
318 MathContext mc = new MathContext(precision, rm);
319
320 System.out.println("\nRounding mode " + rm);
321 System.out.println("\t" + halfWayCase.round(mc) + "\t" + halfWayCase);
322 /*System.out.println("\t" + square.sqrt(mc));*/
323 System.out.println("\t" + BigSquareRoot.sqrt(square, mc));
324
325 failures += equalNumerically(/*square.sqrt(mc),*/
326 BigSquareRoot.sqrt(square, mc),
327 halfWayCase.round(mc),
328 "Rounding halway " + rm);
329 }
330 }
331
332 return failures;
333 }
334
335 private static int compare(BigDecimal a, BigDecimal b, boolean expected, String prefix) {
336 boolean result = a.equals(b);
337 int failed = (result==expected) ? 0 : 1;
338 if (failed == 1) {
339 System.err.println("Testing " + prefix +
340 "(" + a + ").compareTo(" + b + ") => " + result +
341 "\n\tExpected " + expected);
342 }
343 return failed;
344 }
345
346 private static int equalNumerically(BigDecimal a, BigDecimal b,
347 String prefix) {
348 return compareNumerically(a, b, 0, prefix);
349 }
350
351
352 private static int compareNumerically(BigDecimal a, BigDecimal b,
353 int expected, String prefix) {
354 int result = a.compareTo(b);
355 int failed = (result==expected) ? 0 : 1;
356 if (failed == 1) {
357 System.err.println("Testing " + prefix +
358 "(" + a + ").compareTo(" + b + ") => " + result +
359 "\n\tExpected " + expected);
360 }
361 return failed;
362 }
363
364 /**
365 * Alternative implementation of BigDecimal square root which uses
366 * higher-precision for a simpler set of termination conditions
367 * for the Newton iteration.
368 */
369 private static class BigSquareRoot {
370 /**
371 * The value 0.1, with a scale of 1.
372 */
373 private static final BigDecimal ONE_TENTH = valueOf(1L, 1);
374
375 /**
376 * The value 0.5, with a scale of 1.
377 */
378 private static final BigDecimal ONE_HALF = valueOf(5L, 1);
379
380 private static boolean isPowerOfTen(BigDecimal bd) {
381 return BigInteger.ONE.equals(bd.unscaledValue());
382 }
383
384 public static BigDecimal sqrt(BigDecimal bd, MathContext mc) {
385 int signum = bd.signum();
386 if (signum == 1) {
387 /*
388 * The following code draws on the algorithm presented in
389 * "Properly Rounded Variable Precision Square Root," Hull and
390 * Abrham, ACM Transactions on Mathematical Software, Vol 11,
391 * No. 3, September 1985, Pages 229-237.
392 *
393 * The BigDecimal computational model differs from the one
394 * presented in the paper in several ways: first BigDecimal
395 * numbers aren't necessarily normalized, second many more
396 * rounding modes are supported, including UNNECESSARY, and
397 * exact results can be requested.
398 *
399 * The main steps of the algorithm below are as follows,
400 * first argument reduce the value to the numerical range
401 * [1, 10) using the following relations:
402 *
403 * x = y * 10 ^ exp
404 * sqrt(x) = sqrt(y) * 10^(exp / 2) if exp is even
405 * sqrt(x) = sqrt(y/10) * 10 ^((exp+1)/2) is exp is odd
406 *
407 * Then use Newton's iteration on the reduced value to compute
408 * the numerical digits of the desired result.
409 *
410 * Finally, scale back to the desired exponent range and
411 * perform any adjustment to get the preferred scale in the
412 * representation.
413 */
414
415 // The code below favors relative simplicity over checking
416 // for special cases that could run faster.
417
418 int preferredScale = bd.scale()/2;
419 BigDecimal zeroWithFinalPreferredScale =
420 BigDecimal.valueOf(0L, preferredScale);
421
422 // First phase of numerical normalization, strip trailing
423 // zeros and check for even powers of 10.
424 BigDecimal stripped = bd.stripTrailingZeros();
425 int strippedScale = stripped.scale();
426
427 // Numerically sqrt(10^2N) = 10^N
428 if (isPowerOfTen(stripped) &&
429 strippedScale % 2 == 0) {
430 BigDecimal result = BigDecimal.valueOf(1L, strippedScale/2);
431 if (result.scale() != preferredScale) {
432 // Adjust to requested precision and preferred
433 // scale as appropriate.
434 result = result.add(zeroWithFinalPreferredScale, mc);
435 }
436 return result;
437 }
438
439 // After stripTrailingZeros, the representation is normalized as
440 //
441 // unscaledValue * 10^(-scale)
442 //
443 // where unscaledValue is an integer with the mimimum
444 // precision for the cohort of the numerical value. To
445 // allow binary floating-point hardware to be used to get
446 // approximately a 15 digit approximation to the square
447 // root, it is helpful to instead normalize this so that
448 // the significand portion is to right of the decimal
449 // point by roughly (scale() - precision() +1).
450
451 // Now the precision / scale adjustment
452 int scaleAdjust = 0;
453 int scale = stripped.scale() - stripped.precision() + 1;
454 if (scale % 2 == 0) {
455 scaleAdjust = scale;
456 } else {
457 scaleAdjust = scale - 1;
458 }
459
460 BigDecimal working = stripped.scaleByPowerOfTen(scaleAdjust);
461
462 assert // Verify 0.1 <= working < 10
463 ONE_TENTH.compareTo(working) <= 0 && working.compareTo(TEN) < 0;
464
465 // Use good ole' Math.sqrt to get the initial guess for
466 // the Newton iteration, good to at least 15 decimal
467 // digits. This approach does incur the cost of a
468 //
469 // BigDecimal -> double -> BigDecimal
470 //
471 // conversion cycle, but it avoids the need for several
472 // Newton iterations in BigDecimal arithmetic to get the
473 // working answer to 15 digits of precision. If many fewer
474 // than 15 digits were needed, it might be faster to do
475 // the loop entirely in BigDecimal arithmetic.
476 //
477 // (A double value might have as much many as 17 decimal
478 // digits of precision; it depends on the relative density
479 // of binary and decimal numbers at different regions of
480 // the number line.)
481 //
482 // (It would be possible to check for certain special
483 // cases to avoid doing any Newton iterations. For
484 // example, if the BigDecimal -> double conversion was
485 // known to be exact and the rounding mode had a
486 // low-enough precision, the post-Newton rounding logic
487 // could be applied directly.)
488
489 BigDecimal guess = new BigDecimal(Math.sqrt(working.doubleValue()));
490 int guessPrecision = 15;
491 int originalPrecision = mc.getPrecision();
492 int targetPrecision;
493
494 // If an exact value is requested, it must only need about
495 // half of the input digits to represent since multiplying
496 // an N digit number by itself yield a 2N-1 digit or 2N
497 // digit result.
498 if (originalPrecision == 0) {
499 targetPrecision = stripped.precision()/2 + 1;
500 } else {
501 targetPrecision = originalPrecision;
502 }
503
504 // When setting the precision to use inside the Newton
505 // iteration loop, take care to avoid the case where the
506 // precision of the input exceeds the requested precision
507 // and rounding the input value too soon.
508 BigDecimal approx = guess;
509 int workingPrecision = working.precision();
510 // Use "2p + 2" property to guarantee enough
511 // intermediate precision so that a double-rounding
512 // error does not occur when rounded to the final
513 // destination precision.
514 int loopPrecision = Math.max(Math.max(2 * targetPrecision + 2,
515 workingPrecision),
516 34); // Force at least
517 // two Netwon
518 // iterations on the
519 // Math.sqrt result.
520 do {
521 int tmpPrecision = Math.max(Math.max(guessPrecision, targetPrecision + 2),
522 workingPrecision);
523 MathContext mcTmp = new MathContext(loopPrecision, RoundingMode.HALF_EVEN);
524 // approx = 0.5 * (approx + fraction / approx)
525 approx = ONE_HALF.multiply(approx.add(working.divide(approx, mcTmp), mcTmp));
526 guessPrecision *= 2;
527 } while (guessPrecision < loopPrecision);
528
529 BigDecimal result;
530 RoundingMode targetRm = mc.getRoundingMode();
531 if (targetRm == RoundingMode.UNNECESSARY || originalPrecision == 0) {
532 RoundingMode tmpRm =
533 (targetRm == RoundingMode.UNNECESSARY) ? RoundingMode.DOWN : targetRm;
534 MathContext mcTmp = new MathContext(targetPrecision, tmpRm);
535 result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mcTmp);
536
537 // If result*result != this numerically, the square
538 // root isn't exact
539 if (bd.subtract(result.multiply(result)).compareTo(ZERO) != 0) {
540 throw new ArithmeticException("Computed square root not exact.");
541 }
542 } else {
543 result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mc);
544 }
545
546 if (result.scale() != preferredScale) {
547 // The preferred scale of an add is
548 // max(addend.scale(), augend.scale()). Therefore, if
549 // the scale of the result is first minimized using
550 // stripTrailingZeros(), adding a zero of the
551 // preferred scale rounding the correct precision will
552 // perform the proper scale vs precision tradeoffs.
553 result = result.stripTrailingZeros().
554 add(zeroWithFinalPreferredScale,
555 new MathContext(originalPrecision, RoundingMode.UNNECESSARY));
556 }
557 return result;
558 } else {
559 switch (signum) {
560 case -1:
561 throw new ArithmeticException("Attempted square root " +
562 "of negative BigDecimal");
563 case 0:
564 return valueOf(0L, bd.scale()/2);
565
566 default:
567 throw new AssertionError("Bad value from signum");
568 }
569 }
570 }
571 }
572 }