< prev index next >

src/java.base/share/classes/java/math/BigDecimal.java

Print this page

        

*** 344,353 **** --- 344,358 ---- * @since 1.5 */ 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 /** * Trusted package private constructor. * Trusted simply means if val is INFLATED, intVal could not be null and
*** 1993,2002 **** --- 1998,2228 ---- result[0] = lhs.divideToIntegralValue(divisor, mc); result[1] = lhs.subtract(result[0].multiply(divisor)); return result; } + + /** + * Returns an approximation to the square root of {@code this} + * with rounding according to the context settings. + * + * <p>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. + * + * <p>Special case: + * <ul> + * <li> The square root of a number numerically equal to {@code + * ZERO} is equal to {@code ZERO}. + * </ul> + * + * @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 * <code>(this<sup>n</sup>)</code>, The power is computed exactly, to * unlimited precision. *
< prev index next >