2187 // approx = 0.5 * (approx + fraction / approx)
2188 approx = ONE_HALF.multiply(approx.add(working.divide(approx, mcTmp), mcTmp));
2189 guessPrecision *= 2;
2190 } while (guessPrecision < targetPrecision + 2);
2191
2192 BigDecimal result;
2193 RoundingMode targetRm = mc.getRoundingMode();
2194 if (targetRm == RoundingMode.UNNECESSARY || originalPrecision == 0) {
2195 RoundingMode tmpRm =
2196 (targetRm == RoundingMode.UNNECESSARY) ? RoundingMode.DOWN : targetRm;
2197 MathContext mcTmp = new MathContext(targetPrecision, tmpRm);
2198 result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mcTmp);
2199
2200 // If result*result != this numerically, the square
2201 // root isn't exact
2202 if (this.subtract(result.multiply(result)).compareTo(ZERO) != 0) {
2203 throw new ArithmeticException("Computed square root not exact.");
2204 }
2205 } else {
2206 result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mc);
2207 }
2208
2209 if (result.scale() != preferredScale) {
2210 // The preferred scale of an add is
2211 // max(addend.scale(), augend.scale()). Therefore, if
2212 // the scale of the result is first minimized using
2213 // stripTrailingZeros(), adding a zero of the
2214 // preferred scale rounding the correct precision will
2215 // perform the proper scale vs precision tradeoffs.
2216 result = result.stripTrailingZeros().
2217 add(zeroWithFinalPreferredScale,
2218 new MathContext(originalPrecision, RoundingMode.UNNECESSARY));
2219 }
2220 assert squareRootResultAssertions(result, mc);
2221 return result;
2222 } else {
2223 switch (signum) {
2224 case -1:
2225 throw new ArithmeticException("Attempted square root " +
2226 "of negative BigDecimal");
2227 case 0:
2228 return valueOf(0L, scale()/2);
2229
2230 default:
2231 throw new AssertionError("Bad value from signum");
2232 }
2233 }
2234 }
2235
2236 private boolean isPowerOfTen() {
2237 return BigInteger.ONE.equals(this.unscaledValue());
2238 }
2239
2240 /**
2241 * For nonzero values, check numerical correctness properties of
2242 * the computed result for the chosen rounding mode.
2243 *
2244 * For the directed roundings, for DOWN and FLOOR, result^2 must
2245 * be {@code <=} the input and (result+ulp)^2 must be {@code >} the
2246 * input. Conversely, for UP and CEIL, result^2 must be {@code >=} the
2247 * input and (result-ulp)^2 must be {@code <} the input.
2248 */
2249 private boolean squareRootResultAssertions(BigDecimal result, MathContext mc) {
2250 if (result.signum() == 0) {
2251 return squareRootZeroResultAssertions(result, mc);
2252 } else {
2253 RoundingMode rm = mc.getRoundingMode();
2254 BigDecimal ulp = result.ulp();
2255 BigDecimal neighborUp = result.add(ulp);
2256 // Make neighbor down accurate even for powers of ten
2257 if (this.isPowerOfTen()) {
2258 ulp = ulp.divide(TEN);
2259 }
2260 BigDecimal neighborDown = result.subtract(ulp);
2261
2262 // Both the starting value and result should be nonzero and positive.
2263 if (result.signum() != 1 ||
2264 this.signum() != 1) {
2265 return false;
2266 }
2267
2268 switch (rm) {
2269 case DOWN:
2270 case FLOOR:
2271 return
2272 result.multiply(result).compareTo(this) <= 0 &&
2273 neighborUp.multiply(neighborUp).compareTo(this) > 0;
2274
2275 case UP:
2276 case CEILING:
2277 return
2278 result.multiply(result).compareTo(this) >= 0 &&
2279 neighborDown.multiply(neighborDown).compareTo(this) < 0;
2280
2281 case HALF_DOWN:
2282 case HALF_EVEN:
2283 case HALF_UP:
2284 BigDecimal err = result.multiply(result).subtract(this).abs();
2285 BigDecimal errUp = neighborUp.multiply(neighborUp).subtract(this);
2286 BigDecimal errDown = this.subtract(neighborDown.multiply(neighborDown));
2287 // All error values should be positive so don't need to
2288 // compare absolute values.
2289
2290 int err_comp_errUp = err.compareTo(errUp);
2291 int err_comp_errDown = err.compareTo(errDown);
2292
2293 return
2294 errUp.signum() == 1 &&
2295 errDown.signum() == 1 &&
2296
2297 err_comp_errUp <= 0 &&
2298 err_comp_errDown <= 0 &&
2299
2300 ((err_comp_errUp == 0 ) ? err_comp_errDown < 0 : true) &&
2301 ((err_comp_errDown == 0 ) ? err_comp_errUp < 0 : true);
2302 // && could check for digit conditions for ties too
2303
2304 default: // Definition of UNNECESSARY already verified.
2305 return true;
2306 }
2307 }
2308 }
2309
2310 private boolean squareRootZeroResultAssertions(BigDecimal result, MathContext mc) {
2311 return this.compareTo(ZERO) == 0;
2312 }
2313
2314 /**
2315 * Returns a {@code BigDecimal} whose value is
2316 * <code>(this<sup>n</sup>)</code>, The power is computed exactly, to
2317 * unlimited precision.
2318 *
2319 * <p>The parameter {@code n} must be in the range 0 through
2320 * 999999999, inclusive. {@code ZERO.pow(0)} returns {@link
2321 * #ONE}.
2322 *
|
2187 // approx = 0.5 * (approx + fraction / approx)
2188 approx = ONE_HALF.multiply(approx.add(working.divide(approx, mcTmp), mcTmp));
2189 guessPrecision *= 2;
2190 } while (guessPrecision < targetPrecision + 2);
2191
2192 BigDecimal result;
2193 RoundingMode targetRm = mc.getRoundingMode();
2194 if (targetRm == RoundingMode.UNNECESSARY || originalPrecision == 0) {
2195 RoundingMode tmpRm =
2196 (targetRm == RoundingMode.UNNECESSARY) ? RoundingMode.DOWN : targetRm;
2197 MathContext mcTmp = new MathContext(targetPrecision, tmpRm);
2198 result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mcTmp);
2199
2200 // If result*result != this numerically, the square
2201 // root isn't exact
2202 if (this.subtract(result.multiply(result)).compareTo(ZERO) != 0) {
2203 throw new ArithmeticException("Computed square root not exact.");
2204 }
2205 } else {
2206 result = approx.scaleByPowerOfTen(-scaleAdjust/2).round(mc);
2207
2208 switch (targetRm) {
2209 case DOWN:
2210 case FLOOR:
2211 // Check if too big
2212 if (result.multiply(result).compareTo(this) > 0 ) {
2213 BigDecimal ulp = ulp = result.ulp();
2214 // Adjust increment down in case of 1.0 = 10^0
2215 // since the next smaller number is only 1/10
2216 // as far way as the next larger at exponent
2217 // boundaries. Test approx and *not* result to
2218 // avoid having to detect an arbitrary power of ten.
2219 if (approx.compareTo(ONE) == 0) {
2220 ulp = ulp.multiply(ONE_TENTH);
2221 }
2222 result = result.subtract(ulp);
2223 }
2224 break;
2225
2226 case UP:
2227 case CEILING:
2228 // Check if too small
2229 if (result.multiply(result).compareTo(this) < 0 ) {
2230 result = result.add(result.ulp());
2231 }
2232 break;
2233
2234 default:
2235 // Do nothing for half-way cases
2236 // HALF_DOWN, HALF_EVEN, HALF_UP
2237 // See fix-up in paper, up down by *1/2* ulp
2238
2239 break;
2240 }
2241
2242 }
2243
2244 // Test numerical properties at full precision before any
2245 // scale adjustments.
2246 assert squareRootResultAssertions(result, mc);
2247 if (result.scale() != preferredScale) {
2248 // The preferred scale of an add is
2249 // max(addend.scale(), augend.scale()). Therefore, if
2250 // the scale of the result is first minimized using
2251 // stripTrailingZeros(), adding a zero of the
2252 // preferred scale rounding the correct precision will
2253 // perform the proper scale vs precision tradeoffs.
2254 result = result.stripTrailingZeros().
2255 add(zeroWithFinalPreferredScale,
2256 new MathContext(originalPrecision, RoundingMode.UNNECESSARY));
2257 }
2258 return result;
2259 } else {
2260 switch (signum) {
2261 case -1:
2262 throw new ArithmeticException("Attempted square root " +
2263 "of negative BigDecimal");
2264 case 0:
2265 return valueOf(0L, scale()/2);
2266
2267 default:
2268 throw new AssertionError("Bad value from signum");
2269 }
2270 }
2271 }
2272
2273 private boolean isPowerOfTen() {
2274 return BigInteger.ONE.equals(this.unscaledValue());
2275 }
2276
2277 /**
2278 * For nonzero values, check numerical correctness properties of
2279 * the computed result for the chosen rounding mode.
2280 *
2281 * For the directed roundings, for DOWN and FLOOR, result^2 must
2282 * be {@code <=} the input and (result+ulp)^2 must be {@code >} the
2283 * input. Conversely, for UP and CEIL, result^2 must be {@code >=} the
2284 * input and (result-ulp)^2 must be {@code <} the input.
2285 */
2286 private boolean squareRootResultAssertions(BigDecimal result, MathContext mc) {
2287 if (result.signum() == 0) {
2288 return squareRootZeroResultAssertions(result, mc);
2289 } else {
2290 RoundingMode rm = mc.getRoundingMode();
2291 BigDecimal ulp = result.ulp();
2292 BigDecimal neighborUp = result.add(ulp);
2293 // Make neighbor down accurate even for powers of ten
2294 if (result.isPowerOfTen()) {
2295 ulp = ulp.divide(TEN);
2296 }
2297 BigDecimal neighborDown = result.subtract(ulp);
2298
2299 // Both the starting value and result should be nonzero and positive.
2300 if (result.signum() != 1 ||
2301 this.signum() != 1) {
2302 return false;
2303 }
2304
2305 switch (rm) {
2306 case DOWN:
2307 case FLOOR:
2308 assert
2309 result.multiply(result).compareTo(this) <= 0 &&
2310 neighborUp.multiply(neighborUp).compareTo(this) > 0:
2311 "Square of result out for bounds rounding " + rm;
2312 return true;
2313
2314 case UP:
2315 case CEILING:
2316 assert
2317 result.multiply(result).compareTo(this) >= 0 &&
2318 neighborDown.multiply(neighborDown).compareTo(this) < 0:
2319 "Square of result out for bounds rounding " + rm;
2320 return true;
2321
2322
2323 case HALF_DOWN:
2324 case HALF_EVEN:
2325 case HALF_UP:
2326 BigDecimal err = result.multiply(result).subtract(this).abs();
2327 BigDecimal errUp = neighborUp.multiply(neighborUp).subtract(this);
2328 BigDecimal errDown = this.subtract(neighborDown.multiply(neighborDown));
2329 // All error values should be positive so don't need to
2330 // compare absolute values.
2331
2332 int err_comp_errUp = err.compareTo(errUp);
2333 int err_comp_errDown = err.compareTo(errDown);
2334
2335 assert
2336 errUp.signum() == 1 &&
2337 errDown.signum() == 1 :
2338 "Errors of neighbors squared don't have correct signs";
2339
2340 assert
2341 err_comp_errUp <= 0 &&
2342 err_comp_errDown <= 0 :
2343 "Computed square root has larger error than neighbors for " + rm;
2344
2345 assert
2346 ((err_comp_errUp == 0 ) ? err_comp_errDown < 0 : true) &&
2347 ((err_comp_errDown == 0 ) ? err_comp_errUp < 0 : true) :
2348 "Incorrect error relationships";
2349 // && could check for digit conditions for ties too
2350 return true;
2351
2352 default: // Definition of UNNECESSARY already verified.
2353 return true;
2354 }
2355 }
2356 }
2357
2358 private boolean squareRootZeroResultAssertions(BigDecimal result, MathContext mc) {
2359 return this.compareTo(ZERO) == 0;
2360 }
2361
2362 /**
2363 * Returns a {@code BigDecimal} whose value is
2364 * <code>(this<sup>n</sup>)</code>, The power is computed exactly, to
2365 * unlimited precision.
2366 *
2367 * <p>The parameter {@code n} must be in the range 0 through
2368 * 999999999, inclusive. {@code ZERO.pow(0)} returns {@link
2369 * #ONE}.
2370 *
|