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

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

@@ -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,78 @@
         // 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.
+     *
+     * @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() {
+        // Obtain number of bits needed to represent the number.
+        int bitLength = (int)this.bitLength();
+        if (bitLength != this.bitLength()) {
+            throw new ArithmeticException("bitLength() integer overflow");
+        }
+
+        // Initial value length is one more than floor(bitLength()/2).
+        int initBits = bitLength / 2 + 1;
+
+        // The bit length of intermediate values will be no more than 2 bits
+        // larger than the initial bit length, so size the values accordingly
+        // to minimize allocations. The number of ints in the initial guess is
+        // (initBits + 2 + 31)/32.
+        int intLength = (initBits + 33) / 32;
+
+        // Fill the lower initBits bits of the value array with ones.  This
+        // should provide an initial estimate greater than or equal to the
+        // actual integer square root.
+        int[] val = new int[intLength];
+        int initInts = initBits / 32;
+        Arrays.fill(val, 0, initInts, -1);
+        int excessBits = initBits % 32;
+        if (excessBits != 0) {
+            val[initInts] =
+                (int)((0xffffffff & LONG_MASK) >>> (32 - excessBits));
+        }
+
+        // Create the initial estimate from the value array.
+        MutableBigInteger xk = new MutableBigInteger(val);
+
+        // Create intermediate and next value containers and the constant two.
+        MutableBigInteger tmp = new MutableBigInteger(new int[intLength]);
+        MutableBigInteger xk1 = new MutableBigInteger(new int[intLength]);
+        MutableBigInteger two = new MutableBigInteger(2);
+
+        do {
+            // xk1 = (xk + n/xk)/2
+            this.divide(xk, tmp, false);
+            tmp.add(xk);
+            tmp.divide(two, xk1, false);
+
+            if (xk1.compare(xk) >= 0) {
+                return xk;
+            }
+
+            // xk = xk1
+            xk.copyValue(xk1);
+
+            xk1.reset();
+            tmp.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.