--- old/src/java.base/share/classes/java/math/BigDecimal.java 2016-05-12 14:12:23.154960667 -0700
+++ new/src/java.base/share/classes/java/math/BigDecimal.java 2016-05-12 14:12:22.974960670 -0700
@@ -346,6 +346,11 @@
public static final BigDecimal TEN =
ZERO_THROUGH_TEN[10];
+ /**
+ * The value 0.5, with a scale of 1.
+ */
+ private static final BigDecimal ONE_HALF = valueOf(5L, 1);
+
// Constructors
/**
@@ -1995,6 +2000,227 @@
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:
+ *
+ * - The square root of a number numerically equal to {@code
+ * ZERO} is equal to {@code ZERO}.
+ *
+ *
+ * @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
+ // TODO: since max(addend.scale(), augend.scale()) only add if scale appropriate
+ // Result only has one precision digit as using any precision in mc is okay
+ return result.add(zeroWithFinalPreferredScale, mc);
+ }
+
+ System.out.println("\tStarting: " + this + "\tscale: " + this.scale()); // DEBUG
+ System.out.println("\tStripped: " + stripped + "\tscale: " + stripped.scale()); // 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 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.)
+
+ // 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.multiply(TEN); why not??
+ }
+
+ // At least document potential exponent overflow issues here???
+ BigDecimal working = stripped.scaleByPowerOfTen(scaleAdjust);
+ System.out.println("\tWorking2: " + working + "\tscale: " + working.scale() +
+ "\tadjust: " + scaleAdjust); // 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 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
+ // two N digit numbers yield a 2N-1 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();
+ do {
+ // If did more work, could we make sure the
+ // approximation was uniformly above / below the root?
+
+ 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);
+
+ // Need additinal checking for HALF_EVEN since we are only
+ // using two extra precision digits?
+ 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.round(mcTmp);
+ result.scaleByPowerOfTen(-scaleAdjust/2);
+ // If result*result != the starting value numerically,
+ // the square root isn't exact
+ if (this.subtract(result.multiply(result)).compareTo(ZERO) != 0) {
+ throw new ArithmeticException("Computed square root not exact.");
+ }
+ } else {
+ result = approx.round(mc);
+ result = result.scaleByPowerOfTen(-scaleAdjust/2);
+ }
+
+ // 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...
+
+ // Add a zero with the preferred scale with a rounding
+ // precision equal to the original precision.
+
+ // TODO: since max(addend.scale(), augend.scale()) only add if scale appropriate
+ return result.add(zeroWithFinalPreferredScale,
+ new MathContext(originalPrecision, 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