--- old/src/java.base/share/classes/java/math/BigDecimal.java 2016-05-11 10:37:28.126558223 -0700 +++ new/src/java.base/share/classes/java/math/BigDecimal.java 2016-05-11 10:37:27.730558229 -0700 @@ -346,6 +346,16 @@ public static final BigDecimal TEN = ZERO_THROUGH_TEN[10]; + /** + * The value 0.1, with a scale of 1. + */ + private static final BigDecimal ONE_TENTH = valueOf(1L, 1); + + /** + * The value 0.5, with a scale of 1. + */ + private static final BigDecimal ONE_HALF = valueOf(5L, 1); + // Constructors /** @@ -1995,6 +2005,219 @@ return result; } + + /** + * Returns an approximation to the square root of {@code this} + * with rounding according to the context settings. + * + *

The preferred scale of the returned result is equal to + * {@code floor(this.scale()/2.0)}. The value of the returned + * result is always within one ulp of the exact decimal value for + * the precision in question. If the rounding mode is {@link + * RoundingMode#HALF_UP HALF_UP}, {@link RoundingMode#HALF_DOWN + * HALF_DOWN}, or {@link RoundingMode#HALF_EVEN HALF_EVEN}, the + * result is within one half an ulp of the exact decimal value. + * + *

Special case: + *

+ * + * @param mc the context to use. + * @return the square root of {@code this}. + * @throws ArithmeticException if {@code this} is less than zero. + * @throws ArithmeticException if an exact result is requested + * ({@code mc.getPrecision()==0}) and there is no finite decimal + * expansion of the exact result + * @throws ArithmeticException if + * {@code (mc.getRoundingMode()==RoundingMode.UNNECESSARY}) and + * the exact result cannot fit in {@code mc.getPrecision()} + * digits. + * @since 9 + */ + public BigDecimal sqrt(MathContext mc) { + int signum = 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 follow, 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. + + BigDecimal zeroWithFinalPreferredScale = valueOf(0L, this.scale()/2); + + // First phase of numerical normalization, strip trailing + // zeros and check for even powers of 10. + BigDecimal stripped = this.stripTrailingZeros(); + int strippedScale = stripped.scale(); + + // Numerically, sqrt(10^2N) = 10^N + if (BigInteger.ONE.equals(stripped.unscaledValue()) && + strippedScale % 2 == 0) { + BigDecimal result = valueOf(1L, strippedScale/2); + // Adjust to requested precision and preferred scale as possible + return result.add(zeroWithFinalPreferredScale, mc); + } + + System.out.println("\tStarting: " + this); // DEBUG + System.out.println("\tStripped: " + stripped); // DEBUG + + + // 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 as so + // that the significand portion is to right of the decimal + // point, roughly: + // + // (unscaledValue * (10^-precision) * 10^(-scale)) * (10^precision) + // + // so that + // + // sqrt(unscaledValue * (10^-precision) * 10^(-scale) * (10^precision)) = + // + // sqrt(unscaledValue * (10^-precision) * 10^(-scale)) * 10^(precision/2) + // + // Therefore, this adjustment occurs for by 10^-precision is precision is even or + // (adjust as needed, +/-1 + // + // (A double value might have as much as 17 decimal digits + // of precision; it depends on the relative density of + // binary and decimal numbers at different points of the + // number line.) + + // 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; + stripped = stripped.divide(TEN); + } + + // At least document potential exponent overflow issues here??? + BigDecimal working = stripped.scaleByPowerOfTen(scaleAdjust); + System.out.println("\tWorking2: " + working + "\tscale: " + working.scale()); // DEBUG + + // Verify 1 <= working < 10 (verify range) + assert working.compareTo(ONE) >= 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. + BigDecimal guess = new BigDecimal(Math.sqrt(working.doubleValue())); + int guessPrecision = 15; + int targetPrecision = mc.getPrecision(); + RoundingMode rm = mc.getRoundingMode(); + + // When setting the precision to use inside the 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(); + do { + System.out.println(approx); // DEBUG + int tmpPrecision = Math.max(Math.max(guessPrecision, targetPrecision + 2), + workingPrecision); + MathContext mcTmp = new MathContext(tmpPrecision, RoundingMode.HALF_EVEN); + // approx = 0.5 * (approx + fraction / approx) + approx = ONE_HALF.multiply(approx.add(working.divide(approx, mcTmp), mcTmp)); + guessPrecision *= 2; + } while (guessPrecision < targetPrecision + 2); + + + BigDecimal result = approx.round(mc); + + // TODO: need code for Round.UNNECESSARY. Probably also + // need additinal checking for HALF_EVEN since we are only + // using two extra precision digits. + + // For now, skip the rounding up / down fix up + + // Read through "What every ..." for 2p + 2 discussion... + + // Handle 1/2 ulp, 1 ulp, and unnecessary rounding separately? + // UP, DOWN, + // CEILING (equiv to UP), FLOOR (equiv to DOWN), + // + // HALF_UP, HALF_DOWN, HALF_EVEN -- look up 2p + 2 derivation + // -- same for decimal as well as binary? + + // Could use direct definition to test adjacent values to + // target precision -- need to lookup if Newton's iteration + // converges from above or below to avoid checking three values... + + // UNNECESSARY + + // Unnecessary rounding -- in the worst case, only need about + // (1/2 + 1) digits to produce an exact output, sqrt(100) = + // 10, sqrt(10_000) = 100, etc.. How to handle cases that are + // powers of 10 with trailing zeros? + + // Normalize to a range [.1, 1] or [1, 10] -- presunmably this + // is open on one end to account for exact square roots... + // + + + // Convert to double and take double sqrt -- easy way to get + // (at least) 15 digits of a decimal approximation :-) + + // If did more work, could make sure came in from uniformly + // above / below the root? + + // Final answer: rescale fraction and then add a zero with the + // preferred scale with a rounding precision equal to the + // target precision. + + result = result.scaleByPowerOfTen(-scaleAdjust/2); + return result.add(zeroWithFinalPreferredScale, + new MathContext(mc.getPrecision(), RoundingMode.UNNECESSARY)); + } else { + switch (signum) { + case -1: + throw new ArithmeticException("Attempted square root " + + "of negative BigDecimal"); + case 0: + return valueOf(0L, scale()/2); + + default: + throw new AssertionError("Bad value from signum"); + } + } + } + /** * Returns a {@code BigDecimal} whose value is * (thisn), The power is computed exactly, to