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 ****
  /*
!  * Copyright (c) 1999, 2013, 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
--- 1,7 ----
  /*
!  * 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,1874 ****
--- 1865,1942 ----
          // 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.