< prev index next >

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

Print this page

        

*** 126,135 **** --- 126,136 ---- * <tr><th>Operation</th><th>Preferred Scale of Result</th></tr> * <tr><td>Add</td><td>max(addend.scale(), augend.scale())</td> * <tr><td>Subtract</td><td>max(minuend.scale(), subtrahend.scale())</td> * <tr><td>Multiply</td><td>multiplier.scale() + multiplicand.scale()</td> * <tr><td>Divide</td><td>dividend.scale() - divisor.scale()</td> + * <tr><td>Square root</td><td>radicand.scale()/2</td> * </table> * * These scales are the ones used by the methods which return exact * arithmetic results; except that an exact divide may have to use a * larger scale since the exact result may have more digits. For
*** 344,353 **** --- 345,364 ---- * @since 1.5 */ 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 /** * Trusted package private constructor. * Trusted simply means if val is INFLATED, intVal could not be null and
*** 1994,2003 **** --- 2005,2303 ---- 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 this.scale()/2}. 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 numerically equal to {@code ZERO} with a preferred + * scale according to the general rule above. In particular, for + * {@code ZERO}}, {@code ZERO.sqrt(mc).equals(ZERO)} is true with + * any {@code MathContext} as an argument. + * </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. + + int preferredScale = this.scale()/2; + BigDecimal zeroWithFinalPreferredScale = valueOf(0L, preferredScale); + + // 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 (stripped.isPowerOfTen() && + strippedScale % 2 == 0) { + BigDecimal result = valueOf(1L, strippedScale/2); + if (result.scale() != preferredScale) { + // Adjust to requested precision and preferred + // scale as appropriate. + result = result.add(zeroWithFinalPreferredScale, mc); + } + return result; + } + + // 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 by roughly (scale() - precision() +1). + + // 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; + } + + BigDecimal working = stripped.scaleByPowerOfTen(scaleAdjust); + + assert // Verify 0.1 <= working < 10 + ONE_TENTH.compareTo(working) <= 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. This approach does incur the cost of a + // + // BigDecimal -> double -> BigDecimal + // + // conversion cycle, but it avoids the need for several + // Newton iterations in BigDecimal arithmetic to get the + // working answer to 15 digits of precision. If many fewer + // than 15 digits were needed, it might be faster to do + // the loop entirely in BigDecimal arithmetic. + // + // (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.) + // + // (It would be possible to check for certain special + // cases to avoid doing any Newton iterations. For + // example, if the BigDecimal -> double conversion was + // known to be exact and the rounding mode had a + // low-enough precision, the post-Newton rounding logic + // could be applied directly.) + + 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 + // an N digit number by itself yield a 2N-1 digit 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 { + 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; + 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.scaleByPowerOfTen(-scaleAdjust/2).round(mcTmp); + + // If result*result != this 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.scaleByPowerOfTen(-scaleAdjust/2).round(mc); + } + + if (result.scale() != preferredScale) { + // The preferred scale of an add is + // max(addend.scale(), augend.scale()). Therefore, if + // the scale of the result is first minimized using + // stripTrailingZeros(), adding a zero of the + // preferred scale rounding the correct precision will + // perform the proper scale vs precision tradeoffs. + result = result.stripTrailingZeros(). + add(zeroWithFinalPreferredScale, + new MathContext(originalPrecision, RoundingMode.UNNECESSARY)); + } + assert squareRootResultAssertions(result, mc); + return result; + } 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"); + } + } + } + + private boolean isPowerOfTen() { + return BigInteger.ONE.equals(this.unscaledValue()); + } + + /** + * For nonzero values, check numerical correctness properties of + * the computed result for the chosen rounding mode . + * + * For the directed roundings, for DOWN and FLOOR, result^2 must + * be {@code <=} the input and (result+ulp)^2 must be {@code >} the + * input. Conversely, for UP and CEIL, result^2 must be {@code >=} the + * input and (result-ulp)^2 must be {@code <} the input. + */ + private boolean squareRootResultAssertions(BigDecimal result, MathContext mc) { + if (result.signum() == 0) { + return squareRootZeroResultAssertions(result, mc); + } else { + RoundingMode rm = mc.getRoundingMode(); + BigDecimal ulp = result.ulp(); + BigDecimal neighborUp = result.add(ulp); + // Make neighbor down accurate even for powers of ten + if (this.isPowerOfTen()) { + ulp = ulp.divide(TEN); + } + BigDecimal neighborDown = result.subtract(ulp); + + // Both the starting value and result should be nonzero and positive. + if (result.signum() != 1 || + this.signum() != 1) { + return false; + } + + switch (rm) { + case DOWN: + case FLOOR: + return + result.multiply(result).compareTo(this) <= 0 && + neighborUp.multiply(neighborUp).compareTo(this) > 0; + + case UP: + case CEILING: + return + result.multiply(result).compareTo(this) >= 0 && + neighborDown.multiply(neighborDown).compareTo(this) < 0; + + case HALF_DOWN: + case HALF_EVEN: + case HALF_UP: + BigDecimal err = result.multiply(result).subtract(this).abs(); + BigDecimal errUp = neighborUp.multiply(neighborUp).subtract(this); + BigDecimal errDown = this.subtract(neighborDown.multiply(neighborDown)); + // All error values should be positive so don't need to + // compare absolute values. + + int err_comp_errUp = err.compareTo(errUp); + int err_comp_errDown = err.compareTo(errDown); + + return + errUp.signum() == 1 && + errDown.signum() == 1 && + + err_comp_errUp <= 0 && + err_comp_errDown <= 0 && + + ((err_comp_errUp == 0 ) ? err_comp_errDown < 0 : true) && + ((err_comp_errDown == 0 ) ? err_comp_errUp < 0 : true); + // && could check for digit conditions for ties too + + default: // Definition of UNNECESSARY already verified. + return true; + } + } + } + + private boolean squareRootZeroResultAssertions(BigDecimal result, MathContext mc) { + return this.compareTo(ZERO) == 0; + } + + /** * Returns a {@code BigDecimal} whose value is * <code>(this<sup>n</sup>)</code>, The power is computed exactly, to * unlimited precision. * * <p>The parameter {@code n} must be in the range 0 through
< prev index next >