< prev index next >

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

Print this page

        

@@ -344,10 +344,20 @@
      * @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

@@ -1993,10 +2003,223 @@
         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
+                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
      * <code>(this<sup>n</sup>)</code>, The power is computed exactly, to
      * unlimited precision.
      *
< prev index next >