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

rev 13121 : 8032027: Add BigInteger square root methods
Summary: Add sqrt() and sqrtAndReminder() using Newton iteration
Reviewed-by: XXX
rev 13122 : 8032027: Add BigInteger square root methods
Summary: Add sqrt() and sqrtAndReminder() using Newton iteration
Reviewed-by: darcy

@@ -1,7 +1,7 @@
 /*
- * Copyright (c) 1999, 2013, Oracle and/or its affiliates. All rights reserved.
+ * Copyright (c) 1999, 2015, Oracle and/or its affiliates. All rights reserved.
  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
  *
  * This code is free software; you can redistribute it and/or modify it
  * under the terms of the GNU General Public License version 2 only, as
  * published by the Free Software Foundation.  Oracle designates this

@@ -1865,10 +1865,87 @@
         // n - q*dlong == r && 0 <= r <dLong, hence we're done.
         return (r << 32) | (q & LONG_MASK);
     }
 
     /**
+     * Calculate the integer square root {@code floor(sqrt(this))} where
+     * {@code sqrt(.)} denotes the mathematical square root. The contents of
+     * {@code this} are <b>not</b> changed. The value of {@code this} is assumed
+     * to be non-negative.
+     *
+     * @implNote The implementation is based on the material in Henry S. Warren,
+     * Jr., <i>Hacker's Delight (2nd ed.)</i> (Addison Wesley, 2013), 279-282.
+     *
+     * @throws ArithmeticException if the value returned by {@code bitLength()}
+     * overflows the range of {@code int}.
+     * @return the integer square root of {@code this}
+     * @since 1.9
+     */
+    MutableBigInteger sqrt() {
+        // Special cases.
+        if (this.isZero()) {
+            return new MutableBigInteger(0);
+        } else if (this.value.length == 1
+                && (this.value[0] & LONG_MASK) < 4) { // result is unity
+            return ONE;
+        }
+
+        // Set up the initial estimate of the iteration.
+        MutableBigInteger xk;
+        if (bitLength() < 63) {
+            // Use the square root of the positive long value.
+            long v = new BigInteger(this.value, 1).longValueExact();
+            double s = Math.sqrt(v);
+            BigInteger bi = BigInteger.valueOf((long)Math.floor(s));
+            xk = new MutableBigInteger(bi.mag);
+        } else {
+            // Obtain the bitLength > 63.
+            int bitLength = (int) this.bitLength();
+            if (bitLength != this.bitLength()) {
+                throw new ArithmeticException("bitLength() integer overflow");
+            }
+
+            // Determine an even valued right shift into positive long range.
+            int shift = bitLength - 63;
+            if (shift % 2 == 1) {
+                shift++;
+            }
+
+            // Shift the value into positive long range.
+            xk = new MutableBigInteger(this);
+            xk.rightShift(shift);
+            xk.normalize();
+
+            // Use the square root of the shifted value as an approximation.
+            double d = new BigInteger(xk.value, 1).doubleValue();
+            BigInteger bi = BigInteger.valueOf((long)Math.floor(Math.sqrt(d)));
+            xk = new MutableBigInteger(bi.mag);
+
+            // Shift the approximate square root back into the original range.
+            xk.leftShift(shift / 2);
+        }
+
+        // Refine the estimate.
+        MutableBigInteger xk1 = new MutableBigInteger();
+        do {
+            // xk1 = (xk + n/xk)/2
+            this.divide(xk, xk1, false);
+            xk1.add(xk);
+            xk1.rightShift(1);
+
+            if (xk1.compare(xk) >= 0) {
+                return xk;
+            }
+
+            // xk = xk1
+            xk.copyValue(xk1);
+
+            xk1.reset();
+        } while (true);
+    }
+
+    /**
      * Calculate GCD of this and b. This and b are changed by the computation.
      */
     MutableBigInteger hybridGCD(MutableBigInteger b) {
         // Use Euclid's algorithm until the numbers are approximately the
         // same length, then use the binary GCD algorithm to find the GCD.