1 /*
   2  * Copyright (c) 1996, 2013, Oracle and/or its affiliates. All rights reserved.
   3  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
   4  *
   5  * This code is free software; you can redistribute it and/or modify it
   6  * under the terms of the GNU General Public License version 2 only, as
   7  * published by the Free Software Foundation.  Oracle designates this
   8  * particular file as subject to the "Classpath" exception as provided
   9  * by Oracle in the LICENSE file that accompanied this code.
  10  *
  11  * This code is distributed in the hope that it will be useful, but WITHOUT
  12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14  * version 2 for more details (a copy is included in the LICENSE file that
  15  * accompanied this code).
  16  *
  17  * You should have received a copy of the GNU General Public License version
  18  * 2 along with this work; if not, write to the Free Software Foundation,
  19  * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
  20  *
  21  * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
  22  * or visit www.oracle.com if you need additional information or have any
  23  * questions.
  24  */
  25 
  26 /*
  27  * Portions Copyright (c) 1995  Colin Plumb.  All rights reserved.
  28  */
  29 
  30 package java.math;
  31 
  32 import java.io.IOException;
  33 import java.io.ObjectInputStream;
  34 import java.io.ObjectOutputStream;
  35 import java.io.ObjectStreamField;
  36 import java.util.Arrays;
  37 import java.util.Random;
  38 import sun.misc.DoubleConsts;
  39 import sun.misc.FloatConsts;
  40 
  41 /**
  42  * Immutable arbitrary-precision integers.  All operations behave as if
  43  * BigIntegers were represented in two's-complement notation (like Java's
  44  * primitive integer types).  BigInteger provides analogues to all of Java's
  45  * primitive integer operators, and all relevant methods from java.lang.Math.
  46  * Additionally, BigInteger provides operations for modular arithmetic, GCD
  47  * calculation, primality testing, prime generation, bit manipulation,
  48  * and a few other miscellaneous operations.
  49  *
  50  * <p>Semantics of arithmetic operations exactly mimic those of Java's integer
  51  * arithmetic operators, as defined in <i>The Java Language Specification</i>.
  52  * For example, division by zero throws an {@code ArithmeticException}, and
  53  * division of a negative by a positive yields a negative (or zero) remainder.
  54  * All of the details in the Spec concerning overflow are ignored, as
  55  * BigIntegers are made as large as necessary to accommodate the results of an
  56  * operation.
  57  *
  58  * <p>Semantics of shift operations extend those of Java's shift operators
  59  * to allow for negative shift distances.  A right-shift with a negative
  60  * shift distance results in a left shift, and vice-versa.  The unsigned
  61  * right shift operator ({@code >>>}) is omitted, as this operation makes
  62  * little sense in combination with the "infinite word size" abstraction
  63  * provided by this class.
  64  *
  65  * <p>Semantics of bitwise logical operations exactly mimic those of Java's
  66  * bitwise integer operators.  The binary operators ({@code and},
  67  * {@code or}, {@code xor}) implicitly perform sign extension on the shorter
  68  * of the two operands prior to performing the operation.
  69  *
  70  * <p>Comparison operations perform signed integer comparisons, analogous to
  71  * those performed by Java's relational and equality operators.
  72  *
  73  * <p>Modular arithmetic operations are provided to compute residues, perform
  74  * exponentiation, and compute multiplicative inverses.  These methods always
  75  * return a non-negative result, between {@code 0} and {@code (modulus - 1)},
  76  * inclusive.
  77  *
  78  * <p>Bit operations operate on a single bit of the two's-complement
  79  * representation of their operand.  If necessary, the operand is sign-
  80  * extended so that it contains the designated bit.  None of the single-bit
  81  * operations can produce a BigInteger with a different sign from the
  82  * BigInteger being operated on, as they affect only a single bit, and the
  83  * "infinite word size" abstraction provided by this class ensures that there
  84  * are infinitely many "virtual sign bits" preceding each BigInteger.
  85  *
  86  * <p>For the sake of brevity and clarity, pseudo-code is used throughout the
  87  * descriptions of BigInteger methods.  The pseudo-code expression
  88  * {@code (i + j)} is shorthand for "a BigInteger whose value is
  89  * that of the BigInteger {@code i} plus that of the BigInteger {@code j}."
  90  * The pseudo-code expression {@code (i == j)} is shorthand for
  91  * "{@code true} if and only if the BigInteger {@code i} represents the same
  92  * value as the BigInteger {@code j}."  Other pseudo-code expressions are
  93  * interpreted similarly.
  94  *
  95  * <p>All methods and constructors in this class throw
  96  * {@code NullPointerException} when passed
  97  * a null object reference for any input parameter.
  98  *
  99  * @see     BigDecimal
 100  * @author  Josh Bloch
 101  * @author  Michael McCloskey
 102  * @author  Alan Eliasen
 103  * @author  Timothy Buktu
 104  * @since JDK1.1
 105  */
 106 
 107 public class BigInteger extends Number implements Comparable<BigInteger> {
 108     /**
 109      * The signum of this BigInteger: -1 for negative, 0 for zero, or
 110      * 1 for positive.  Note that the BigInteger zero <i>must</i> have
 111      * a signum of 0.  This is necessary to ensures that there is exactly one
 112      * representation for each BigInteger value.
 113      *
 114      * @serial
 115      */
 116     final int signum;
 117 
 118     /**
 119      * The magnitude of this BigInteger, in <i>big-endian</i> order: the
 120      * zeroth element of this array is the most-significant int of the
 121      * magnitude.  The magnitude must be "minimal" in that the most-significant
 122      * int ({@code mag[0]}) must be non-zero.  This is necessary to
 123      * ensure that there is exactly one representation for each BigInteger
 124      * value.  Note that this implies that the BigInteger zero has a
 125      * zero-length mag array.
 126      */
 127     final int[] mag;
 128 
 129     // These "redundant fields" are initialized with recognizable nonsense
 130     // values, and cached the first time they are needed (or never, if they
 131     // aren't needed).
 132 
 133      /**
 134      * One plus the bitCount of this BigInteger. Zeros means unitialized.
 135      *
 136      * @serial
 137      * @see #bitCount
 138      * @deprecated Deprecated since logical value is offset from stored
 139      * value and correction factor is applied in accessor method.
 140      */
 141     @Deprecated
 142     private int bitCount;
 143 
 144     /**
 145      * One plus the bitLength of this BigInteger. Zeros means unitialized.
 146      * (either value is acceptable).
 147      *
 148      * @serial
 149      * @see #bitLength()
 150      * @deprecated Deprecated since logical value is offset from stored
 151      * value and correction factor is applied in accessor method.
 152      */
 153     @Deprecated
 154     private int bitLength;
 155 
 156     /**
 157      * Two plus the lowest set bit of this BigInteger, as returned by
 158      * getLowestSetBit().
 159      *
 160      * @serial
 161      * @see #getLowestSetBit
 162      * @deprecated Deprecated since logical value is offset from stored
 163      * value and correction factor is applied in accessor method.
 164      */
 165     @Deprecated
 166     private int lowestSetBit;
 167 
 168     /**
 169      * Two plus the index of the lowest-order int in the magnitude of this
 170      * BigInteger that contains a nonzero int, or -2 (either value is acceptable).
 171      * The least significant int has int-number 0, the next int in order of
 172      * increasing significance has int-number 1, and so forth.
 173      * @deprecated Deprecated since logical value is offset from stored
 174      * value and correction factor is applied in accessor method.
 175      */
 176     @Deprecated
 177     private int firstNonzeroIntNum;
 178 
 179     /**
 180      * This mask is used to obtain the value of an int as if it were unsigned.
 181      */
 182     final static long LONG_MASK = 0xffffffffL;
 183 
 184     /**
 185      * The threshold value for using Karatsuba multiplication.  If the number
 186      * of ints in both mag arrays are greater than this number, then
 187      * Karatsuba multiplication will be used.   This value is found
 188      * experimentally to work well.
 189      */
 190     private static final int KARATSUBA_THRESHOLD = 50;
 191 
 192     /**
 193      * The threshold value for using 3-way Toom-Cook multiplication.
 194      * If the number of ints in each mag array is greater than the
 195      * Karatsuba threshold, and the number of ints in at least one of
 196      * the mag arrays is greater than this threshold, then Toom-Cook
 197      * multiplication will be used.
 198      */
 199     private static final int TOOM_COOK_THRESHOLD = 75;
 200 
 201     /**
 202      * The threshold value for using Karatsuba squaring.  If the number
 203      * of ints in the number are larger than this value,
 204      * Karatsuba squaring will be used.   This value is found
 205      * experimentally to work well.
 206      */
 207     private static final int KARATSUBA_SQUARE_THRESHOLD = 90;
 208 
 209     /**
 210      * The threshold value for using Toom-Cook squaring.  If the number
 211      * of ints in the number are larger than this value,
 212      * Toom-Cook squaring will be used.   This value is found
 213      * experimentally to work well.
 214      */
 215     private static final int TOOM_COOK_SQUARE_THRESHOLD = 140;
 216 
 217     /**
 218      * The threshold value for using Burnikel-Ziegler division.  If the number
 219      * of ints in the number are larger than this value,
 220      * Burnikel-Ziegler division will be used.   This value is found
 221      * experimentally to work well.
 222      */
 223     static final int BURNIKEL_ZIEGLER_THRESHOLD = 50;
 224 
 225     /**
 226      * The threshold value for using Schoenhage recursive base conversion. If
 227      * the number of ints in the number are larger than this value,
 228      * the Schoenhage algorithm will be used.  In practice, it appears that the
 229      * Schoenhage routine is faster for any threshold down to 2, and is
 230      * relatively flat for thresholds between 2-25, so this choice may be
 231      * varied within this range for very small effect.
 232      */
 233     private static final int SCHOENHAGE_BASE_CONVERSION_THRESHOLD = 8;
 234 
 235     //Constructors
 236 
 237     /**
 238      * Translates a byte array containing the two's-complement binary
 239      * representation of a BigInteger into a BigInteger.  The input array is
 240      * assumed to be in <i>big-endian</i> byte-order: the most significant
 241      * byte is in the zeroth element.
 242      *
 243      * @param  val big-endian two's-complement binary representation of
 244      *         BigInteger.
 245      * @throws NumberFormatException {@code val} is zero bytes long.
 246      */
 247     public BigInteger(byte[] val) {
 248         if (val.length == 0)
 249             throw new NumberFormatException("Zero length BigInteger");
 250 
 251         if (val[0] < 0) {
 252             mag = makePositive(val);
 253             signum = -1;
 254         } else {
 255             mag = stripLeadingZeroBytes(val);
 256             signum = (mag.length == 0 ? 0 : 1);
 257         }
 258     }
 259 
 260     /**
 261      * This private constructor translates an int array containing the
 262      * two's-complement binary representation of a BigInteger into a
 263      * BigInteger. The input array is assumed to be in <i>big-endian</i>
 264      * int-order: the most significant int is in the zeroth element.
 265      */
 266     private BigInteger(int[] val) {
 267         if (val.length == 0)
 268             throw new NumberFormatException("Zero length BigInteger");
 269 
 270         if (val[0] < 0) {
 271             mag = makePositive(val);
 272             signum = -1;
 273         } else {
 274             mag = trustedStripLeadingZeroInts(val);
 275             signum = (mag.length == 0 ? 0 : 1);
 276         }
 277     }
 278 
 279     /**
 280      * Translates the sign-magnitude representation of a BigInteger into a
 281      * BigInteger.  The sign is represented as an integer signum value: -1 for
 282      * negative, 0 for zero, or 1 for positive.  The magnitude is a byte array
 283      * in <i>big-endian</i> byte-order: the most significant byte is in the
 284      * zeroth element.  A zero-length magnitude array is permissible, and will
 285      * result in a BigInteger value of 0, whether signum is -1, 0 or 1.
 286      *
 287      * @param  signum signum of the number (-1 for negative, 0 for zero, 1
 288      *         for positive).
 289      * @param  magnitude big-endian binary representation of the magnitude of
 290      *         the number.
 291      * @throws NumberFormatException {@code signum} is not one of the three
 292      *         legal values (-1, 0, and 1), or {@code signum} is 0 and
 293      *         {@code magnitude} contains one or more non-zero bytes.
 294      */
 295     public BigInteger(int signum, byte[] magnitude) {
 296         this.mag = stripLeadingZeroBytes(magnitude);
 297 
 298         if (signum < -1 || signum > 1)
 299             throw(new NumberFormatException("Invalid signum value"));
 300 
 301         if (this.mag.length==0) {
 302             this.signum = 0;
 303         } else {
 304             if (signum == 0)
 305                 throw(new NumberFormatException("signum-magnitude mismatch"));
 306             this.signum = signum;
 307         }
 308     }
 309 
 310     /**
 311      * A constructor for internal use that translates the sign-magnitude
 312      * representation of a BigInteger into a BigInteger. It checks the
 313      * arguments and copies the magnitude so this constructor would be
 314      * safe for external use.
 315      */
 316     private BigInteger(int signum, int[] magnitude) {
 317         this.mag = stripLeadingZeroInts(magnitude);
 318 
 319         if (signum < -1 || signum > 1)
 320             throw(new NumberFormatException("Invalid signum value"));
 321 
 322         if (this.mag.length==0) {
 323             this.signum = 0;
 324         } else {
 325             if (signum == 0)
 326                 throw(new NumberFormatException("signum-magnitude mismatch"));
 327             this.signum = signum;
 328         }
 329     }
 330 
 331     /**
 332      * Translates the String representation of a BigInteger in the
 333      * specified radix into a BigInteger.  The String representation
 334      * consists of an optional minus or plus sign followed by a
 335      * sequence of one or more digits in the specified radix.  The
 336      * character-to-digit mapping is provided by {@code
 337      * Character.digit}.  The String may not contain any extraneous
 338      * characters (whitespace, for example).
 339      *
 340      * @param val String representation of BigInteger.
 341      * @param radix radix to be used in interpreting {@code val}.
 342      * @throws NumberFormatException {@code val} is not a valid representation
 343      *         of a BigInteger in the specified radix, or {@code radix} is
 344      *         outside the range from {@link Character#MIN_RADIX} to
 345      *         {@link Character#MAX_RADIX}, inclusive.
 346      * @see    Character#digit
 347      */
 348     public BigInteger(String val, int radix) {
 349         int cursor = 0, numDigits;
 350         final int len = val.length();
 351 
 352         if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
 353             throw new NumberFormatException("Radix out of range");
 354         if (len == 0)
 355             throw new NumberFormatException("Zero length BigInteger");
 356 
 357         // Check for at most one leading sign
 358         int sign = 1;
 359         int index1 = val.lastIndexOf('-');
 360         int index2 = val.lastIndexOf('+');
 361         if ((index1 + index2) <= -1) {
 362             // No leading sign character or at most one leading sign character
 363             if (index1 == 0 || index2 == 0) {
 364                 cursor = 1;
 365                 if (len == 1)
 366                     throw new NumberFormatException("Zero length BigInteger");
 367             }
 368             if (index1 == 0)
 369                 sign = -1;
 370         } else
 371             throw new NumberFormatException("Illegal embedded sign character");
 372 
 373         // Skip leading zeros and compute number of digits in magnitude
 374         while (cursor < len &&
 375                Character.digit(val.charAt(cursor), radix) == 0)
 376             cursor++;
 377         if (cursor == len) {
 378             signum = 0;
 379             mag = ZERO.mag;
 380             return;
 381         }
 382 
 383         numDigits = len - cursor;
 384         signum = sign;
 385 
 386         // Pre-allocate array of expected size. May be too large but can
 387         // never be too small. Typically exact.
 388         int numBits = (int)(((numDigits * bitsPerDigit[radix]) >>> 10) + 1);
 389         int numWords = (numBits + 31) >>> 5;
 390         int[] magnitude = new int[numWords];
 391 
 392         // Process first (potentially short) digit group
 393         int firstGroupLen = numDigits % digitsPerInt[radix];
 394         if (firstGroupLen == 0)
 395             firstGroupLen = digitsPerInt[radix];
 396         String group = val.substring(cursor, cursor += firstGroupLen);
 397         magnitude[numWords - 1] = Integer.parseInt(group, radix);
 398         if (magnitude[numWords - 1] < 0)
 399             throw new NumberFormatException("Illegal digit");
 400 
 401         // Process remaining digit groups
 402         int superRadix = intRadix[radix];
 403         int groupVal = 0;
 404         while (cursor < len) {
 405             group = val.substring(cursor, cursor += digitsPerInt[radix]);
 406             groupVal = Integer.parseInt(group, radix);
 407             if (groupVal < 0)
 408                 throw new NumberFormatException("Illegal digit");
 409             destructiveMulAdd(magnitude, superRadix, groupVal);
 410         }
 411         // Required for cases where the array was overallocated.
 412         mag = trustedStripLeadingZeroInts(magnitude);
 413     }
 414 
 415     /*
 416      * Constructs a new BigInteger using a char array with radix=10.
 417      * Sign is precalculated outside and not allowed in the val.
 418      */
 419     BigInteger(char[] val, int sign, int len) {
 420         int cursor = 0, numDigits;
 421 
 422         // Skip leading zeros and compute number of digits in magnitude
 423         while (cursor < len && Character.digit(val[cursor], 10) == 0) {
 424             cursor++;
 425         }
 426         if (cursor == len) {
 427             signum = 0;
 428             mag = ZERO.mag;
 429             return;
 430         }
 431 
 432         numDigits = len - cursor;
 433         signum = sign;
 434         // Pre-allocate array of expected size
 435         int numWords;
 436         if (len < 10) {
 437             numWords = 1;
 438         } else {
 439             int numBits = (int)(((numDigits * bitsPerDigit[10]) >>> 10) + 1);
 440             numWords = (numBits + 31) >>> 5;
 441         }
 442         int[] magnitude = new int[numWords];
 443 
 444         // Process first (potentially short) digit group
 445         int firstGroupLen = numDigits % digitsPerInt[10];
 446         if (firstGroupLen == 0)
 447             firstGroupLen = digitsPerInt[10];
 448         magnitude[numWords - 1] = parseInt(val, cursor,  cursor += firstGroupLen);
 449 
 450         // Process remaining digit groups
 451         while (cursor < len) {
 452             int groupVal = parseInt(val, cursor, cursor += digitsPerInt[10]);
 453             destructiveMulAdd(magnitude, intRadix[10], groupVal);
 454         }
 455         mag = trustedStripLeadingZeroInts(magnitude);
 456     }
 457 
 458     // Create an integer with the digits between the two indexes
 459     // Assumes start < end. The result may be negative, but it
 460     // is to be treated as an unsigned value.
 461     private int parseInt(char[] source, int start, int end) {
 462         int result = Character.digit(source[start++], 10);
 463         if (result == -1)
 464             throw new NumberFormatException(new String(source));
 465 
 466         for (int index = start; index<end; index++) {
 467             int nextVal = Character.digit(source[index], 10);
 468             if (nextVal == -1)
 469                 throw new NumberFormatException(new String(source));
 470             result = 10*result + nextVal;
 471         }
 472 
 473         return result;
 474     }
 475 
 476     // bitsPerDigit in the given radix times 1024
 477     // Rounded up to avoid underallocation.
 478     private static long bitsPerDigit[] = { 0, 0,
 479         1024, 1624, 2048, 2378, 2648, 2875, 3072, 3247, 3402, 3543, 3672,
 480         3790, 3899, 4001, 4096, 4186, 4271, 4350, 4426, 4498, 4567, 4633,
 481         4696, 4756, 4814, 4870, 4923, 4975, 5025, 5074, 5120, 5166, 5210,
 482                                            5253, 5295};
 483 
 484     // Multiply x array times word y in place, and add word z
 485     private static void destructiveMulAdd(int[] x, int y, int z) {
 486         // Perform the multiplication word by word
 487         long ylong = y & LONG_MASK;
 488         long zlong = z & LONG_MASK;
 489         int len = x.length;
 490 
 491         long product = 0;
 492         long carry = 0;
 493         for (int i = len-1; i >= 0; i--) {
 494             product = ylong * (x[i] & LONG_MASK) + carry;
 495             x[i] = (int)product;
 496             carry = product >>> 32;
 497         }
 498 
 499         // Perform the addition
 500         long sum = (x[len-1] & LONG_MASK) + zlong;
 501         x[len-1] = (int)sum;
 502         carry = sum >>> 32;
 503         for (int i = len-2; i >= 0; i--) {
 504             sum = (x[i] & LONG_MASK) + carry;
 505             x[i] = (int)sum;
 506             carry = sum >>> 32;
 507         }
 508     }
 509 
 510     /**
 511      * Translates the decimal String representation of a BigInteger into a
 512      * BigInteger.  The String representation consists of an optional minus
 513      * sign followed by a sequence of one or more decimal digits.  The
 514      * character-to-digit mapping is provided by {@code Character.digit}.
 515      * The String may not contain any extraneous characters (whitespace, for
 516      * example).
 517      *
 518      * @param val decimal String representation of BigInteger.
 519      * @throws NumberFormatException {@code val} is not a valid representation
 520      *         of a BigInteger.
 521      * @see    Character#digit
 522      */
 523     public BigInteger(String val) {
 524         this(val, 10);
 525     }
 526 
 527     /**
 528      * Constructs a randomly generated BigInteger, uniformly distributed over
 529      * the range 0 to (2<sup>{@code numBits}</sup> - 1), inclusive.
 530      * The uniformity of the distribution assumes that a fair source of random
 531      * bits is provided in {@code rnd}.  Note that this constructor always
 532      * constructs a non-negative BigInteger.
 533      *
 534      * @param  numBits maximum bitLength of the new BigInteger.
 535      * @param  rnd source of randomness to be used in computing the new
 536      *         BigInteger.
 537      * @throws IllegalArgumentException {@code numBits} is negative.
 538      * @see #bitLength()
 539      */
 540     public BigInteger(int numBits, Random rnd) {
 541         this(1, randomBits(numBits, rnd));
 542     }
 543 
 544     private static byte[] randomBits(int numBits, Random rnd) {
 545         if (numBits < 0)
 546             throw new IllegalArgumentException("numBits must be non-negative");
 547         int numBytes = (int)(((long)numBits+7)/8); // avoid overflow
 548         byte[] randomBits = new byte[numBytes];
 549 
 550         // Generate random bytes and mask out any excess bits
 551         if (numBytes > 0) {
 552             rnd.nextBytes(randomBits);
 553             int excessBits = 8*numBytes - numBits;
 554             randomBits[0] &= (1 << (8-excessBits)) - 1;
 555         }
 556         return randomBits;
 557     }
 558 
 559     /**
 560      * Constructs a randomly generated positive BigInteger that is probably
 561      * prime, with the specified bitLength.
 562      *
 563      * <p>It is recommended that the {@link #probablePrime probablePrime}
 564      * method be used in preference to this constructor unless there
 565      * is a compelling need to specify a certainty.
 566      *
 567      * @param  bitLength bitLength of the returned BigInteger.
 568      * @param  certainty a measure of the uncertainty that the caller is
 569      *         willing to tolerate.  The probability that the new BigInteger
 570      *         represents a prime number will exceed
 571      *         (1 - 1/2<sup>{@code certainty}</sup>).  The execution time of
 572      *         this constructor is proportional to the value of this parameter.
 573      * @param  rnd source of random bits used to select candidates to be
 574      *         tested for primality.
 575      * @throws ArithmeticException {@code bitLength < 2}.
 576      * @see    #bitLength()
 577      */
 578     public BigInteger(int bitLength, int certainty, Random rnd) {
 579         BigInteger prime;
 580 
 581         if (bitLength < 2)
 582             throw new ArithmeticException("bitLength < 2");
 583         prime = (bitLength < SMALL_PRIME_THRESHOLD
 584                                 ? smallPrime(bitLength, certainty, rnd)
 585                                 : largePrime(bitLength, certainty, rnd));
 586         signum = 1;
 587         mag = prime.mag;
 588     }
 589 
 590     // Minimum size in bits that the requested prime number has
 591     // before we use the large prime number generating algorithms.
 592     // The cutoff of 95 was chosen empirically for best performance.
 593     private static final int SMALL_PRIME_THRESHOLD = 95;
 594 
 595     // Certainty required to meet the spec of probablePrime
 596     private static final int DEFAULT_PRIME_CERTAINTY = 100;
 597 
 598     /**
 599      * Returns a positive BigInteger that is probably prime, with the
 600      * specified bitLength. The probability that a BigInteger returned
 601      * by this method is composite does not exceed 2<sup>-100</sup>.
 602      *
 603      * @param  bitLength bitLength of the returned BigInteger.
 604      * @param  rnd source of random bits used to select candidates to be
 605      *         tested for primality.
 606      * @return a BigInteger of {@code bitLength} bits that is probably prime
 607      * @throws ArithmeticException {@code bitLength < 2}.
 608      * @see    #bitLength()
 609      * @since 1.4
 610      */
 611     public static BigInteger probablePrime(int bitLength, Random rnd) {
 612         if (bitLength < 2)
 613             throw new ArithmeticException("bitLength < 2");
 614 
 615         return (bitLength < SMALL_PRIME_THRESHOLD ?
 616                 smallPrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd) :
 617                 largePrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd));
 618     }
 619 
 620     /**
 621      * Find a random number of the specified bitLength that is probably prime.
 622      * This method is used for smaller primes, its performance degrades on
 623      * larger bitlengths.
 624      *
 625      * This method assumes bitLength > 1.
 626      */
 627     private static BigInteger smallPrime(int bitLength, int certainty, Random rnd) {
 628         int magLen = (bitLength + 31) >>> 5;
 629         int temp[] = new int[magLen];
 630         int highBit = 1 << ((bitLength+31) & 0x1f);  // High bit of high int
 631         int highMask = (highBit << 1) - 1;  // Bits to keep in high int
 632 
 633         while(true) {
 634             // Construct a candidate
 635             for (int i=0; i<magLen; i++)
 636                 temp[i] = rnd.nextInt();
 637             temp[0] = (temp[0] & highMask) | highBit;  // Ensure exact length
 638             if (bitLength > 2)
 639                 temp[magLen-1] |= 1;  // Make odd if bitlen > 2
 640 
 641             BigInteger p = new BigInteger(temp, 1);
 642 
 643             // Do cheap "pre-test" if applicable
 644             if (bitLength > 6) {
 645                 long r = p.remainder(SMALL_PRIME_PRODUCT).longValue();
 646                 if ((r%3==0)  || (r%5==0)  || (r%7==0)  || (r%11==0) ||
 647                     (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) ||
 648                     (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0))
 649                     continue; // Candidate is composite; try another
 650             }
 651 
 652             // All candidates of bitLength 2 and 3 are prime by this point
 653             if (bitLength < 4)
 654                 return p;
 655 
 656             // Do expensive test if we survive pre-test (or it's inapplicable)
 657             if (p.primeToCertainty(certainty, rnd))
 658                 return p;
 659         }
 660     }
 661 
 662     private static final BigInteger SMALL_PRIME_PRODUCT
 663                        = valueOf(3L*5*7*11*13*17*19*23*29*31*37*41);
 664 
 665     /**
 666      * Find a random number of the specified bitLength that is probably prime.
 667      * This method is more appropriate for larger bitlengths since it uses
 668      * a sieve to eliminate most composites before using a more expensive
 669      * test.
 670      */
 671     private static BigInteger largePrime(int bitLength, int certainty, Random rnd) {
 672         BigInteger p;
 673         p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
 674         p.mag[p.mag.length-1] &= 0xfffffffe;
 675 
 676         // Use a sieve length likely to contain the next prime number
 677         int searchLen = (bitLength / 20) * 64;
 678         BitSieve searchSieve = new BitSieve(p, searchLen);
 679         BigInteger candidate = searchSieve.retrieve(p, certainty, rnd);
 680 
 681         while ((candidate == null) || (candidate.bitLength() != bitLength)) {
 682             p = p.add(BigInteger.valueOf(2*searchLen));
 683             if (p.bitLength() != bitLength)
 684                 p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
 685             p.mag[p.mag.length-1] &= 0xfffffffe;
 686             searchSieve = new BitSieve(p, searchLen);
 687             candidate = searchSieve.retrieve(p, certainty, rnd);
 688         }
 689         return candidate;
 690     }
 691 
 692    /**
 693     * Returns the first integer greater than this {@code BigInteger} that
 694     * is probably prime.  The probability that the number returned by this
 695     * method is composite does not exceed 2<sup>-100</sup>. This method will
 696     * never skip over a prime when searching: if it returns {@code p}, there
 697     * is no prime {@code q} such that {@code this < q < p}.
 698     *
 699     * @return the first integer greater than this {@code BigInteger} that
 700     *         is probably prime.
 701     * @throws ArithmeticException {@code this < 0}.
 702     * @since 1.5
 703     */
 704     public BigInteger nextProbablePrime() {
 705         if (this.signum < 0)
 706             throw new ArithmeticException("start < 0: " + this);
 707 
 708         // Handle trivial cases
 709         if ((this.signum == 0) || this.equals(ONE))
 710             return TWO;
 711 
 712         BigInteger result = this.add(ONE);
 713 
 714         // Fastpath for small numbers
 715         if (result.bitLength() < SMALL_PRIME_THRESHOLD) {
 716 
 717             // Ensure an odd number
 718             if (!result.testBit(0))
 719                 result = result.add(ONE);
 720 
 721             while(true) {
 722                 // Do cheap "pre-test" if applicable
 723                 if (result.bitLength() > 6) {
 724                     long r = result.remainder(SMALL_PRIME_PRODUCT).longValue();
 725                     if ((r%3==0)  || (r%5==0)  || (r%7==0)  || (r%11==0) ||
 726                         (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) ||
 727                         (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0)) {
 728                         result = result.add(TWO);
 729                         continue; // Candidate is composite; try another
 730                     }
 731                 }
 732 
 733                 // All candidates of bitLength 2 and 3 are prime by this point
 734                 if (result.bitLength() < 4)
 735                     return result;
 736 
 737                 // The expensive test
 738                 if (result.primeToCertainty(DEFAULT_PRIME_CERTAINTY, null))
 739                     return result;
 740 
 741                 result = result.add(TWO);
 742             }
 743         }
 744 
 745         // Start at previous even number
 746         if (result.testBit(0))
 747             result = result.subtract(ONE);
 748 
 749         // Looking for the next large prime
 750         int searchLen = (result.bitLength() / 20) * 64;
 751 
 752         while(true) {
 753            BitSieve searchSieve = new BitSieve(result, searchLen);
 754            BigInteger candidate = searchSieve.retrieve(result,
 755                                                  DEFAULT_PRIME_CERTAINTY, null);
 756            if (candidate != null)
 757                return candidate;
 758            result = result.add(BigInteger.valueOf(2 * searchLen));
 759         }
 760     }
 761 
 762     /**
 763      * Returns {@code true} if this BigInteger is probably prime,
 764      * {@code false} if it's definitely composite.
 765      *
 766      * This method assumes bitLength > 2.
 767      *
 768      * @param  certainty a measure of the uncertainty that the caller is
 769      *         willing to tolerate: if the call returns {@code true}
 770      *         the probability that this BigInteger is prime exceeds
 771      *         {@code (1 - 1/2<sup>certainty</sup>)}.  The execution time of
 772      *         this method is proportional to the value of this parameter.
 773      * @return {@code true} if this BigInteger is probably prime,
 774      *         {@code false} if it's definitely composite.
 775      */
 776     boolean primeToCertainty(int certainty, Random random) {
 777         int rounds = 0;
 778         int n = (Math.min(certainty, Integer.MAX_VALUE-1)+1)/2;
 779 
 780         // The relationship between the certainty and the number of rounds
 781         // we perform is given in the draft standard ANSI X9.80, "PRIME
 782         // NUMBER GENERATION, PRIMALITY TESTING, AND PRIMALITY CERTIFICATES".
 783         int sizeInBits = this.bitLength();
 784         if (sizeInBits < 100) {
 785             rounds = 50;
 786             rounds = n < rounds ? n : rounds;
 787             return passesMillerRabin(rounds, random);
 788         }
 789 
 790         if (sizeInBits < 256) {
 791             rounds = 27;
 792         } else if (sizeInBits < 512) {
 793             rounds = 15;
 794         } else if (sizeInBits < 768) {
 795             rounds = 8;
 796         } else if (sizeInBits < 1024) {
 797             rounds = 4;
 798         } else {
 799             rounds = 2;
 800         }
 801         rounds = n < rounds ? n : rounds;
 802 
 803         return passesMillerRabin(rounds, random) && passesLucasLehmer();
 804     }
 805 
 806     /**
 807      * Returns true iff this BigInteger is a Lucas-Lehmer probable prime.
 808      *
 809      * The following assumptions are made:
 810      * This BigInteger is a positive, odd number.
 811      */
 812     private boolean passesLucasLehmer() {
 813         BigInteger thisPlusOne = this.add(ONE);
 814 
 815         // Step 1
 816         int d = 5;
 817         while (jacobiSymbol(d, this) != -1) {
 818             // 5, -7, 9, -11, ...
 819             d = (d<0) ? Math.abs(d)+2 : -(d+2);
 820         }
 821 
 822         // Step 2
 823         BigInteger u = lucasLehmerSequence(d, thisPlusOne, this);
 824 
 825         // Step 3
 826         return u.mod(this).equals(ZERO);
 827     }
 828 
 829     /**
 830      * Computes Jacobi(p,n).
 831      * Assumes n positive, odd, n>=3.
 832      */
 833     private static int jacobiSymbol(int p, BigInteger n) {
 834         if (p == 0)
 835             return 0;
 836 
 837         // Algorithm and comments adapted from Colin Plumb's C library.
 838         int j = 1;
 839         int u = n.mag[n.mag.length-1];
 840 
 841         // Make p positive
 842         if (p < 0) {
 843             p = -p;
 844             int n8 = u & 7;
 845             if ((n8 == 3) || (n8 == 7))
 846                 j = -j; // 3 (011) or 7 (111) mod 8
 847         }
 848 
 849         // Get rid of factors of 2 in p
 850         while ((p & 3) == 0)
 851             p >>= 2;
 852         if ((p & 1) == 0) {
 853             p >>= 1;
 854             if (((u ^ (u>>1)) & 2) != 0)
 855                 j = -j; // 3 (011) or 5 (101) mod 8
 856         }
 857         if (p == 1)
 858             return j;
 859         // Then, apply quadratic reciprocity
 860         if ((p & u & 2) != 0)   // p = u = 3 (mod 4)?
 861             j = -j;
 862         // And reduce u mod p
 863         u = n.mod(BigInteger.valueOf(p)).intValue();
 864 
 865         // Now compute Jacobi(u,p), u < p
 866         while (u != 0) {
 867             while ((u & 3) == 0)
 868                 u >>= 2;
 869             if ((u & 1) == 0) {
 870                 u >>= 1;
 871                 if (((p ^ (p>>1)) & 2) != 0)
 872                     j = -j;     // 3 (011) or 5 (101) mod 8
 873             }
 874             if (u == 1)
 875                 return j;
 876             // Now both u and p are odd, so use quadratic reciprocity
 877             assert (u < p);
 878             int t = u; u = p; p = t;
 879             if ((u & p & 2) != 0) // u = p = 3 (mod 4)?
 880                 j = -j;
 881             // Now u >= p, so it can be reduced
 882             u %= p;
 883         }
 884         return 0;
 885     }
 886 
 887     private static BigInteger lucasLehmerSequence(int z, BigInteger k, BigInteger n) {
 888         BigInteger d = BigInteger.valueOf(z);
 889         BigInteger u = ONE; BigInteger u2;
 890         BigInteger v = ONE; BigInteger v2;
 891 
 892         for (int i=k.bitLength()-2; i>=0; i--) {
 893             u2 = u.multiply(v).mod(n);
 894 
 895             v2 = v.square().add(d.multiply(u.square())).mod(n);
 896             if (v2.testBit(0))
 897                 v2 = v2.subtract(n);
 898 
 899             v2 = v2.shiftRight(1);
 900 
 901             u = u2; v = v2;
 902             if (k.testBit(i)) {
 903                 u2 = u.add(v).mod(n);
 904                 if (u2.testBit(0))
 905                     u2 = u2.subtract(n);
 906 
 907                 u2 = u2.shiftRight(1);
 908                 v2 = v.add(d.multiply(u)).mod(n);
 909                 if (v2.testBit(0))
 910                     v2 = v2.subtract(n);
 911                 v2 = v2.shiftRight(1);
 912 
 913                 u = u2; v = v2;
 914             }
 915         }
 916         return u;
 917     }
 918 
 919     private static volatile Random staticRandom;
 920 
 921     private static Random getSecureRandom() {
 922         if (staticRandom == null) {
 923             staticRandom = new java.security.SecureRandom();
 924         }
 925         return staticRandom;
 926     }
 927 
 928     /**
 929      * Returns true iff this BigInteger passes the specified number of
 930      * Miller-Rabin tests. This test is taken from the DSA spec (NIST FIPS
 931      * 186-2).
 932      *
 933      * The following assumptions are made:
 934      * This BigInteger is a positive, odd number greater than 2.
 935      * iterations<=50.
 936      */
 937     private boolean passesMillerRabin(int iterations, Random rnd) {
 938         // Find a and m such that m is odd and this == 1 + 2**a * m
 939         BigInteger thisMinusOne = this.subtract(ONE);
 940         BigInteger m = thisMinusOne;
 941         int a = m.getLowestSetBit();
 942         m = m.shiftRight(a);
 943 
 944         // Do the tests
 945         if (rnd == null) {
 946             rnd = getSecureRandom();
 947         }
 948         for (int i=0; i<iterations; i++) {
 949             // Generate a uniform random on (1, this)
 950             BigInteger b;
 951             do {
 952                 b = new BigInteger(this.bitLength(), rnd);
 953             } while (b.compareTo(ONE) <= 0 || b.compareTo(this) >= 0);
 954 
 955             int j = 0;
 956             BigInteger z = b.modPow(m, this);
 957             while(!((j==0 && z.equals(ONE)) || z.equals(thisMinusOne))) {
 958                 if (j>0 && z.equals(ONE) || ++j==a)
 959                     return false;
 960                 z = z.modPow(TWO, this);
 961             }
 962         }
 963         return true;
 964     }
 965 
 966     /**
 967      * This internal constructor differs from its public cousin
 968      * with the arguments reversed in two ways: it assumes that its
 969      * arguments are correct, and it doesn't copy the magnitude array.
 970      */
 971     BigInteger(int[] magnitude, int signum) {
 972         this.signum = (magnitude.length==0 ? 0 : signum);
 973         this.mag = magnitude;
 974     }
 975 
 976     /**
 977      * This private constructor is for internal use and assumes that its
 978      * arguments are correct.
 979      */
 980     private BigInteger(byte[] magnitude, int signum) {
 981         this.signum = (magnitude.length==0 ? 0 : signum);
 982         this.mag = stripLeadingZeroBytes(magnitude);
 983     }
 984 
 985     //Static Factory Methods
 986 
 987     /**
 988      * Returns a BigInteger whose value is equal to that of the
 989      * specified {@code long}.  This "static factory method" is
 990      * provided in preference to a ({@code long}) constructor
 991      * because it allows for reuse of frequently used BigIntegers.
 992      *
 993      * @param  val value of the BigInteger to return.
 994      * @return a BigInteger with the specified value.
 995      */
 996     public static BigInteger valueOf(long val) {
 997         // If -MAX_CONSTANT < val < MAX_CONSTANT, return stashed constant
 998         if (val == 0)
 999             return ZERO;
1000         if (val > 0 && val <= MAX_CONSTANT)
1001             return posConst[(int) val];
1002         else if (val < 0 && val >= -MAX_CONSTANT)
1003             return negConst[(int) -val];
1004 
1005         return new BigInteger(val);
1006     }
1007 
1008     /**
1009      * Constructs a BigInteger with the specified value, which may not be zero.
1010      */
1011     private BigInteger(long val) {
1012         if (val < 0) {
1013             val = -val;
1014             signum = -1;
1015         } else {
1016             signum = 1;
1017         }
1018 
1019         int highWord = (int)(val >>> 32);
1020         if (highWord==0) {
1021             mag = new int[1];
1022             mag[0] = (int)val;
1023         } else {
1024             mag = new int[2];
1025             mag[0] = highWord;
1026             mag[1] = (int)val;
1027         }
1028     }
1029 
1030     /**
1031      * Returns a BigInteger with the given two's complement representation.
1032      * Assumes that the input array will not be modified (the returned
1033      * BigInteger will reference the input array if feasible).
1034      */
1035     private static BigInteger valueOf(int val[]) {
1036         return (val[0]>0 ? new BigInteger(val, 1) : new BigInteger(val));
1037     }
1038 
1039     // Constants
1040 
1041     /**
1042      * Initialize static constant array when class is loaded.
1043      */
1044     private final static int MAX_CONSTANT = 16;
1045     private static BigInteger posConst[] = new BigInteger[MAX_CONSTANT+1];
1046     private static BigInteger negConst[] = new BigInteger[MAX_CONSTANT+1];
1047 
1048     /**
1049      * The cache of powers of each radix.  This allows us to not have to
1050      * recalculate powers of radix^(2^n) more than once.  This speeds
1051      * Schoenhage recursive base conversion significantly.
1052      */
1053     private static volatile BigInteger[][] powerCache;
1054 
1055     /** The cache of logarithms of radices for base conversion. */
1056     private static final double[] logCache;
1057 
1058     /** The natural log of 2.  This is used in computing cache indices. */
1059     private static final double LOG_TWO = Math.log(2.0);
1060 
1061     static {
1062         for (int i = 1; i <= MAX_CONSTANT; i++) {
1063             int[] magnitude = new int[1];
1064             magnitude[0] = i;
1065             posConst[i] = new BigInteger(magnitude,  1);
1066             negConst[i] = new BigInteger(magnitude, -1);
1067         }
1068 
1069         /*
1070          * Initialize the cache of radix^(2^x) values used for base conversion
1071          * with just the very first value.  Additional values will be created
1072          * on demand.
1073          */
1074         powerCache = new BigInteger[Character.MAX_RADIX+1][];
1075         logCache = new double[Character.MAX_RADIX+1];
1076 
1077         for (int i=Character.MIN_RADIX; i<=Character.MAX_RADIX; i++)
1078         {
1079             powerCache[i] = new BigInteger[] { BigInteger.valueOf(i) };
1080             logCache[i] = Math.log(i);
1081         }
1082     }
1083 
1084     /**
1085      * The BigInteger constant zero.
1086      *
1087      * @since   1.2
1088      */
1089     public static final BigInteger ZERO = new BigInteger(new int[0], 0);
1090 
1091     /**
1092      * The BigInteger constant one.
1093      *
1094      * @since   1.2
1095      */
1096     public static final BigInteger ONE = valueOf(1);
1097 
1098     /**
1099      * The BigInteger constant two.  (Not exported.)
1100      */
1101     private static final BigInteger TWO = valueOf(2);
1102 
1103     /**
1104      * The BigInteger constant -1.  (Not exported.)
1105      */
1106     private static final BigInteger NEGATIVE_ONE = valueOf(-1);
1107 
1108     /**
1109      * The BigInteger constant ten.
1110      *
1111      * @since   1.5
1112      */
1113     public static final BigInteger TEN = valueOf(10);
1114 
1115     // Arithmetic Operations
1116 
1117     /**
1118      * Returns a BigInteger whose value is {@code (this + val)}.
1119      *
1120      * @param  val value to be added to this BigInteger.
1121      * @return {@code this + val}
1122      */
1123     public BigInteger add(BigInteger val) {
1124         if (val.signum == 0)
1125             return this;
1126         if (signum == 0)
1127             return val;
1128         if (val.signum == signum)
1129             return new BigInteger(add(mag, val.mag), signum);
1130 
1131         int cmp = compareMagnitude(val);
1132         if (cmp == 0)
1133             return ZERO;
1134         int[] resultMag = (cmp > 0 ? subtract(mag, val.mag)
1135                            : subtract(val.mag, mag));
1136         resultMag = trustedStripLeadingZeroInts(resultMag);
1137 
1138         return new BigInteger(resultMag, cmp == signum ? 1 : -1);
1139     }
1140 
1141     /**
1142      * Package private methods used by BigDecimal code to add a BigInteger
1143      * with a long. Assumes val is not equal to INFLATED.
1144      */
1145     BigInteger add(long val) {
1146         if (val == 0)
1147             return this;
1148         if (signum == 0)
1149             return valueOf(val);
1150         if (Long.signum(val) == signum)
1151             return new BigInteger(add(mag, Math.abs(val)), signum);
1152         int cmp = compareMagnitude(val);
1153         if (cmp == 0)
1154             return ZERO;
1155         int[] resultMag = (cmp > 0 ? subtract(mag, Math.abs(val)) : subtract(Math.abs(val), mag));
1156         resultMag = trustedStripLeadingZeroInts(resultMag);
1157         return new BigInteger(resultMag, cmp == signum ? 1 : -1);
1158     }
1159 
1160     /**
1161      * Adds the contents of the int array x and long value val. This
1162      * method allocates a new int array to hold the answer and returns
1163      * a reference to that array.  Assumes x.length &gt; 0 and val is
1164      * non-negative
1165      */
1166     private static int[] add(int[] x, long val) {
1167         int[] y;
1168         long sum = 0;
1169         int xIndex = x.length;
1170         int[] result;
1171         int highWord = (int)(val >>> 32);
1172         if (highWord==0) {
1173             result = new int[xIndex];
1174             sum = (x[--xIndex] & LONG_MASK) + val;
1175             result[xIndex] = (int)sum;
1176         } else {
1177             if (xIndex == 1) {
1178                 result = new int[2];
1179                 sum = val  + (x[0] & LONG_MASK);
1180                 result[1] = (int)sum;
1181                 result[0] = (int)(sum >>> 32);
1182                 return result;
1183             } else {
1184                 result = new int[xIndex];
1185                 sum = (x[--xIndex] & LONG_MASK) + (val & LONG_MASK);
1186                 result[xIndex] = (int)sum;
1187                 sum = (x[--xIndex] & LONG_MASK) + (highWord & LONG_MASK) + (sum >>> 32);
1188                 result[xIndex] = (int)sum;
1189             }
1190         }
1191         // Copy remainder of longer number while carry propagation is required
1192         boolean carry = (sum >>> 32 != 0);
1193         while (xIndex > 0 && carry)
1194             carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
1195         // Copy remainder of longer number
1196         while (xIndex > 0)
1197             result[--xIndex] = x[xIndex];
1198         // Grow result if necessary
1199         if (carry) {
1200             int bigger[] = new int[result.length + 1];
1201             System.arraycopy(result, 0, bigger, 1, result.length);
1202             bigger[0] = 0x01;
1203             return bigger;
1204         }
1205         return result;
1206     }
1207 
1208     /**
1209      * Adds the contents of the int arrays x and y. This method allocates
1210      * a new int array to hold the answer and returns a reference to that
1211      * array.
1212      */
1213     private static int[] add(int[] x, int[] y) {
1214         // If x is shorter, swap the two arrays
1215         if (x.length < y.length) {
1216             int[] tmp = x;
1217             x = y;
1218             y = tmp;
1219         }
1220 
1221         int xIndex = x.length;
1222         int yIndex = y.length;
1223         int result[] = new int[xIndex];
1224         long sum = 0;
1225         if(yIndex==1) {
1226             sum = (x[--xIndex] & LONG_MASK) + (y[0] & LONG_MASK) ;
1227             result[xIndex] = (int)sum;
1228         } else {
1229             // Add common parts of both numbers
1230             while(yIndex > 0) {
1231                 sum = (x[--xIndex] & LONG_MASK) +
1232                       (y[--yIndex] & LONG_MASK) + (sum >>> 32);
1233                 result[xIndex] = (int)sum;
1234             }
1235         }
1236         // Copy remainder of longer number while carry propagation is required
1237         boolean carry = (sum >>> 32 != 0);
1238         while (xIndex > 0 && carry)
1239             carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
1240 
1241         // Copy remainder of longer number
1242         while (xIndex > 0)
1243             result[--xIndex] = x[xIndex];
1244 
1245         // Grow result if necessary
1246         if (carry) {
1247             int bigger[] = new int[result.length + 1];
1248             System.arraycopy(result, 0, bigger, 1, result.length);
1249             bigger[0] = 0x01;
1250             return bigger;
1251         }
1252         return result;
1253     }
1254 
1255     private static int[] subtract(long val, int[] little) {
1256         int highWord = (int)(val >>> 32);
1257         if (highWord==0) {
1258             int result[] = new int[1];
1259             result[0] = (int)(val - (little[0] & LONG_MASK));
1260             return result;
1261         } else {
1262             int result[] = new int[2];
1263             if(little.length==1) {
1264                 long difference = ((int)val & LONG_MASK) - (little[0] & LONG_MASK);
1265                 result[1] = (int)difference;
1266                 // Subtract remainder of longer number while borrow propagates
1267                 boolean borrow = (difference >> 32 != 0);
1268                 if(borrow) {
1269                     result[0] = highWord - 1;
1270                 } else {        // Copy remainder of longer number
1271                     result[0] = highWord;
1272                 }
1273                 return result;
1274             } else { // little.length==2
1275                 long difference = ((int)val & LONG_MASK) - (little[1] & LONG_MASK);
1276                 result[1] = (int)difference;
1277                 difference = (highWord & LONG_MASK) - (little[0] & LONG_MASK) + (difference >> 32);
1278                 result[0] = (int)difference;
1279                 return result;
1280             }
1281         }
1282     }
1283 
1284     /**
1285      * Subtracts the contents of the second argument (val) from the
1286      * first (big).  The first int array (big) must represent a larger number
1287      * than the second.  This method allocates the space necessary to hold the
1288      * answer.
1289      * assumes val &gt;= 0
1290      */
1291     private static int[] subtract(int[] big, long val) {
1292         int highWord = (int)(val >>> 32);
1293         int bigIndex = big.length;
1294         int result[] = new int[bigIndex];
1295         long difference = 0;
1296 
1297         if (highWord==0) {
1298             difference = (big[--bigIndex] & LONG_MASK) - val;
1299             result[bigIndex] = (int)difference;
1300         } else {
1301             difference = (big[--bigIndex] & LONG_MASK) - (val & LONG_MASK);
1302             result[bigIndex] = (int)difference;
1303             difference = (big[--bigIndex] & LONG_MASK) - (highWord & LONG_MASK) + (difference >> 32);
1304             result[bigIndex] = (int)difference;
1305         }
1306 
1307 
1308         // Subtract remainder of longer number while borrow propagates
1309         boolean borrow = (difference >> 32 != 0);
1310         while (bigIndex > 0 && borrow)
1311             borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
1312 
1313         // Copy remainder of longer number
1314         while (bigIndex > 0)
1315             result[--bigIndex] = big[bigIndex];
1316 
1317         return result;
1318     }
1319 
1320     /**
1321      * Returns a BigInteger whose value is {@code (this - val)}.
1322      *
1323      * @param  val value to be subtracted from this BigInteger.
1324      * @return {@code this - val}
1325      */
1326     public BigInteger subtract(BigInteger val) {
1327         if (val.signum == 0)
1328             return this;
1329         if (signum == 0)
1330             return val.negate();
1331         if (val.signum != signum)
1332             return new BigInteger(add(mag, val.mag), signum);
1333 
1334         int cmp = compareMagnitude(val);
1335         if (cmp == 0)
1336             return ZERO;
1337         int[] resultMag = (cmp > 0 ? subtract(mag, val.mag)
1338                            : subtract(val.mag, mag));
1339         resultMag = trustedStripLeadingZeroInts(resultMag);
1340         return new BigInteger(resultMag, cmp == signum ? 1 : -1);
1341     }
1342 
1343     /**
1344      * Subtracts the contents of the second int arrays (little) from the
1345      * first (big).  The first int array (big) must represent a larger number
1346      * than the second.  This method allocates the space necessary to hold the
1347      * answer.
1348      */
1349     private static int[] subtract(int[] big, int[] little) {
1350         int bigIndex = big.length;
1351         int result[] = new int[bigIndex];
1352         int littleIndex = little.length;
1353         long difference = 0;
1354 
1355         // Subtract common parts of both numbers
1356         while(littleIndex > 0) {
1357             difference = (big[--bigIndex] & LONG_MASK) -
1358                          (little[--littleIndex] & LONG_MASK) +
1359                          (difference >> 32);
1360             result[bigIndex] = (int)difference;
1361         }
1362 
1363         // Subtract remainder of longer number while borrow propagates
1364         boolean borrow = (difference >> 32 != 0);
1365         while (bigIndex > 0 && borrow)
1366             borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
1367 
1368         // Copy remainder of longer number
1369         while (bigIndex > 0)
1370             result[--bigIndex] = big[bigIndex];
1371 
1372         return result;
1373     }
1374 
1375     /**
1376      * Returns a BigInteger whose value is {@code (this * val)}.
1377      *
1378      * @param  val value to be multiplied by this BigInteger.
1379      * @return {@code this * val}
1380      */
1381     public BigInteger multiply(BigInteger val) {
1382         if (val.signum == 0 || signum == 0)
1383             return ZERO;
1384 
1385         int xlen = mag.length;
1386         int ylen = val.mag.length;
1387 
1388         if ((xlen < KARATSUBA_THRESHOLD) || (ylen < KARATSUBA_THRESHOLD))
1389         {
1390             int resultSign = signum == val.signum ? 1 : -1;
1391             if (val.mag.length == 1) {
1392                 return multiplyByInt(mag,val.mag[0], resultSign);
1393             }
1394             if(mag.length == 1) {
1395                 return multiplyByInt(val.mag,mag[0], resultSign);
1396             }
1397             int[] result = multiplyToLen(mag, xlen,
1398                                          val.mag, ylen, null);
1399             result = trustedStripLeadingZeroInts(result);
1400             return new BigInteger(result, resultSign);
1401         }
1402         else
1403             if ((xlen < TOOM_COOK_THRESHOLD) && (ylen < TOOM_COOK_THRESHOLD))
1404                 return multiplyKaratsuba(this, val);
1405             else
1406                 return multiplyToomCook3(this, val);
1407     }
1408 
1409     private static BigInteger multiplyByInt(int[] x, int y, int sign) {
1410         if(Integer.bitCount(y)==1) {
1411             return new BigInteger(shiftLeft(x,Integer.numberOfTrailingZeros(y)), sign);
1412         }
1413         int xlen = x.length;
1414         int[] rmag =  new int[xlen + 1];
1415         long carry = 0;
1416         long yl = y & LONG_MASK;
1417         int rstart = rmag.length - 1;
1418         for (int i = xlen - 1; i >= 0; i--) {
1419             long product = (x[i] & LONG_MASK) * yl + carry;
1420             rmag[rstart--] = (int)product;
1421             carry = product >>> 32;
1422         }
1423         if (carry == 0L) {
1424             rmag = java.util.Arrays.copyOfRange(rmag, 1, rmag.length);
1425         } else {
1426             rmag[rstart] = (int)carry;
1427         }
1428         return new BigInteger(rmag, sign);
1429     }
1430 
1431     /**
1432      * Package private methods used by BigDecimal code to multiply a BigInteger
1433      * with a long. Assumes v is not equal to INFLATED.
1434      */
1435     BigInteger multiply(long v) {
1436         if (v == 0 || signum == 0)
1437           return ZERO;
1438         if (v == BigDecimal.INFLATED)
1439             return multiply(BigInteger.valueOf(v));
1440         int rsign = (v > 0 ? signum : -signum);
1441         if (v < 0)
1442             v = -v;
1443         long dh = v >>> 32;      // higher order bits
1444         long dl = v & LONG_MASK; // lower order bits
1445 
1446         int xlen = mag.length;
1447         int[] value = mag;
1448         int[] rmag = (dh == 0L) ? (new int[xlen + 1]) : (new int[xlen + 2]);
1449         long carry = 0;
1450         int rstart = rmag.length - 1;
1451         for (int i = xlen - 1; i >= 0; i--) {
1452             long product = (value[i] & LONG_MASK) * dl + carry;
1453             rmag[rstart--] = (int)product;
1454             carry = product >>> 32;
1455         }
1456         rmag[rstart] = (int)carry;
1457         if (dh != 0L) {
1458             carry = 0;
1459             rstart = rmag.length - 2;
1460             for (int i = xlen - 1; i >= 0; i--) {
1461                 long product = (value[i] & LONG_MASK) * dh +
1462                     (rmag[rstart] & LONG_MASK) + carry;
1463                 rmag[rstart--] = (int)product;
1464                 carry = product >>> 32;
1465             }
1466             rmag[0] = (int)carry;
1467         }
1468         if (carry == 0L)
1469             rmag = java.util.Arrays.copyOfRange(rmag, 1, rmag.length);
1470         return new BigInteger(rmag, rsign);
1471     }
1472 
1473     /**
1474      * Multiplies int arrays x and y to the specified lengths and places
1475      * the result into z. There will be no leading zeros in the resultant array.
1476      */
1477     private int[] multiplyToLen(int[] x, int xlen, int[] y, int ylen, int[] z) {
1478         int xstart = xlen - 1;
1479         int ystart = ylen - 1;
1480 
1481         if (z == null || z.length < (xlen+ ylen))
1482             z = new int[xlen+ylen];
1483 
1484         long carry = 0;
1485         for (int j=ystart, k=ystart+1+xstart; j>=0; j--, k--) {
1486             long product = (y[j] & LONG_MASK) *
1487                            (x[xstart] & LONG_MASK) + carry;
1488             z[k] = (int)product;
1489             carry = product >>> 32;
1490         }
1491         z[xstart] = (int)carry;
1492 
1493         for (int i = xstart-1; i >= 0; i--) {
1494             carry = 0;
1495             for (int j=ystart, k=ystart+1+i; j>=0; j--, k--) {
1496                 long product = (y[j] & LONG_MASK) *
1497                                (x[i] & LONG_MASK) +
1498                                (z[k] & LONG_MASK) + carry;
1499                 z[k] = (int)product;
1500                 carry = product >>> 32;
1501             }
1502             z[i] = (int)carry;
1503         }
1504         return z;
1505     }
1506 
1507     /**
1508      * Multiplies two BigIntegers using the Karatsuba multiplication
1509      * algorithm.  This is a recursive divide-and-conquer algorithm which is
1510      * more efficient for large numbers than what is commonly called the
1511      * "grade-school" algorithm used in multiplyToLen.  If the numbers to be
1512      * multiplied have length n, the "grade-school" algorithm has an
1513      * asymptotic complexity of O(n^2).  In contrast, the Karatsuba algorithm
1514      * has complexity of O(n^(log2(3))), or O(n^1.585).  It achieves this
1515      * increased performance by doing 3 multiplies instead of 4 when
1516      * evaluating the product.  As it has some overhead, should be used when
1517      * both numbers are larger than a certain threshold (found
1518      * experimentally).
1519      *
1520      * See:  http://en.wikipedia.org/wiki/Karatsuba_algorithm
1521      */
1522     private static BigInteger multiplyKaratsuba(BigInteger x, BigInteger y)
1523     {
1524         int xlen = x.mag.length;
1525         int ylen = y.mag.length;
1526 
1527         // The number of ints in each half of the number.
1528         int half = (Math.max(xlen, ylen)+1) / 2;
1529 
1530         // xl and yl are the lower halves of x and y respectively,
1531         // xh and yh are the upper halves.
1532         BigInteger xl = x.getLower(half);
1533         BigInteger xh = x.getUpper(half);
1534         BigInteger yl = y.getLower(half);
1535         BigInteger yh = y.getUpper(half);
1536 
1537         BigInteger p1 = xh.multiply(yh);  // p1 = xh*yh
1538         BigInteger p2 = xl.multiply(yl);  // p2 = xl*yl
1539 
1540         // p3=(xh+xl)*(yh+yl)
1541         BigInteger p3 = xh.add(xl).multiply(yh.add(yl));
1542 
1543         // result = p1 * 2^(32*2*half) + (p3 - p1 - p2) * 2^(32*half) + p2
1544         BigInteger result = p1.shiftLeft(32*half).add(p3.subtract(p1).subtract(p2)).shiftLeft(32*half).add(p2);
1545 
1546         if (x.signum != y.signum)
1547             return result.negate();
1548         else
1549             return result;
1550     }
1551 
1552     /**
1553      * Multiplies two BigIntegers using a 3-way Toom-Cook multiplication
1554      * algorithm.  This is a recursive divide-and-conquer algorithm which is
1555      * more efficient for large numbers than what is commonly called the
1556      * "grade-school" algorithm used in multiplyToLen.  If the numbers to be
1557      * multiplied have length n, the "grade-school" algorithm has an
1558      * asymptotic complexity of O(n^2).  In contrast, 3-way Toom-Cook has a
1559      * complexity of about O(n^1.465).  It achieves this increased asymptotic
1560      * performance by breaking each number into three parts and by doing 5
1561      * multiplies instead of 9 when evaluating the product.  Due to overhead
1562      * (additions, shifts, and one division) in the Toom-Cook algorithm, it
1563      * should only be used when both numbers are larger than a certain
1564      * threshold (found experimentally).  This threshold is generally larger
1565      * than that for Karatsuba multiplication, so this algorithm is generally
1566      * only used when numbers become significantly larger.
1567      *
1568      * The algorithm used is the "optimal" 3-way Toom-Cook algorithm outlined
1569      * by Marco Bodrato.
1570      *
1571      *  See: http://bodrato.it/toom-cook/
1572      *       http://bodrato.it/papers/#WAIFI2007
1573      *
1574      * "Towards Optimal Toom-Cook Multiplication for Univariate and
1575      * Multivariate Polynomials in Characteristic 2 and 0." by Marco BODRATO;
1576      * In C.Carlet and B.Sunar, Eds., "WAIFI'07 proceedings", p. 116-133,
1577      * LNCS #4547. Springer, Madrid, Spain, June 21-22, 2007.
1578      *
1579      */
1580     private static BigInteger multiplyToomCook3(BigInteger a, BigInteger b)
1581     {
1582         int alen = a.mag.length;
1583         int blen = b.mag.length;
1584 
1585         int largest = Math.max(alen, blen);
1586 
1587         // k is the size (in ints) of the lower-order slices.
1588         int k = (largest+2)/3;   // Equal to ceil(largest/3)
1589 
1590         // r is the size (in ints) of the highest-order slice.
1591         int r = largest - 2*k;
1592 
1593         // Obtain slices of the numbers. a2 and b2 are the most significant
1594         // bits of the numbers a and b, and a0 and b0 the least significant.
1595         BigInteger a0, a1, a2, b0, b1, b2;
1596         a2 = a.getToomSlice(k, r, 0, largest);
1597         a1 = a.getToomSlice(k, r, 1, largest);
1598         a0 = a.getToomSlice(k, r, 2, largest);
1599         b2 = b.getToomSlice(k, r, 0, largest);
1600         b1 = b.getToomSlice(k, r, 1, largest);
1601         b0 = b.getToomSlice(k, r, 2, largest);
1602 
1603         BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1, db1;
1604 
1605         v0 = a0.multiply(b0);
1606         da1 = a2.add(a0);
1607         db1 = b2.add(b0);
1608         vm1 = da1.subtract(a1).multiply(db1.subtract(b1));
1609         da1 = da1.add(a1);
1610         db1 = db1.add(b1);
1611         v1 = da1.multiply(db1);
1612         v2 = da1.add(a2).shiftLeft(1).subtract(a0).multiply(
1613              db1.add(b2).shiftLeft(1).subtract(b0));
1614         vinf = a2.multiply(b2);
1615 
1616         /* The algorithm requires two divisions by 2 and one by 3.
1617            All divisions are known to be exact, that is, they do not produce
1618            remainders, and all results are positive.  The divisions by 2 are
1619            implemented as right shifts which are relatively efficient, leaving
1620            only an exact division by 3, which is done by a specialized
1621            linear-time algorithm. */
1622         t2 = v2.subtract(vm1).exactDivideBy3();
1623         tm1 = v1.subtract(vm1).shiftRight(1);
1624         t1 = v1.subtract(v0);
1625         t2 = t2.subtract(t1).shiftRight(1);
1626         t1 = t1.subtract(tm1).subtract(vinf);
1627         t2 = t2.subtract(vinf.shiftLeft(1));
1628         tm1 = tm1.subtract(t2);
1629 
1630         // Number of bits to shift left.
1631         int ss = k*32;
1632 
1633         BigInteger result = vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
1634 
1635         if (a.signum != b.signum)
1636             return result.negate();
1637         else
1638             return result;
1639     }
1640 
1641 
1642     /**
1643      * Returns a slice of a BigInteger for use in Toom-Cook multiplication.
1644      *
1645      * @param lowerSize The size of the lower-order bit slices.
1646      * @param upperSize The size of the higher-order bit slices.
1647      * @param slice The index of which slice is requested, which must be a
1648      * number from 0 to size-1. Slice 0 is the highest-order bits, and slice
1649      * size-1 are the lowest-order bits. Slice 0 may be of different size than
1650      * the other slices.
1651      * @param fullsize The size of the larger integer array, used to align
1652      * slices to the appropriate position when multiplying different-sized
1653      * numbers.
1654      */
1655     private BigInteger getToomSlice(int lowerSize, int upperSize, int slice,
1656                                     int fullsize)
1657     {
1658         int start, end, sliceSize, len, offset;
1659 
1660         len = mag.length;
1661         offset = fullsize - len;
1662 
1663         if (slice == 0)
1664         {
1665             start = 0 - offset;
1666             end = upperSize - 1 - offset;
1667         }
1668         else
1669         {
1670             start = upperSize + (slice-1)*lowerSize - offset;
1671             end = start + lowerSize - 1;
1672         }
1673 
1674         if (start < 0)
1675             start = 0;
1676         if (end < 0)
1677            return ZERO;
1678 
1679         sliceSize = (end-start) + 1;
1680 
1681         if (sliceSize <= 0)
1682             return ZERO;
1683 
1684         // While performing Toom-Cook, all slices are positive and
1685         // the sign is adjusted when the final number is composed.
1686         if (start==0 && sliceSize >= len)
1687             return this.abs();
1688 
1689         int intSlice[] = new int[sliceSize];
1690         System.arraycopy(mag, start, intSlice, 0, sliceSize);
1691 
1692         return new BigInteger(trustedStripLeadingZeroInts(intSlice), 1);
1693     }
1694 
1695     /**
1696      * Does an exact division (that is, the remainder is known to be zero)
1697      * of the specified number by 3.  This is used in Toom-Cook
1698      * multiplication.  This is an efficient algorithm that runs in linear
1699      * time.  If the argument is not exactly divisible by 3, results are
1700      * undefined.  Note that this is expected to be called with positive
1701      * arguments only.
1702      */
1703     private BigInteger exactDivideBy3()
1704     {
1705         int len = mag.length;
1706         int[] result = new int[len];
1707         long x, w, q, borrow;
1708         borrow = 0L;
1709         for (int i=len-1; i>=0; i--)
1710         {
1711             x = (mag[i] & LONG_MASK);
1712             w = x - borrow;
1713             if (borrow > x)       // Did we make the number go negative?
1714                 borrow = 1L;
1715             else
1716                 borrow = 0L;
1717 
1718             // 0xAAAAAAAB is the modular inverse of 3 (mod 2^32).  Thus,
1719             // the effect of this is to divide by 3 (mod 2^32).
1720             // This is much faster than division on most architectures.
1721             q = (w * 0xAAAAAAABL) & LONG_MASK;
1722             result[i] = (int) q;
1723 
1724             // Now check the borrow. The second check can of course be
1725             // eliminated if the first fails.
1726             if (q >= 0x55555556L)
1727             {
1728                 borrow++;
1729                 if (q >= 0xAAAAAAABL)
1730                     borrow++;
1731             }
1732         }
1733         result = trustedStripLeadingZeroInts(result);
1734         return new BigInteger(result, signum);
1735     }
1736 
1737     /**
1738      * Returns a new BigInteger representing n lower ints of the number.
1739      * This is used by Karatsuba multiplication and Karatsuba squaring.
1740      */
1741     private BigInteger getLower(int n) {
1742         int len = mag.length;
1743 
1744         if (len <= n)
1745             return this;
1746 
1747         int lowerInts[] = new int[n];
1748         System.arraycopy(mag, len-n, lowerInts, 0, n);
1749 
1750         return new BigInteger(trustedStripLeadingZeroInts(lowerInts), 1);
1751     }
1752 
1753     /**
1754      * Returns a new BigInteger representing mag.length-n upper
1755      * ints of the number.  This is used by Karatsuba multiplication and
1756      * Karatsuba squaring.
1757      */
1758     private BigInteger getUpper(int n) {
1759         int len = mag.length;
1760 
1761         if (len <= n)
1762             return ZERO;
1763 
1764         int upperLen = len - n;
1765         int upperInts[] = new int[upperLen];
1766         System.arraycopy(mag, 0, upperInts, 0, upperLen);
1767 
1768         return new BigInteger(trustedStripLeadingZeroInts(upperInts), 1);
1769     }
1770 
1771     // Squaring
1772 
1773     /**
1774      * Returns a BigInteger whose value is {@code (this<sup>2</sup>)}.
1775      *
1776      * @return {@code this<sup>2</sup>}
1777      */
1778     private BigInteger square() {
1779         if (signum == 0)
1780             return ZERO;
1781         int len = mag.length;
1782 
1783         if (len < KARATSUBA_SQUARE_THRESHOLD)
1784         {
1785             int[] z = squareToLen(mag, len, null);
1786             return new BigInteger(trustedStripLeadingZeroInts(z), 1);
1787         }
1788         else
1789             if (len < TOOM_COOK_SQUARE_THRESHOLD)
1790                 return squareKaratsuba();
1791             else
1792                 return squareToomCook3();
1793     }
1794 
1795     /**
1796      * Squares the contents of the int array x. The result is placed into the
1797      * int array z.  The contents of x are not changed.
1798      */
1799     private static final int[] squareToLen(int[] x, int len, int[] z) {
1800         /*
1801          * The algorithm used here is adapted from Colin Plumb's C library.
1802          * Technique: Consider the partial products in the multiplication
1803          * of "abcde" by itself:
1804          *
1805          *               a  b  c  d  e
1806          *            *  a  b  c  d  e
1807          *          ==================
1808          *              ae be ce de ee
1809          *           ad bd cd dd de
1810          *        ac bc cc cd ce
1811          *     ab bb bc bd be
1812          *  aa ab ac ad ae
1813          *
1814          * Note that everything above the main diagonal:
1815          *              ae be ce de = (abcd) * e
1816          *           ad bd cd       = (abc) * d
1817          *        ac bc             = (ab) * c
1818          *     ab                   = (a) * b
1819          *
1820          * is a copy of everything below the main diagonal:
1821          *                       de
1822          *                 cd ce
1823          *           bc bd be
1824          *     ab ac ad ae
1825          *
1826          * Thus, the sum is 2 * (off the diagonal) + diagonal.
1827          *
1828          * This is accumulated beginning with the diagonal (which
1829          * consist of the squares of the digits of the input), which is then
1830          * divided by two, the off-diagonal added, and multiplied by two
1831          * again.  The low bit is simply a copy of the low bit of the
1832          * input, so it doesn't need special care.
1833          */
1834         int zlen = len << 1;
1835         if (z == null || z.length < zlen)
1836             z = new int[zlen];
1837 
1838         // Store the squares, right shifted one bit (i.e., divided by 2)
1839         int lastProductLowWord = 0;
1840         for (int j=0, i=0; j<len; j++) {
1841             long piece = (x[j] & LONG_MASK);
1842             long product = piece * piece;
1843             z[i++] = (lastProductLowWord << 31) | (int)(product >>> 33);
1844             z[i++] = (int)(product >>> 1);
1845             lastProductLowWord = (int)product;
1846         }
1847 
1848         // Add in off-diagonal sums
1849         for (int i=len, offset=1; i>0; i--, offset+=2) {
1850             int t = x[i-1];
1851             t = mulAdd(z, x, offset, i-1, t);
1852             addOne(z, offset-1, i, t);
1853         }
1854 
1855         // Shift back up and set low bit
1856         primitiveLeftShift(z, zlen, 1);
1857         z[zlen-1] |= x[len-1] & 1;
1858 
1859         return z;
1860     }
1861 
1862     /**
1863      * Squares a BigInteger using the Karatsuba squaring algorithm.  It should
1864      * be used when both numbers are larger than a certain threshold (found
1865      * experimentally).  It is a recursive divide-and-conquer algorithm that
1866      * has better asymptotic performance than the algorithm used in
1867      * squareToLen.
1868      */
1869     private BigInteger squareKaratsuba()
1870     {
1871         int half = (mag.length+1) / 2;
1872 
1873         BigInteger xl = getLower(half);
1874         BigInteger xh = getUpper(half);
1875 
1876         BigInteger xhs = xh.square();  // xhs = xh^2
1877         BigInteger xls = xl.square();  // xls = xl^2
1878 
1879         // xh^2 << 64  +  (((xl+xh)^2 - (xh^2 + xl^2)) << 32) + xl^2
1880         return xhs.shiftLeft(half*32).add(xl.add(xh).square().subtract(xhs.add(xls))).shiftLeft(half*32).add(xls);
1881     }
1882 
1883     /**
1884      * Squares a BigInteger using the 3-way Toom-Cook squaring algorithm.  It
1885      * should be used when both numbers are larger than a certain threshold
1886      * (found experimentally).  It is a recursive divide-and-conquer algorithm
1887      * that has better asymptotic performance than the algorithm used in
1888      * squareToLen or squareKaratsuba.
1889      */
1890     private BigInteger squareToomCook3()
1891     {
1892         int len = mag.length;
1893 
1894         // k is the size (in ints) of the lower-order slices.
1895         int k = (len+2)/3;   // Equal to ceil(largest/3)
1896 
1897         // r is the size (in ints) of the highest-order slice.
1898         int r = len - 2*k;
1899 
1900         // Obtain slices of the numbers. a2 is the most significant
1901         // bits of the number, and a0 the least significant.
1902         BigInteger a0, a1, a2;
1903         a2 = getToomSlice(k, r, 0, len);
1904         a1 = getToomSlice(k, r, 1, len);
1905         a0 = getToomSlice(k, r, 2, len);
1906         BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1;
1907 
1908         v0 = a0.square();
1909         da1 = a2.add(a0);
1910         vm1 = da1.subtract(a1).square();
1911         da1 = da1.add(a1);
1912         v1 = da1.square();
1913         vinf = a2.square();
1914         v2 = da1.add(a2).shiftLeft(1).subtract(a0).square();
1915 
1916         /* The algorithm requires two divisions by 2 and one by 3.
1917            All divisions are known to be exact, that is, they do not produce
1918            remainders, and all results are positive.  The divisions by 2 are
1919            implemented as right shifts which are relatively efficient, leaving
1920            only a division by 3.
1921            The division by 3 is done by an optimized algorithm for this case.
1922         */
1923         t2 = v2.subtract(vm1).exactDivideBy3();
1924         tm1 = v1.subtract(vm1).shiftRight(1);
1925         t1 = v1.subtract(v0);
1926         t2 = t2.subtract(t1).shiftRight(1);
1927         t1 = t1.subtract(tm1).subtract(vinf);
1928         t2 = t2.subtract(vinf.shiftLeft(1));
1929         tm1 = tm1.subtract(t2);
1930 
1931         // Number of bits to shift left.
1932         int ss = k*32;
1933 
1934         return vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
1935     }
1936 
1937     // Division
1938 
1939     /**
1940      * Returns a BigInteger whose value is {@code (this / val)}.
1941      *
1942      * @param  val value by which this BigInteger is to be divided.
1943      * @return {@code this / val}
1944      * @throws ArithmeticException if {@code val} is zero.
1945      */
1946     public BigInteger divide(BigInteger val) {
1947         if (mag.length<BURNIKEL_ZIEGLER_THRESHOLD || val.mag.length<BURNIKEL_ZIEGLER_THRESHOLD)
1948             return divideKnuth(val);
1949         else
1950             return divideBurnikelZiegler(val);
1951     }
1952 
1953     /**
1954      * Returns a BigInteger whose value is {@code (this / val)} using an O(n^2) algorithm from Knuth.
1955      *
1956      * @param  val value by which this BigInteger is to be divided.
1957      * @return {@code this / val}
1958      * @throws ArithmeticException if {@code val} is zero.
1959      * @see MutableBigInteger#divideKnuth(MutableBigInteger, MutableBigInteger, boolean)
1960      */
1961     private BigInteger divideKnuth(BigInteger val) {
1962         MutableBigInteger q = new MutableBigInteger(),
1963                           a = new MutableBigInteger(this.mag),
1964                           b = new MutableBigInteger(val.mag);
1965 
1966         a.divideKnuth(b, q, false);
1967         return q.toBigInteger(this.signum * val.signum);
1968     }
1969 
1970     /**
1971      * Returns an array of two BigIntegers containing {@code (this / val)}
1972      * followed by {@code (this % val)}.
1973      *
1974      * @param  val value by which this BigInteger is to be divided, and the
1975      *         remainder computed.
1976      * @return an array of two BigIntegers: the quotient {@code (this / val)}
1977      *         is the initial element, and the remainder {@code (this % val)}
1978      *         is the final element.
1979      * @throws ArithmeticException if {@code val} is zero.
1980      */
1981     public BigInteger[] divideAndRemainder(BigInteger val) {
1982         if (mag.length<BURNIKEL_ZIEGLER_THRESHOLD || val.mag.length<BURNIKEL_ZIEGLER_THRESHOLD)
1983             return divideAndRemainderKnuth(val);
1984         else
1985             return divideAndRemainderBurnikelZiegler(val);
1986     }
1987 
1988     /** Long division */
1989     private BigInteger[] divideAndRemainderKnuth(BigInteger val) {
1990         BigInteger[] result = new BigInteger[2];
1991         MutableBigInteger q = new MutableBigInteger(),
1992                           a = new MutableBigInteger(this.mag),
1993                           b = new MutableBigInteger(val.mag);
1994         MutableBigInteger r = a.divideKnuth(b, q);
1995         result[0] = q.toBigInteger(this.signum == val.signum ? 1 : -1);
1996         result[1] = r.toBigInteger(this.signum);
1997         return result;
1998     }
1999 
2000     /**
2001      * Returns a BigInteger whose value is {@code (this % val)}.
2002      *
2003      * @param  val value by which this BigInteger is to be divided, and the
2004      *         remainder computed.
2005      * @return {@code this % val}
2006      * @throws ArithmeticException if {@code val} is zero.
2007      */
2008     public BigInteger remainder(BigInteger val) {
2009         if (mag.length<BURNIKEL_ZIEGLER_THRESHOLD || val.mag.length<BURNIKEL_ZIEGLER_THRESHOLD)
2010             return remainderKnuth(val);
2011         else
2012             return remainderBurnikelZiegler(val);
2013     }
2014 
2015     /** Long division */
2016     private BigInteger remainderKnuth(BigInteger val) {
2017         MutableBigInteger q = new MutableBigInteger(),
2018                           a = new MutableBigInteger(this.mag),
2019                           b = new MutableBigInteger(val.mag);
2020 
2021         return a.divideKnuth(b, q).toBigInteger(this.signum);
2022     }
2023 
2024     /**
2025      * Calculates {@code this / val} using the Burnikel-Ziegler algorithm.
2026      * @param  val the divisor
2027      * @return {@code this / val}
2028      */
2029     private BigInteger divideBurnikelZiegler(BigInteger val) {
2030         return divideAndRemainderBurnikelZiegler(val)[0];
2031     }
2032 
2033     /**
2034      * Calculates {@code this % val} using the Burnikel-Ziegler algorithm.
2035      * @param val the divisor
2036      * @return {@code this % val}
2037      */
2038     private BigInteger remainderBurnikelZiegler(BigInteger val) {
2039         return divideAndRemainderBurnikelZiegler(val)[1];
2040     }
2041 
2042     /**
2043      * Computes {@code this / val} and {@code this % val} using the
2044      * Burnikel-Ziegler algorithm.
2045      * @param val the divisor
2046      * @return an array containing the quotient and remainder
2047      */
2048     private BigInteger[] divideAndRemainderBurnikelZiegler(BigInteger val) {
2049         MutableBigInteger q = new MutableBigInteger();
2050         MutableBigInteger r = new MutableBigInteger(this).divideAndRemainderBurnikelZiegler(new MutableBigInteger(val), q);
2051         BigInteger qBigInt = q.isZero() ? ZERO : q.toBigInteger(signum*val.signum);
2052         BigInteger rBigInt = r.isZero() ? ZERO : r.toBigInteger(signum);
2053         return new BigInteger[] {qBigInt, rBigInt};
2054     }
2055 
2056     /**
2057      * Returns a BigInteger whose value is <tt>(this<sup>exponent</sup>)</tt>.
2058      * Note that {@code exponent} is an integer rather than a BigInteger.
2059      *
2060      * @param  exponent exponent to which this BigInteger is to be raised.
2061      * @return <tt>this<sup>exponent</sup></tt>
2062      * @throws ArithmeticException {@code exponent} is negative.  (This would
2063      *         cause the operation to yield a non-integer value.)
2064      */
2065     public BigInteger pow(int exponent) {
2066         if (exponent < 0)
2067             throw new ArithmeticException("Negative exponent");
2068         if (signum==0)
2069             return (exponent==0 ? ONE : this);
2070 
2071         BigInteger partToSquare = this.abs();
2072 
2073         // Factor out powers of two from the base, as the exponentiation of
2074         // these can be done by left shifts only.
2075         // The remaining part can then be exponentiated faster.  The
2076         // powers of two will be multiplied back at the end.
2077         int powersOfTwo = partToSquare.getLowestSetBit();
2078 
2079         int remainingBits;
2080 
2081         // Factor the powers of two out quickly by shifting right, if needed.
2082         if (powersOfTwo > 0)
2083         {
2084             partToSquare = partToSquare.shiftRight(powersOfTwo);
2085             remainingBits = partToSquare.bitLength();
2086             if (remainingBits == 1)  // Nothing left but +/- 1?
2087                 if (signum<0 && (exponent&1)==1)
2088                     return NEGATIVE_ONE.shiftLeft(powersOfTwo*exponent);
2089                 else
2090                     return ONE.shiftLeft(powersOfTwo*exponent);
2091         }
2092         else
2093         {
2094             remainingBits = partToSquare.bitLength();
2095             if (remainingBits == 1)  // Nothing left but +/- 1?
2096                 if (signum<0 && (exponent&1)==1)
2097                     return NEGATIVE_ONE;
2098                 else
2099                     return ONE;
2100         }
2101 
2102         // This is a quick way to approximate the size of the result,
2103         // similar to doing log2[n] * exponent.  This will give an upper bound
2104         // of how big the result can be, and which algorithm to use.
2105         int scaleFactor = remainingBits * exponent;
2106 
2107         // Use slightly different algorithms for small and large operands.
2108         // See if the result will safely fit into a long. (Largest 2^63-1)
2109         if (partToSquare.mag.length==1 && scaleFactor <= 62)
2110         {
2111             // Small number algorithm.  Everything fits into a long.
2112             int newSign = (signum<0 && (exponent&1)==1 ? -1 : 1);
2113             long result = 1;
2114             long baseToPow2 = partToSquare.mag[0] & LONG_MASK;
2115 
2116             int workingExponent = exponent;
2117 
2118             // Perform exponentiation using repeated squaring trick
2119             while (workingExponent != 0) {
2120                 if ((workingExponent & 1)==1)
2121                     result = result * baseToPow2;
2122 
2123                 if ((workingExponent >>>= 1) != 0)
2124                     baseToPow2 = baseToPow2 * baseToPow2;
2125             }
2126 
2127             // Multiply back the powers of two (quickly, by shifting left)
2128             if (powersOfTwo > 0)
2129             {
2130                 int bitsToShift = powersOfTwo*exponent;
2131                 if (bitsToShift + scaleFactor <= 62) // Fits in long?
2132                     return valueOf((result << bitsToShift) * newSign);
2133                 else
2134                     return valueOf(result*newSign).shiftLeft(bitsToShift);
2135             }
2136             else
2137                 return valueOf(result*newSign);
2138         }
2139         else
2140         {
2141             // Large number algorithm.  This is basically identical to
2142             // the algorithm above, but calls multiply() and square()
2143             // which may use more efficient algorithms for large numbers.
2144             BigInteger answer = ONE;
2145 
2146             int workingExponent = exponent;
2147             // Perform exponentiation using repeated squaring trick
2148             while (workingExponent != 0) {
2149                 if ((workingExponent & 1)==1)
2150                     answer = answer.multiply(partToSquare);
2151 
2152                 if ((workingExponent >>>= 1) != 0)
2153                     partToSquare = partToSquare.square();
2154             }
2155             // Multiply back the (exponentiated) powers of two (quickly,
2156             // by shifting left)
2157             if (powersOfTwo > 0)
2158                 answer = answer.shiftLeft(powersOfTwo*exponent);
2159 
2160             if (signum<0 && (exponent&1)==1)
2161                 return answer.negate();
2162             else
2163                 return answer;
2164         }
2165     }
2166 
2167     /**
2168      * Returns a BigInteger whose value is the greatest common divisor of
2169      * {@code abs(this)} and {@code abs(val)}.  Returns 0 if
2170      * {@code this==0 && val==0}.
2171      *
2172      * @param  val value with which the GCD is to be computed.
2173      * @return {@code GCD(abs(this), abs(val))}
2174      */
2175     public BigInteger gcd(BigInteger val) {
2176         if (val.signum == 0)
2177             return this.abs();
2178         else if (this.signum == 0)
2179             return val.abs();
2180 
2181         MutableBigInteger a = new MutableBigInteger(this);
2182         MutableBigInteger b = new MutableBigInteger(val);
2183 
2184         MutableBigInteger result = a.hybridGCD(b);
2185 
2186         return result.toBigInteger(1);
2187     }
2188 
2189     /**
2190      * Package private method to return bit length for an integer.
2191      */
2192     static int bitLengthForInt(int n) {
2193         return 32 - Integer.numberOfLeadingZeros(n);
2194     }
2195 
2196     /**
2197      * Left shift int array a up to len by n bits. Returns the array that
2198      * results from the shift since space may have to be reallocated.
2199      */
2200     private static int[] leftShift(int[] a, int len, int n) {
2201         int nInts = n >>> 5;
2202         int nBits = n&0x1F;
2203         int bitsInHighWord = bitLengthForInt(a[0]);
2204 
2205         // If shift can be done without recopy, do so
2206         if (n <= (32-bitsInHighWord)) {
2207             primitiveLeftShift(a, len, nBits);
2208             return a;
2209         } else { // Array must be resized
2210             if (nBits <= (32-bitsInHighWord)) {
2211                 int result[] = new int[nInts+len];
2212                 System.arraycopy(a, 0, result, 0, len);
2213                 primitiveLeftShift(result, result.length, nBits);
2214                 return result;
2215             } else {
2216                 int result[] = new int[nInts+len+1];
2217                 System.arraycopy(a, 0, result, 0, len);
2218                 primitiveRightShift(result, result.length, 32 - nBits);
2219                 return result;
2220             }
2221         }
2222     }
2223 
2224     // shifts a up to len right n bits assumes no leading zeros, 0<n<32
2225     static void primitiveRightShift(int[] a, int len, int n) {
2226         int n2 = 32 - n;
2227         for (int i=len-1, c=a[i]; i>0; i--) {
2228             int b = c;
2229             c = a[i-1];
2230             a[i] = (c << n2) | (b >>> n);
2231         }
2232         a[0] >>>= n;
2233     }
2234 
2235     // shifts a up to len left n bits assumes no leading zeros, 0<=n<32
2236     static void primitiveLeftShift(int[] a, int len, int n) {
2237         if (len == 0 || n == 0)
2238             return;
2239 
2240         int n2 = 32 - n;
2241         for (int i=0, c=a[i], m=i+len-1; i<m; i++) {
2242             int b = c;
2243             c = a[i+1];
2244             a[i] = (b << n) | (c >>> n2);
2245         }
2246         a[len-1] <<= n;
2247     }
2248 
2249     /**
2250      * Calculate bitlength of contents of the first len elements an int array,
2251      * assuming there are no leading zero ints.
2252      */
2253     private static int bitLength(int[] val, int len) {
2254         if (len == 0)
2255             return 0;
2256         return ((len - 1) << 5) + bitLengthForInt(val[0]);
2257     }
2258 
2259     /**
2260      * Returns a BigInteger whose value is the absolute value of this
2261      * BigInteger.
2262      *
2263      * @return {@code abs(this)}
2264      */
2265     public BigInteger abs() {
2266         return (signum >= 0 ? this : this.negate());
2267     }
2268 
2269     /**
2270      * Returns a BigInteger whose value is {@code (-this)}.
2271      *
2272      * @return {@code -this}
2273      */
2274     public BigInteger negate() {
2275         return new BigInteger(this.mag, -this.signum);
2276     }
2277 
2278     /**
2279      * Returns the signum function of this BigInteger.
2280      *
2281      * @return -1, 0 or 1 as the value of this BigInteger is negative, zero or
2282      *         positive.
2283      */
2284     public int signum() {
2285         return this.signum;
2286     }
2287 
2288     // Modular Arithmetic Operations
2289 
2290     /**
2291      * Returns a BigInteger whose value is {@code (this mod m}).  This method
2292      * differs from {@code remainder} in that it always returns a
2293      * <i>non-negative</i> BigInteger.
2294      *
2295      * @param  m the modulus.
2296      * @return {@code this mod m}
2297      * @throws ArithmeticException {@code m} &le; 0
2298      * @see    #remainder
2299      */
2300     public BigInteger mod(BigInteger m) {
2301         if (m.signum <= 0)
2302             throw new ArithmeticException("BigInteger: modulus not positive");
2303 
2304         BigInteger result = this.remainder(m);
2305         return (result.signum >= 0 ? result : result.add(m));
2306     }
2307 
2308     /**
2309      * Returns a BigInteger whose value is
2310      * <tt>(this<sup>exponent</sup> mod m)</tt>.  (Unlike {@code pow}, this
2311      * method permits negative exponents.)
2312      *
2313      * @param  exponent the exponent.
2314      * @param  m the modulus.
2315      * @return <tt>this<sup>exponent</sup> mod m</tt>
2316      * @throws ArithmeticException {@code m} &le; 0 or the exponent is
2317      *         negative and this BigInteger is not <i>relatively
2318      *         prime</i> to {@code m}.
2319      * @see    #modInverse
2320      */
2321     public BigInteger modPow(BigInteger exponent, BigInteger m) {
2322         if (m.signum <= 0)
2323             throw new ArithmeticException("BigInteger: modulus not positive");
2324 
2325         // Trivial cases
2326         if (exponent.signum == 0)
2327             return (m.equals(ONE) ? ZERO : ONE);
2328 
2329         if (this.equals(ONE))
2330             return (m.equals(ONE) ? ZERO : ONE);
2331 
2332         if (this.equals(ZERO) && exponent.signum >= 0)
2333             return ZERO;
2334 
2335         if (this.equals(negConst[1]) && (!exponent.testBit(0)))
2336             return (m.equals(ONE) ? ZERO : ONE);
2337 
2338         boolean invertResult;
2339         if ((invertResult = (exponent.signum < 0)))
2340             exponent = exponent.negate();
2341 
2342         BigInteger base = (this.signum < 0 || this.compareTo(m) >= 0
2343                            ? this.mod(m) : this);
2344         BigInteger result;
2345         if (m.testBit(0)) { // odd modulus
2346             result = base.oddModPow(exponent, m);
2347         } else {
2348             /*
2349              * Even modulus.  Tear it into an "odd part" (m1) and power of two
2350              * (m2), exponentiate mod m1, manually exponentiate mod m2, and
2351              * use Chinese Remainder Theorem to combine results.
2352              */
2353 
2354             // Tear m apart into odd part (m1) and power of 2 (m2)
2355             int p = m.getLowestSetBit();   // Max pow of 2 that divides m
2356 
2357             BigInteger m1 = m.shiftRight(p);  // m/2**p
2358             BigInteger m2 = ONE.shiftLeft(p); // 2**p
2359 
2360             // Calculate new base from m1
2361             BigInteger base2 = (this.signum < 0 || this.compareTo(m1) >= 0
2362                                 ? this.mod(m1) : this);
2363 
2364             // Caculate (base ** exponent) mod m1.
2365             BigInteger a1 = (m1.equals(ONE) ? ZERO :
2366                              base2.oddModPow(exponent, m1));
2367 
2368             // Calculate (this ** exponent) mod m2
2369             BigInteger a2 = base.modPow2(exponent, p);
2370 
2371             // Combine results using Chinese Remainder Theorem
2372             BigInteger y1 = m2.modInverse(m1);
2373             BigInteger y2 = m1.modInverse(m2);
2374 
2375             result = a1.multiply(m2).multiply(y1).add
2376                      (a2.multiply(m1).multiply(y2)).mod(m);
2377         }
2378 
2379         return (invertResult ? result.modInverse(m) : result);
2380     }
2381 
2382     static int[] bnExpModThreshTable = {7, 25, 81, 241, 673, 1793,
2383                                                 Integer.MAX_VALUE}; // Sentinel
2384 
2385     /**
2386      * Returns a BigInteger whose value is x to the power of y mod z.
2387      * Assumes: z is odd && x < z.
2388      */
2389     private BigInteger oddModPow(BigInteger y, BigInteger z) {
2390     /*
2391      * The algorithm is adapted from Colin Plumb's C library.
2392      *
2393      * The window algorithm:
2394      * The idea is to keep a running product of b1 = n^(high-order bits of exp)
2395      * and then keep appending exponent bits to it.  The following patterns
2396      * apply to a 3-bit window (k = 3):
2397      * To append   0: square
2398      * To append   1: square, multiply by n^1
2399      * To append  10: square, multiply by n^1, square
2400      * To append  11: square, square, multiply by n^3
2401      * To append 100: square, multiply by n^1, square, square
2402      * To append 101: square, square, square, multiply by n^5
2403      * To append 110: square, square, multiply by n^3, square
2404      * To append 111: square, square, square, multiply by n^7
2405      *
2406      * Since each pattern involves only one multiply, the longer the pattern
2407      * the better, except that a 0 (no multiplies) can be appended directly.
2408      * We precompute a table of odd powers of n, up to 2^k, and can then
2409      * multiply k bits of exponent at a time.  Actually, assuming random
2410      * exponents, there is on average one zero bit between needs to
2411      * multiply (1/2 of the time there's none, 1/4 of the time there's 1,
2412      * 1/8 of the time, there's 2, 1/32 of the time, there's 3, etc.), so
2413      * you have to do one multiply per k+1 bits of exponent.
2414      *
2415      * The loop walks down the exponent, squaring the result buffer as
2416      * it goes.  There is a wbits+1 bit lookahead buffer, buf, that is
2417      * filled with the upcoming exponent bits.  (What is read after the
2418      * end of the exponent is unimportant, but it is filled with zero here.)
2419      * When the most-significant bit of this buffer becomes set, i.e.
2420      * (buf & tblmask) != 0, we have to decide what pattern to multiply
2421      * by, and when to do it.  We decide, remember to do it in future
2422      * after a suitable number of squarings have passed (e.g. a pattern
2423      * of "100" in the buffer requires that we multiply by n^1 immediately;
2424      * a pattern of "110" calls for multiplying by n^3 after one more
2425      * squaring), clear the buffer, and continue.
2426      *
2427      * When we start, there is one more optimization: the result buffer
2428      * is implcitly one, so squaring it or multiplying by it can be
2429      * optimized away.  Further, if we start with a pattern like "100"
2430      * in the lookahead window, rather than placing n into the buffer
2431      * and then starting to square it, we have already computed n^2
2432      * to compute the odd-powers table, so we can place that into
2433      * the buffer and save a squaring.
2434      *
2435      * This means that if you have a k-bit window, to compute n^z,
2436      * where z is the high k bits of the exponent, 1/2 of the time
2437      * it requires no squarings.  1/4 of the time, it requires 1
2438      * squaring, ... 1/2^(k-1) of the time, it reqires k-2 squarings.
2439      * And the remaining 1/2^(k-1) of the time, the top k bits are a
2440      * 1 followed by k-1 0 bits, so it again only requires k-2
2441      * squarings, not k-1.  The average of these is 1.  Add that
2442      * to the one squaring we have to do to compute the table,
2443      * and you'll see that a k-bit window saves k-2 squarings
2444      * as well as reducing the multiplies.  (It actually doesn't
2445      * hurt in the case k = 1, either.)
2446      */
2447         // Special case for exponent of one
2448         if (y.equals(ONE))
2449             return this;
2450 
2451         // Special case for base of zero
2452         if (signum==0)
2453             return ZERO;
2454 
2455         int[] base = mag.clone();
2456         int[] exp = y.mag;
2457         int[] mod = z.mag;
2458         int modLen = mod.length;
2459 
2460         // Select an appropriate window size
2461         int wbits = 0;
2462         int ebits = bitLength(exp, exp.length);
2463         // if exponent is 65537 (0x10001), use minimum window size
2464         if ((ebits != 17) || (exp[0] != 65537)) {
2465             while (ebits > bnExpModThreshTable[wbits]) {
2466                 wbits++;
2467             }
2468         }
2469 
2470         // Calculate appropriate table size
2471         int tblmask = 1 << wbits;
2472 
2473         // Allocate table for precomputed odd powers of base in Montgomery form
2474         int[][] table = new int[tblmask][];
2475         for (int i=0; i<tblmask; i++)
2476             table[i] = new int[modLen];
2477 
2478         // Compute the modular inverse
2479         int inv = -MutableBigInteger.inverseMod32(mod[modLen-1]);
2480 
2481         // Convert base to Montgomery form
2482         int[] a = leftShift(base, base.length, modLen << 5);
2483 
2484         MutableBigInteger q = new MutableBigInteger(),
2485                           a2 = new MutableBigInteger(a),
2486                           b2 = new MutableBigInteger(mod);
2487 
2488         MutableBigInteger r= a2.divide(b2, q);
2489         table[0] = r.toIntArray();
2490 
2491         // Pad table[0] with leading zeros so its length is at least modLen
2492         if (table[0].length < modLen) {
2493            int offset = modLen - table[0].length;
2494            int[] t2 = new int[modLen];
2495            for (int i=0; i<table[0].length; i++)
2496                t2[i+offset] = table[0][i];
2497            table[0] = t2;
2498         }
2499 
2500         // Set b to the square of the base
2501         int[] b = squareToLen(table[0], modLen, null);
2502         b = montReduce(b, mod, modLen, inv);
2503 
2504         // Set t to high half of b
2505         int[] t = Arrays.copyOf(b, modLen);
2506 
2507         // Fill in the table with odd powers of the base
2508         for (int i=1; i<tblmask; i++) {
2509             int[] prod = multiplyToLen(t, modLen, table[i-1], modLen, null);
2510             table[i] = montReduce(prod, mod, modLen, inv);
2511         }
2512 
2513         // Pre load the window that slides over the exponent
2514         int bitpos = 1 << ((ebits-1) & (32-1));
2515 
2516         int buf = 0;
2517         int elen = exp.length;
2518         int eIndex = 0;
2519         for (int i = 0; i <= wbits; i++) {
2520             buf = (buf << 1) | (((exp[eIndex] & bitpos) != 0)?1:0);
2521             bitpos >>>= 1;
2522             if (bitpos == 0) {
2523                 eIndex++;
2524                 bitpos = 1 << (32-1);
2525                 elen--;
2526             }
2527         }
2528 
2529         int multpos = ebits;
2530 
2531         // The first iteration, which is hoisted out of the main loop
2532         ebits--;
2533         boolean isone = true;
2534 
2535         multpos = ebits - wbits;
2536         while ((buf & 1) == 0) {
2537             buf >>>= 1;
2538             multpos++;
2539         }
2540 
2541         int[] mult = table[buf >>> 1];
2542 
2543         buf = 0;
2544         if (multpos == ebits)
2545             isone = false;
2546 
2547         // The main loop
2548         while(true) {
2549             ebits--;
2550             // Advance the window
2551             buf <<= 1;
2552 
2553             if (elen != 0) {
2554                 buf |= ((exp[eIndex] & bitpos) != 0) ? 1 : 0;
2555                 bitpos >>>= 1;
2556                 if (bitpos == 0) {
2557                     eIndex++;
2558                     bitpos = 1 << (32-1);
2559                     elen--;
2560                 }
2561             }
2562 
2563             // Examine the window for pending multiplies
2564             if ((buf & tblmask) != 0) {
2565                 multpos = ebits - wbits;
2566                 while ((buf & 1) == 0) {
2567                     buf >>>= 1;
2568                     multpos++;
2569                 }
2570                 mult = table[buf >>> 1];
2571                 buf = 0;
2572             }
2573 
2574             // Perform multiply
2575             if (ebits == multpos) {
2576                 if (isone) {
2577                     b = mult.clone();
2578                     isone = false;
2579                 } else {
2580                     t = b;
2581                     a = multiplyToLen(t, modLen, mult, modLen, a);
2582                     a = montReduce(a, mod, modLen, inv);
2583                     t = a; a = b; b = t;
2584                 }
2585             }
2586 
2587             // Check if done
2588             if (ebits == 0)
2589                 break;
2590 
2591             // Square the input
2592             if (!isone) {
2593                 t = b;
2594                 a = squareToLen(t, modLen, a);
2595                 a = montReduce(a, mod, modLen, inv);
2596                 t = a; a = b; b = t;
2597             }
2598         }
2599 
2600         // Convert result out of Montgomery form and return
2601         int[] t2 = new int[2*modLen];
2602         System.arraycopy(b, 0, t2, modLen, modLen);
2603 
2604         b = montReduce(t2, mod, modLen, inv);
2605 
2606         t2 = Arrays.copyOf(b, modLen);
2607 
2608         return new BigInteger(1, t2);
2609     }
2610 
2611     /**
2612      * Montgomery reduce n, modulo mod.  This reduces modulo mod and divides
2613      * by 2^(32*mlen). Adapted from Colin Plumb's C library.
2614      */
2615     private static int[] montReduce(int[] n, int[] mod, int mlen, int inv) {
2616         int c=0;
2617         int len = mlen;
2618         int offset=0;
2619 
2620         do {
2621             int nEnd = n[n.length-1-offset];
2622             int carry = mulAdd(n, mod, offset, mlen, inv * nEnd);
2623             c += addOne(n, offset, mlen, carry);
2624             offset++;
2625         } while(--len > 0);
2626 
2627         while(c>0)
2628             c += subN(n, mod, mlen);
2629 
2630         while (intArrayCmpToLen(n, mod, mlen) >= 0)
2631             subN(n, mod, mlen);
2632 
2633         return n;
2634     }
2635 
2636 
2637     /*
2638      * Returns -1, 0 or +1 as big-endian unsigned int array arg1 is less than,
2639      * equal to, or greater than arg2 up to length len.
2640      */
2641     private static int intArrayCmpToLen(int[] arg1, int[] arg2, int len) {
2642         for (int i=0; i<len; i++) {
2643             long b1 = arg1[i] & LONG_MASK;
2644             long b2 = arg2[i] & LONG_MASK;
2645             if (b1 < b2)
2646                 return -1;
2647             if (b1 > b2)
2648                 return 1;
2649         }
2650         return 0;
2651     }
2652 
2653     /**
2654      * Subtracts two numbers of same length, returning borrow.
2655      */
2656     private static int subN(int[] a, int[] b, int len) {
2657         long sum = 0;
2658 
2659         while(--len >= 0) {
2660             sum = (a[len] & LONG_MASK) -
2661                  (b[len] & LONG_MASK) + (sum >> 32);
2662             a[len] = (int)sum;
2663         }
2664 
2665         return (int)(sum >> 32);
2666     }
2667 
2668     /**
2669      * Multiply an array by one word k and add to result, return the carry
2670      */
2671     static int mulAdd(int[] out, int[] in, int offset, int len, int k) {
2672         long kLong = k & LONG_MASK;
2673         long carry = 0;
2674 
2675         offset = out.length-offset - 1;
2676         for (int j=len-1; j >= 0; j--) {
2677             long product = (in[j] & LONG_MASK) * kLong +
2678                            (out[offset] & LONG_MASK) + carry;
2679             out[offset--] = (int)product;
2680             carry = product >>> 32;
2681         }
2682         return (int)carry;
2683     }
2684 
2685     /**
2686      * Add one word to the number a mlen words into a. Return the resulting
2687      * carry.
2688      */
2689     static int addOne(int[] a, int offset, int mlen, int carry) {
2690         offset = a.length-1-mlen-offset;
2691         long t = (a[offset] & LONG_MASK) + (carry & LONG_MASK);
2692 
2693         a[offset] = (int)t;
2694         if ((t >>> 32) == 0)
2695             return 0;
2696         while (--mlen >= 0) {
2697             if (--offset < 0) { // Carry out of number
2698                 return 1;
2699             } else {
2700                 a[offset]++;
2701                 if (a[offset] != 0)
2702                     return 0;
2703             }
2704         }
2705         return 1;
2706     }
2707 
2708     /**
2709      * Returns a BigInteger whose value is (this ** exponent) mod (2**p)
2710      */
2711     private BigInteger modPow2(BigInteger exponent, int p) {
2712         /*
2713          * Perform exponentiation using repeated squaring trick, chopping off
2714          * high order bits as indicated by modulus.
2715          */
2716         BigInteger result = ONE;
2717         BigInteger baseToPow2 = this.mod2(p);
2718         int expOffset = 0;
2719 
2720         int limit = exponent.bitLength();
2721 
2722         if (this.testBit(0))
2723            limit = (p-1) < limit ? (p-1) : limit;
2724 
2725         while (expOffset < limit) {
2726             if (exponent.testBit(expOffset))
2727                 result = result.multiply(baseToPow2).mod2(p);
2728             expOffset++;
2729             if (expOffset < limit)
2730                 baseToPow2 = baseToPow2.square().mod2(p);
2731         }
2732 
2733         return result;
2734     }
2735 
2736     /**
2737      * Returns a BigInteger whose value is this mod(2**p).
2738      * Assumes that this {@code BigInteger >= 0} and {@code p > 0}.
2739      */
2740     private BigInteger mod2(int p) {
2741         if (bitLength() <= p)
2742             return this;
2743 
2744         // Copy remaining ints of mag
2745         int numInts = (p + 31) >>> 5;
2746         int[] mag = new int[numInts];
2747         System.arraycopy(this.mag, (this.mag.length - numInts), mag, 0, numInts);
2748 
2749         // Mask out any excess bits
2750         int excessBits = (numInts << 5) - p;
2751         mag[0] &= (1L << (32-excessBits)) - 1;
2752 
2753         return (mag[0]==0 ? new BigInteger(1, mag) : new BigInteger(mag, 1));
2754     }
2755 
2756     /**
2757      * Returns a BigInteger whose value is {@code (this}<sup>-1</sup> {@code mod m)}.
2758      *
2759      * @param  m the modulus.
2760      * @return {@code this}<sup>-1</sup> {@code mod m}.
2761      * @throws ArithmeticException {@code  m} &le; 0, or this BigInteger
2762      *         has no multiplicative inverse mod m (that is, this BigInteger
2763      *         is not <i>relatively prime</i> to m).
2764      */
2765     public BigInteger modInverse(BigInteger m) {
2766         if (m.signum != 1)
2767             throw new ArithmeticException("BigInteger: modulus not positive");
2768 
2769         if (m.equals(ONE))
2770             return ZERO;
2771 
2772         // Calculate (this mod m)
2773         BigInteger modVal = this;
2774         if (signum < 0 || (this.compareMagnitude(m) >= 0))
2775             modVal = this.mod(m);
2776 
2777         if (modVal.equals(ONE))
2778             return ONE;
2779 
2780         MutableBigInteger a = new MutableBigInteger(modVal);
2781         MutableBigInteger b = new MutableBigInteger(m);
2782 
2783         MutableBigInteger result = a.mutableModInverse(b);
2784         return result.toBigInteger(1);
2785     }
2786 
2787     // Shift Operations
2788 
2789     /**
2790      * Returns a BigInteger whose value is {@code (this << n)}.
2791      * The shift distance, {@code n}, may be negative, in which case
2792      * this method performs a right shift.
2793      * (Computes <tt>floor(this * 2<sup>n</sup>)</tt>.)
2794      *
2795      * @param  n shift distance, in bits.
2796      * @return {@code this << n}
2797      * @throws ArithmeticException if the shift distance is {@code
2798      *         Integer.MIN_VALUE}.
2799      * @see #shiftRight
2800      */
2801     public BigInteger shiftLeft(int n) {
2802         if (signum == 0)
2803             return ZERO;
2804         if (n==0)
2805             return this;
2806         if (n<0) {
2807             if (n == Integer.MIN_VALUE) {
2808                 throw new ArithmeticException("Shift distance of Integer.MIN_VALUE not supported.");
2809             } else {
2810                 return shiftRight(-n);
2811             }
2812         }
2813         int[] newMag = shiftLeft(mag, n);
2814 
2815         return new BigInteger(newMag, signum);
2816     }
2817 
2818     private static int[] shiftLeft(int[] mag, int n) {
2819         int nInts = n >>> 5;
2820         int nBits = n & 0x1f;
2821         int magLen = mag.length;
2822         int newMag[] = null;
2823 
2824         if (nBits == 0) {
2825             newMag = new int[magLen + nInts];
2826             System.arraycopy(mag, 0, newMag, 0, magLen);
2827         } else {
2828             int i = 0;
2829             int nBits2 = 32 - nBits;
2830             int highBits = mag[0] >>> nBits2;
2831             if (highBits != 0) {
2832                 newMag = new int[magLen + nInts + 1];
2833                 newMag[i++] = highBits;
2834             } else {
2835                 newMag = new int[magLen + nInts];
2836             }
2837             int j=0;
2838             while (j < magLen-1)
2839                 newMag[i++] = mag[j++] << nBits | mag[j] >>> nBits2;
2840             newMag[i] = mag[j] << nBits;
2841         }
2842         return newMag;
2843     }
2844 
2845     /**
2846      * Returns a BigInteger whose value is {@code (this >> n)}.  Sign
2847      * extension is performed.  The shift distance, {@code n}, may be
2848      * negative, in which case this method performs a left shift.
2849      * (Computes <tt>floor(this / 2<sup>n</sup>)</tt>.)
2850      *
2851      * @param  n shift distance, in bits.
2852      * @return {@code this >> n}
2853      * @throws ArithmeticException if the shift distance is {@code
2854      *         Integer.MIN_VALUE}.
2855      * @see #shiftLeft
2856      */
2857     public BigInteger shiftRight(int n) {
2858         if (n==0)
2859             return this;
2860         if (n<0) {
2861             if (n == Integer.MIN_VALUE) {
2862                 throw new ArithmeticException("Shift distance of Integer.MIN_VALUE not supported.");
2863             } else {
2864                 return shiftLeft(-n);
2865             }
2866         }
2867 
2868         int nInts = n >>> 5;
2869         int nBits = n & 0x1f;
2870         int magLen = mag.length;
2871         int newMag[] = null;
2872 
2873         // Special case: entire contents shifted off the end
2874         if (nInts >= magLen)
2875             return (signum >= 0 ? ZERO : negConst[1]);
2876 
2877         if (nBits == 0) {
2878             int newMagLen = magLen - nInts;
2879             newMag = Arrays.copyOf(mag, newMagLen);
2880         } else {
2881             int i = 0;
2882             int highBits = mag[0] >>> nBits;
2883             if (highBits != 0) {
2884                 newMag = new int[magLen - nInts];
2885                 newMag[i++] = highBits;
2886             } else {
2887                 newMag = new int[magLen - nInts -1];
2888             }
2889 
2890             int nBits2 = 32 - nBits;
2891             int j=0;
2892             while (j < magLen - nInts - 1)
2893                 newMag[i++] = (mag[j++] << nBits2) | (mag[j] >>> nBits);
2894         }
2895 
2896         if (signum < 0) {
2897             // Find out whether any one-bits were shifted off the end.
2898             boolean onesLost = false;
2899             for (int i=magLen-1, j=magLen-nInts; i>=j && !onesLost; i--)
2900                 onesLost = (mag[i] != 0);
2901             if (!onesLost && nBits != 0)
2902                 onesLost = (mag[magLen - nInts - 1] << (32 - nBits) != 0);
2903 
2904             if (onesLost)
2905                 newMag = javaIncrement(newMag);
2906         }
2907 
2908         return new BigInteger(newMag, signum);
2909     }
2910 
2911     int[] javaIncrement(int[] val) {
2912         int lastSum = 0;
2913         for (int i=val.length-1;  i >= 0 && lastSum == 0; i--)
2914             lastSum = (val[i] += 1);
2915         if (lastSum == 0) {
2916             val = new int[val.length+1];
2917             val[0] = 1;
2918         }
2919         return val;
2920     }
2921 
2922     // Bitwise Operations
2923 
2924     /**
2925      * Returns a BigInteger whose value is {@code (this & val)}.  (This
2926      * method returns a negative BigInteger if and only if this and val are
2927      * both negative.)
2928      *
2929      * @param val value to be AND'ed with this BigInteger.
2930      * @return {@code this & val}
2931      */
2932     public BigInteger and(BigInteger val) {
2933         int[] result = new int[Math.max(intLength(), val.intLength())];
2934         for (int i=0; i<result.length; i++)
2935             result[i] = (getInt(result.length-i-1)
2936                          & val.getInt(result.length-i-1));
2937 
2938         return valueOf(result);
2939     }
2940 
2941     /**
2942      * Returns a BigInteger whose value is {@code (this | val)}.  (This method
2943      * returns a negative BigInteger if and only if either this or val is
2944      * negative.)
2945      *
2946      * @param val value to be OR'ed with this BigInteger.
2947      * @return {@code this | val}
2948      */
2949     public BigInteger or(BigInteger val) {
2950         int[] result = new int[Math.max(intLength(), val.intLength())];
2951         for (int i=0; i<result.length; i++)
2952             result[i] = (getInt(result.length-i-1)
2953                          | val.getInt(result.length-i-1));
2954 
2955         return valueOf(result);
2956     }
2957 
2958     /**
2959      * Returns a BigInteger whose value is {@code (this ^ val)}.  (This method
2960      * returns a negative BigInteger if and only if exactly one of this and
2961      * val are negative.)
2962      *
2963      * @param val value to be XOR'ed with this BigInteger.
2964      * @return {@code this ^ val}
2965      */
2966     public BigInteger xor(BigInteger val) {
2967         int[] result = new int[Math.max(intLength(), val.intLength())];
2968         for (int i=0; i<result.length; i++)
2969             result[i] = (getInt(result.length-i-1)
2970                          ^ val.getInt(result.length-i-1));
2971 
2972         return valueOf(result);
2973     }
2974 
2975     /**
2976      * Returns a BigInteger whose value is {@code (~this)}.  (This method
2977      * returns a negative value if and only if this BigInteger is
2978      * non-negative.)
2979      *
2980      * @return {@code ~this}
2981      */
2982     public BigInteger not() {
2983         int[] result = new int[intLength()];
2984         for (int i=0; i<result.length; i++)
2985             result[i] = ~getInt(result.length-i-1);
2986 
2987         return valueOf(result);
2988     }
2989 
2990     /**
2991      * Returns a BigInteger whose value is {@code (this & ~val)}.  This
2992      * method, which is equivalent to {@code and(val.not())}, is provided as
2993      * a convenience for masking operations.  (This method returns a negative
2994      * BigInteger if and only if {@code this} is negative and {@code val} is
2995      * positive.)
2996      *
2997      * @param val value to be complemented and AND'ed with this BigInteger.
2998      * @return {@code this & ~val}
2999      */
3000     public BigInteger andNot(BigInteger val) {
3001         int[] result = new int[Math.max(intLength(), val.intLength())];
3002         for (int i=0; i<result.length; i++)
3003             result[i] = (getInt(result.length-i-1)
3004                          & ~val.getInt(result.length-i-1));
3005 
3006         return valueOf(result);
3007     }
3008 
3009 
3010     // Single Bit Operations
3011 
3012     /**
3013      * Returns {@code true} if and only if the designated bit is set.
3014      * (Computes {@code ((this & (1<<n)) != 0)}.)
3015      *
3016      * @param  n index of bit to test.
3017      * @return {@code true} if and only if the designated bit is set.
3018      * @throws ArithmeticException {@code n} is negative.
3019      */
3020     public boolean testBit(int n) {
3021         if (n<0)
3022             throw new ArithmeticException("Negative bit address");
3023 
3024         return (getInt(n >>> 5) & (1 << (n & 31))) != 0;
3025     }
3026 
3027     /**
3028      * Returns a BigInteger whose value is equivalent to this BigInteger
3029      * with the designated bit set.  (Computes {@code (this | (1<<n))}.)
3030      *
3031      * @param  n index of bit to set.
3032      * @return {@code this | (1<<n)}
3033      * @throws ArithmeticException {@code n} is negative.
3034      */
3035     public BigInteger setBit(int n) {
3036         if (n<0)
3037             throw new ArithmeticException("Negative bit address");
3038 
3039         int intNum = n >>> 5;
3040         int[] result = new int[Math.max(intLength(), intNum+2)];
3041 
3042         for (int i=0; i<result.length; i++)
3043             result[result.length-i-1] = getInt(i);
3044 
3045         result[result.length-intNum-1] |= (1 << (n & 31));
3046 
3047         return valueOf(result);
3048     }
3049 
3050     /**
3051      * Returns a BigInteger whose value is equivalent to this BigInteger
3052      * with the designated bit cleared.
3053      * (Computes {@code (this & ~(1<<n))}.)
3054      *
3055      * @param  n index of bit to clear.
3056      * @return {@code this & ~(1<<n)}
3057      * @throws ArithmeticException {@code n} is negative.
3058      */
3059     public BigInteger clearBit(int n) {
3060         if (n<0)
3061             throw new ArithmeticException("Negative bit address");
3062 
3063         int intNum = n >>> 5;
3064         int[] result = new int[Math.max(intLength(), ((n + 1) >>> 5) + 1)];
3065 
3066         for (int i=0; i<result.length; i++)
3067             result[result.length-i-1] = getInt(i);
3068 
3069         result[result.length-intNum-1] &= ~(1 << (n & 31));
3070 
3071         return valueOf(result);
3072     }
3073 
3074     /**
3075      * Returns a BigInteger whose value is equivalent to this BigInteger
3076      * with the designated bit flipped.
3077      * (Computes {@code (this ^ (1<<n))}.)
3078      *
3079      * @param  n index of bit to flip.
3080      * @return {@code this ^ (1<<n)}
3081      * @throws ArithmeticException {@code n} is negative.
3082      */
3083     public BigInteger flipBit(int n) {
3084         if (n<0)
3085             throw new ArithmeticException("Negative bit address");
3086 
3087         int intNum = n >>> 5;
3088         int[] result = new int[Math.max(intLength(), intNum+2)];
3089 
3090         for (int i=0; i<result.length; i++)
3091             result[result.length-i-1] = getInt(i);
3092 
3093         result[result.length-intNum-1] ^= (1 << (n & 31));
3094 
3095         return valueOf(result);
3096     }
3097 
3098     /**
3099      * Returns the index of the rightmost (lowest-order) one bit in this
3100      * BigInteger (the number of zero bits to the right of the rightmost
3101      * one bit).  Returns -1 if this BigInteger contains no one bits.
3102      * (Computes {@code (this==0? -1 : log2(this & -this))}.)
3103      *
3104      * @return index of the rightmost one bit in this BigInteger.
3105      */
3106     public int getLowestSetBit() {
3107         @SuppressWarnings("deprecation") int lsb = lowestSetBit - 2;
3108         if (lsb == -2) {  // lowestSetBit not initialized yet
3109             lsb = 0;
3110             if (signum == 0) {
3111                 lsb -= 1;
3112             } else {
3113                 // Search for lowest order nonzero int
3114                 int i,b;
3115                 for (i=0; (b = getInt(i))==0; i++)
3116                     ;
3117                 lsb += (i << 5) + Integer.numberOfTrailingZeros(b);
3118             }
3119             lowestSetBit = lsb + 2;
3120         }
3121         return lsb;
3122     }
3123 
3124 
3125     // Miscellaneous Bit Operations
3126 
3127     /**
3128      * Returns the number of bits in the minimal two's-complement
3129      * representation of this BigInteger, <i>excluding</i> a sign bit.
3130      * For positive BigIntegers, this is equivalent to the number of bits in
3131      * the ordinary binary representation.  (Computes
3132      * {@code (ceil(log2(this < 0 ? -this : this+1)))}.)
3133      *
3134      * @return number of bits in the minimal two's-complement
3135      *         representation of this BigInteger, <i>excluding</i> a sign bit.
3136      */
3137     public int bitLength() {
3138         @SuppressWarnings("deprecation") int n = bitLength - 1;
3139         if (n == -1) { // bitLength not initialized yet
3140             int[] m = mag;
3141             int len = m.length;
3142             if (len == 0) {
3143                 n = 0; // offset by one to initialize
3144             }  else {
3145                 // Calculate the bit length of the magnitude
3146                 int magBitLength = ((len - 1) << 5) + bitLengthForInt(mag[0]);
3147                  if (signum < 0) {
3148                      // Check if magnitude is a power of two
3149                      boolean pow2 = (Integer.bitCount(mag[0]) == 1);
3150                      for (int i=1; i< len && pow2; i++)
3151                          pow2 = (mag[i] == 0);
3152 
3153                      n = (pow2 ? magBitLength -1 : magBitLength);
3154                  } else {
3155                      n = magBitLength;
3156                  }
3157             }
3158             bitLength = n + 1;
3159         }
3160         return n;
3161     }
3162 
3163     /**
3164      * Returns the number of bits in the two's complement representation
3165      * of this BigInteger that differ from its sign bit.  This method is
3166      * useful when implementing bit-vector style sets atop BigIntegers.
3167      *
3168      * @return number of bits in the two's complement representation
3169      *         of this BigInteger that differ from its sign bit.
3170      */
3171     public int bitCount() {
3172         @SuppressWarnings("deprecation") int bc = bitCount - 1;
3173         if (bc == -1) {  // bitCount not initialized yet
3174             bc = 0;      // offset by one to initialize
3175             // Count the bits in the magnitude
3176             for (int i=0; i<mag.length; i++)
3177                 bc += Integer.bitCount(mag[i]);
3178             if (signum < 0) {
3179                 // Count the trailing zeros in the magnitude
3180                 int magTrailingZeroCount = 0, j;
3181                 for (j=mag.length-1; mag[j]==0; j--)
3182                     magTrailingZeroCount += 32;
3183                 magTrailingZeroCount += Integer.numberOfTrailingZeros(mag[j]);
3184                 bc += magTrailingZeroCount - 1;
3185             }
3186             bitCount = bc + 1;
3187         }
3188         return bc;
3189     }
3190 
3191     // Primality Testing
3192 
3193     /**
3194      * Returns {@code true} if this BigInteger is probably prime,
3195      * {@code false} if it's definitely composite.  If
3196      * {@code certainty} is &le; 0, {@code true} is
3197      * returned.
3198      *
3199      * @param  certainty a measure of the uncertainty that the caller is
3200      *         willing to tolerate: if the call returns {@code true}
3201      *         the probability that this BigInteger is prime exceeds
3202      *         (1 - 1/2<sup>{@code certainty}</sup>).  The execution time of
3203      *         this method is proportional to the value of this parameter.
3204      * @return {@code true} if this BigInteger is probably prime,
3205      *         {@code false} if it's definitely composite.
3206      */
3207     public boolean isProbablePrime(int certainty) {
3208         if (certainty <= 0)
3209             return true;
3210         BigInteger w = this.abs();
3211         if (w.equals(TWO))
3212             return true;
3213         if (!w.testBit(0) || w.equals(ONE))
3214             return false;
3215 
3216         return w.primeToCertainty(certainty, null);
3217     }
3218 
3219     // Comparison Operations
3220 
3221     /**
3222      * Compares this BigInteger with the specified BigInteger.  This
3223      * method is provided in preference to individual methods for each
3224      * of the six boolean comparison operators ({@literal <}, ==,
3225      * {@literal >}, {@literal >=}, !=, {@literal <=}).  The suggested
3226      * idiom for performing these comparisons is: {@code
3227      * (x.compareTo(y)} &lt;<i>op</i>&gt; {@code 0)}, where
3228      * &lt;<i>op</i>&gt; is one of the six comparison operators.
3229      *
3230      * @param  val BigInteger to which this BigInteger is to be compared.
3231      * @return -1, 0 or 1 as this BigInteger is numerically less than, equal
3232      *         to, or greater than {@code val}.
3233      */
3234     public int compareTo(BigInteger val) {
3235         if (signum == val.signum) {
3236             switch (signum) {
3237             case 1:
3238                 return compareMagnitude(val);
3239             case -1:
3240                 return val.compareMagnitude(this);
3241             default:
3242                 return 0;
3243             }
3244         }
3245         return signum > val.signum ? 1 : -1;
3246     }
3247 
3248     /**
3249      * Compares the magnitude array of this BigInteger with the specified
3250      * BigInteger's. This is the version of compareTo ignoring sign.
3251      *
3252      * @param val BigInteger whose magnitude array to be compared.
3253      * @return -1, 0 or 1 as this magnitude array is less than, equal to or
3254      *         greater than the magnitude aray for the specified BigInteger's.
3255      */
3256     final int compareMagnitude(BigInteger val) {
3257         int[] m1 = mag;
3258         int len1 = m1.length;
3259         int[] m2 = val.mag;
3260         int len2 = m2.length;
3261         if (len1 < len2)
3262             return -1;
3263         if (len1 > len2)
3264             return 1;
3265         for (int i = 0; i < len1; i++) {
3266             int a = m1[i];
3267             int b = m2[i];
3268             if (a != b)
3269                 return ((a & LONG_MASK) < (b & LONG_MASK)) ? -1 : 1;
3270         }
3271         return 0;
3272     }
3273 
3274     /**
3275      * Version of compareMagnitude that compares magnitude with long value.
3276      * val can't be Long.MIN_VALUE.
3277      */
3278     final int compareMagnitude(long val) {
3279         assert val != Long.MIN_VALUE;
3280         int[] m1 = mag;
3281         int len = m1.length;
3282         if(len > 2) {
3283             return 1;
3284         }
3285         if (val < 0) {
3286             val = -val;
3287         }
3288         int highWord = (int)(val >>> 32);
3289         if (highWord==0) {
3290             if (len < 1)
3291                 return -1;
3292             if (len > 1)
3293                 return 1;
3294             int a = m1[0];
3295             int b = (int)val;
3296             if (a != b) {
3297                 return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
3298             }
3299             return 0;
3300         } else {
3301             if (len < 2)
3302                 return -1;
3303             int a = m1[0];
3304             int b = highWord;
3305             if (a != b) {
3306                 return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
3307             }
3308             a = m1[1];
3309             b = (int)val;
3310             if (a != b) {
3311                 return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
3312             }
3313             return 0;
3314         }
3315     }
3316 
3317     /**
3318      * Compares this BigInteger with the specified Object for equality.
3319      *
3320      * @param  x Object to which this BigInteger is to be compared.
3321      * @return {@code true} if and only if the specified Object is a
3322      *         BigInteger whose value is numerically equal to this BigInteger.
3323      */
3324     public boolean equals(Object x) {
3325         // This test is just an optimization, which may or may not help
3326         if (x == this)
3327             return true;
3328 
3329         if (!(x instanceof BigInteger))
3330             return false;
3331 
3332         BigInteger xInt = (BigInteger) x;
3333         if (xInt.signum != signum)
3334             return false;
3335 
3336         int[] m = mag;
3337         int len = m.length;
3338         int[] xm = xInt.mag;
3339         if (len != xm.length)
3340             return false;
3341 
3342         for (int i = 0; i < len; i++)
3343             if (xm[i] != m[i])
3344                 return false;
3345 
3346         return true;
3347     }
3348 
3349     /**
3350      * Returns the minimum of this BigInteger and {@code val}.
3351      *
3352      * @param  val value with which the minimum is to be computed.
3353      * @return the BigInteger whose value is the lesser of this BigInteger and
3354      *         {@code val}.  If they are equal, either may be returned.
3355      */
3356     public BigInteger min(BigInteger val) {
3357         return (compareTo(val)<0 ? this : val);
3358     }
3359 
3360     /**
3361      * Returns the maximum of this BigInteger and {@code val}.
3362      *
3363      * @param  val value with which the maximum is to be computed.
3364      * @return the BigInteger whose value is the greater of this and
3365      *         {@code val}.  If they are equal, either may be returned.
3366      */
3367     public BigInteger max(BigInteger val) {
3368         return (compareTo(val)>0 ? this : val);
3369     }
3370 
3371 
3372     // Hash Function
3373 
3374     /**
3375      * Returns the hash code for this BigInteger.
3376      *
3377      * @return hash code for this BigInteger.
3378      */
3379     public int hashCode() {
3380         int hashCode = 0;
3381 
3382         for (int i=0; i<mag.length; i++)
3383             hashCode = (int)(31*hashCode + (mag[i] & LONG_MASK));
3384 
3385         return hashCode * signum;
3386     }
3387 
3388     /**
3389      * Returns the String representation of this BigInteger in the
3390      * given radix.  If the radix is outside the range from {@link
3391      * Character#MIN_RADIX} to {@link Character#MAX_RADIX} inclusive,
3392      * it will default to 10 (as is the case for
3393      * {@code Integer.toString}).  The digit-to-character mapping
3394      * provided by {@code Character.forDigit} is used, and a minus
3395      * sign is prepended if appropriate.  (This representation is
3396      * compatible with the {@link #BigInteger(String, int) (String,
3397      * int)} constructor.)
3398      *
3399      * @param  radix  radix of the String representation.
3400      * @return String representation of this BigInteger in the given radix.
3401      * @see    Integer#toString
3402      * @see    Character#forDigit
3403      * @see    #BigInteger(java.lang.String, int)
3404      */
3405     public String toString(int radix) {
3406         if (signum == 0)
3407             return "0";
3408         if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
3409             radix = 10;
3410 
3411         // If it's small enough, use smallToString.
3412         if (mag.length <= SCHOENHAGE_BASE_CONVERSION_THRESHOLD)
3413            return smallToString(radix);
3414 
3415         // Otherwise use recursive toString, which requires positive arguments.
3416         // The results will be concatenated into this StringBuilder
3417         StringBuilder sb = new StringBuilder();
3418         if (signum < 0) {
3419             toString(this.negate(), sb, radix, 0);
3420             sb.insert(0, '-');
3421         }
3422         else
3423             toString(this, sb, radix, 0);
3424 
3425         return sb.toString();
3426     }
3427 
3428     /** This method is used to perform toString when arguments are small. */
3429     private String smallToString(int radix) {
3430         if (signum == 0)
3431             return "0";
3432 
3433         // Compute upper bound on number of digit groups and allocate space
3434         int maxNumDigitGroups = (4*mag.length + 6)/7;
3435         String digitGroup[] = new String[maxNumDigitGroups];
3436 
3437         // Translate number to string, a digit group at a time
3438         BigInteger tmp = this.abs();
3439         int numGroups = 0;
3440         while (tmp.signum != 0) {
3441             BigInteger d = longRadix[radix];
3442 
3443             MutableBigInteger q = new MutableBigInteger(),
3444                               a = new MutableBigInteger(tmp.mag),
3445                               b = new MutableBigInteger(d.mag);
3446             MutableBigInteger r = a.divide(b, q);
3447             BigInteger q2 = q.toBigInteger(tmp.signum * d.signum);
3448             BigInteger r2 = r.toBigInteger(tmp.signum * d.signum);
3449 
3450             digitGroup[numGroups++] = Long.toString(r2.longValue(), radix);
3451             tmp = q2;
3452         }
3453 
3454         // Put sign (if any) and first digit group into result buffer
3455         StringBuilder buf = new StringBuilder(numGroups*digitsPerLong[radix]+1);
3456         if (signum<0)
3457             buf.append('-');
3458         buf.append(digitGroup[numGroups-1]);
3459 
3460         // Append remaining digit groups padded with leading zeros
3461         for (int i=numGroups-2; i>=0; i--) {
3462             // Prepend (any) leading zeros for this digit group
3463             int numLeadingZeros = digitsPerLong[radix]-digitGroup[i].length();
3464             if (numLeadingZeros != 0)
3465                 buf.append(zeros[numLeadingZeros]);
3466             buf.append(digitGroup[i]);
3467         }
3468         return buf.toString();
3469     }
3470 
3471     /**
3472      * Converts the specified BigInteger to a string and appends to
3473      * {@code sb}.  This implements the recursive Schoenhage algorithm
3474      * for base conversions.
3475      * <p/>
3476      * See Knuth, Donald,  _The Art of Computer Programming_, Vol. 2,
3477      * Answers to Exercises (4.4) Question 14.
3478      *
3479      * @param u      The number to convert to a string.
3480      * @param sb     The StringBuilder that will be appended to in place.
3481      * @param radix  The base to convert to.
3482      * @param digits The minimum number of digits to pad to.
3483      */
3484     private static void toString(BigInteger u, StringBuilder sb, int radix,
3485                                  int digits) {
3486         /* If we're smaller than a certain threshold, use the smallToString
3487            method, padding with leading zeroes when necessary. */
3488         if (u.mag.length <= SCHOENHAGE_BASE_CONVERSION_THRESHOLD) {
3489             String s = u.smallToString(radix);
3490 
3491             // Pad with internal zeros if necessary.
3492             // Don't pad if we're at the beginning of the string.
3493             if ((s.length() < digits) && (sb.length() > 0))
3494                 for (int i=s.length(); i<digits; i++) // May be a faster way to
3495                     sb.append('0');                    // do this?
3496 
3497             sb.append(s);
3498             return;
3499         }
3500 
3501         int b, n;
3502         b = u.bitLength();
3503 
3504         // Calculate a value for n in the equation radix^(2^n) = u
3505         // and subtract 1 from that value.  This is used to find the
3506         // cache index that contains the best value to divide u.
3507         n = (int) Math.round(Math.log(b * LOG_TWO / logCache[radix]) / LOG_TWO - 1.0);
3508         BigInteger v = getRadixConversionCache(radix, n);
3509         BigInteger[] results;
3510         results = u.divideAndRemainder(v);
3511 
3512         int expectedDigits = 1 << n;
3513 
3514         // Now recursively build the two halves of each number.
3515         toString(results[0], sb, radix, digits-expectedDigits);
3516         toString(results[1], sb, radix, expectedDigits);
3517     }
3518 
3519     /**
3520      * Returns the value radix^(2^exponent) from the cache.
3521      * If this value doesn't already exist in the cache, it is added.
3522      * <p/>
3523      * This could be changed to a more complicated caching method using
3524      * {@code Future}.
3525      */
3526     private static BigInteger getRadixConversionCache(int radix, int exponent) {
3527         BigInteger[] cacheLine = powerCache[radix]; // volatile read
3528         if (exponent < cacheLine.length) {
3529             return cacheLine[exponent];
3530         }
3531 
3532         int oldLength = cacheLine.length;
3533         cacheLine = Arrays.copyOf(cacheLine, exponent + 1);
3534         for (int i = oldLength; i <= exponent; i++) {
3535             cacheLine[i] = cacheLine[i - 1].pow(2);
3536         }
3537 
3538         BigInteger[][] pc = powerCache; // volatile read again
3539         if (exponent >= pc[radix].length) {
3540             pc = pc.clone();
3541             pc[radix] = cacheLine;
3542             powerCache = pc; // volatile write, publish
3543         }
3544         return cacheLine[exponent];
3545     }
3546 
3547     /* zero[i] is a string of i consecutive zeros. */
3548     private static String zeros[] = new String[64];
3549     static {
3550         zeros[63] =
3551             "000000000000000000000000000000000000000000000000000000000000000";
3552         for (int i=0; i<63; i++)
3553             zeros[i] = zeros[63].substring(0, i);
3554     }
3555 
3556     /**
3557      * Returns the decimal String representation of this BigInteger.
3558      * The digit-to-character mapping provided by
3559      * {@code Character.forDigit} is used, and a minus sign is
3560      * prepended if appropriate.  (This representation is compatible
3561      * with the {@link #BigInteger(String) (String)} constructor, and
3562      * allows for String concatenation with Java's + operator.)
3563      *
3564      * @return decimal String representation of this BigInteger.
3565      * @see    Character#forDigit
3566      * @see    #BigInteger(java.lang.String)
3567      */
3568     public String toString() {
3569         return toString(10);
3570     }
3571 
3572     /**
3573      * Returns a byte array containing the two's-complement
3574      * representation of this BigInteger.  The byte array will be in
3575      * <i>big-endian</i> byte-order: the most significant byte is in
3576      * the zeroth element.  The array will contain the minimum number
3577      * of bytes required to represent this BigInteger, including at
3578      * least one sign bit, which is {@code (ceil((this.bitLength() +
3579      * 1)/8))}.  (This representation is compatible with the
3580      * {@link #BigInteger(byte[]) (byte[])} constructor.)
3581      *
3582      * @return a byte array containing the two's-complement representation of
3583      *         this BigInteger.
3584      * @see    #BigInteger(byte[])
3585      */
3586     public byte[] toByteArray() {
3587         int byteLen = bitLength()/8 + 1;
3588         byte[] byteArray = new byte[byteLen];
3589 
3590         for (int i=byteLen-1, bytesCopied=4, nextInt=0, intIndex=0; i>=0; i--) {
3591             if (bytesCopied == 4) {
3592                 nextInt = getInt(intIndex++);
3593                 bytesCopied = 1;
3594             } else {
3595                 nextInt >>>= 8;
3596                 bytesCopied++;
3597             }
3598             byteArray[i] = (byte)nextInt;
3599         }
3600         return byteArray;
3601     }
3602 
3603     /**
3604      * Converts this BigInteger to an {@code int}.  This
3605      * conversion is analogous to a
3606      * <i>narrowing primitive conversion</i> from {@code long} to
3607      * {@code int} as defined in section 5.1.3 of
3608      * <cite>The Java&trade; Language Specification</cite>:
3609      * if this BigInteger is too big to fit in an
3610      * {@code int}, only the low-order 32 bits are returned.
3611      * Note that this conversion can lose information about the
3612      * overall magnitude of the BigInteger value as well as return a
3613      * result with the opposite sign.
3614      *
3615      * @return this BigInteger converted to an {@code int}.
3616      * @see #intValueExact()
3617      */
3618     public int intValue() {
3619         int result = 0;
3620         result = getInt(0);
3621         return result;
3622     }
3623 
3624     /**
3625      * Converts this BigInteger to a {@code long}.  This
3626      * conversion is analogous to a
3627      * <i>narrowing primitive conversion</i> from {@code long} to
3628      * {@code int} as defined in section 5.1.3 of
3629      * <cite>The Java&trade; Language Specification</cite>:
3630      * if this BigInteger is too big to fit in a
3631      * {@code long}, only the low-order 64 bits are returned.
3632      * Note that this conversion can lose information about the
3633      * overall magnitude of the BigInteger value as well as return a
3634      * result with the opposite sign.
3635      *
3636      * @return this BigInteger converted to a {@code long}.
3637      * @see #longValueExact()
3638      */
3639     public long longValue() {
3640         long result = 0;
3641 
3642         for (int i=1; i>=0; i--)
3643             result = (result << 32) + (getInt(i) & LONG_MASK);
3644         return result;
3645     }
3646 
3647     /**
3648      * Converts this BigInteger to a {@code float}.  This
3649      * conversion is similar to the
3650      * <i>narrowing primitive conversion</i> from {@code double} to
3651      * {@code float} as defined in section 5.1.3 of
3652      * <cite>The Java&trade; Language Specification</cite>:
3653      * if this BigInteger has too great a magnitude
3654      * to represent as a {@code float}, it will be converted to
3655      * {@link Float#NEGATIVE_INFINITY} or {@link
3656      * Float#POSITIVE_INFINITY} as appropriate.  Note that even when
3657      * the return value is finite, this conversion can lose
3658      * information about the precision of the BigInteger value.
3659      *
3660      * @return this BigInteger converted to a {@code float}.
3661      */
3662     public float floatValue() {
3663         if (signum == 0) {
3664             return 0.0f;
3665         }
3666 
3667         int exponent = ((mag.length - 1) << 5) + bitLengthForInt(mag[0]) - 1;
3668 
3669         // exponent == floor(log2(abs(this)))
3670         if (exponent < Long.SIZE - 1) {
3671             return longValue();
3672         } else if (exponent > Float.MAX_EXPONENT) {
3673             return signum > 0 ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3674         }
3675 
3676         /*
3677          * We need the top SIGNIFICAND_WIDTH bits, including the "implicit"
3678          * one bit. To make rounding easier, we pick out the top
3679          * SIGNIFICAND_WIDTH + 1 bits, so we have one to help us round up or
3680          * down. twiceSignifFloor will contain the top SIGNIFICAND_WIDTH + 1
3681          * bits, and signifFloor the top SIGNIFICAND_WIDTH.
3682          *
3683          * It helps to consider the real number signif = abs(this) *
3684          * 2^(SIGNIFICAND_WIDTH - 1 - exponent).
3685          */
3686         int shift = exponent - FloatConsts.SIGNIFICAND_WIDTH;
3687 
3688         int twiceSignifFloor;
3689         // twiceSignifFloor will be == abs().shiftRight(shift).intValue()
3690         // We do the shift into an int directly to improve performance.
3691 
3692         int nBits = shift & 0x1f;
3693         int nBits2 = 32 - nBits;
3694 
3695         if (nBits == 0) {
3696             twiceSignifFloor = mag[0];
3697         } else {
3698             twiceSignifFloor = mag[0] >>> nBits;
3699             if (twiceSignifFloor == 0) {
3700                 twiceSignifFloor = (mag[0] << nBits2) | (mag[1] >>> nBits);
3701             }
3702         }
3703 
3704         int signifFloor = twiceSignifFloor >> 1;
3705         signifFloor &= FloatConsts.SIGNIF_BIT_MASK; // remove the implied bit
3706 
3707         /*
3708          * We round up if either the fractional part of signif is strictly
3709          * greater than 0.5 (which is true if the 0.5 bit is set and any lower
3710          * bit is set), or if the fractional part of signif is >= 0.5 and
3711          * signifFloor is odd (which is true if both the 0.5 bit and the 1 bit
3712          * are set). This is equivalent to the desired HALF_EVEN rounding.
3713          */
3714         boolean increment = (twiceSignifFloor & 1) != 0
3715                 && ((signifFloor & 1) != 0 || abs().getLowestSetBit() < shift);
3716         int signifRounded = increment ? signifFloor + 1 : signifFloor;
3717         int bits = ((exponent + FloatConsts.EXP_BIAS))
3718                 << (FloatConsts.SIGNIFICAND_WIDTH - 1);
3719         bits += signifRounded;
3720         /*
3721          * If signifRounded == 2^24, we'd need to set all of the significand
3722          * bits to zero and add 1 to the exponent. This is exactly the behavior
3723          * we get from just adding signifRounded to bits directly. If the
3724          * exponent is Float.MAX_EXPONENT, we round up (correctly) to
3725          * Float.POSITIVE_INFINITY.
3726          */
3727         bits |= signum & FloatConsts.SIGN_BIT_MASK;
3728         return Float.intBitsToFloat(bits);
3729     }
3730 
3731     /**
3732      * Converts this BigInteger to a {@code double}.  This
3733      * conversion is similar to the
3734      * <i>narrowing primitive conversion</i> from {@code double} to
3735      * {@code float} as defined in section 5.1.3 of
3736      * <cite>The Java&trade; Language Specification</cite>:
3737      * if this BigInteger has too great a magnitude
3738      * to represent as a {@code double}, it will be converted to
3739      * {@link Double#NEGATIVE_INFINITY} or {@link
3740      * Double#POSITIVE_INFINITY} as appropriate.  Note that even when
3741      * the return value is finite, this conversion can lose
3742      * information about the precision of the BigInteger value.
3743      *
3744      * @return this BigInteger converted to a {@code double}.
3745      */
3746     public double doubleValue() {
3747         if (signum == 0) {
3748             return 0.0;
3749         }
3750 
3751         int exponent = ((mag.length - 1) << 5) + bitLengthForInt(mag[0]) - 1;
3752 
3753         // exponent == floor(log2(abs(this))Double)
3754         if (exponent < Long.SIZE - 1) {
3755             return longValue();
3756         } else if (exponent > Double.MAX_EXPONENT) {
3757             return signum > 0 ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3758         }
3759 
3760         /*
3761          * We need the top SIGNIFICAND_WIDTH bits, including the "implicit"
3762          * one bit. To make rounding easier, we pick out the top
3763          * SIGNIFICAND_WIDTH + 1 bits, so we have one to help us round up or
3764          * down. twiceSignifFloor will contain the top SIGNIFICAND_WIDTH + 1
3765          * bits, and signifFloor the top SIGNIFICAND_WIDTH.
3766          *
3767          * It helps to consider the real number signif = abs(this) *
3768          * 2^(SIGNIFICAND_WIDTH - 1 - exponent).
3769          */
3770         int shift = exponent - DoubleConsts.SIGNIFICAND_WIDTH;
3771 
3772         long twiceSignifFloor;
3773         // twiceSignifFloor will be == abs().shiftRight(shift).longValue()
3774         // We do the shift into a long directly to improve performance.
3775 
3776         int nBits = shift & 0x1f;
3777         int nBits2 = 32 - nBits;
3778 
3779         int highBits;
3780         int lowBits;
3781         if (nBits == 0) {
3782             highBits = mag[0];
3783             lowBits = mag[1];
3784         } else {
3785             highBits = mag[0] >>> nBits;
3786             lowBits = (mag[0] << nBits2) | (mag[1] >>> nBits);
3787             if (highBits == 0) {
3788                 highBits = lowBits;
3789                 lowBits = (mag[1] << nBits2) | (mag[2] >>> nBits);
3790             }
3791         }
3792 
3793         twiceSignifFloor = ((highBits & LONG_MASK) << 32)
3794                 | (lowBits & LONG_MASK);
3795 
3796         long signifFloor = twiceSignifFloor >> 1;
3797         signifFloor &= DoubleConsts.SIGNIF_BIT_MASK; // remove the implied bit
3798 
3799         /*
3800          * We round up if either the fractional part of signif is strictly
3801          * greater than 0.5 (which is true if the 0.5 bit is set and any lower
3802          * bit is set), or if the fractional part of signif is >= 0.5 and
3803          * signifFloor is odd (which is true if both the 0.5 bit and the 1 bit
3804          * are set). This is equivalent to the desired HALF_EVEN rounding.
3805          */
3806         boolean increment = (twiceSignifFloor & 1) != 0
3807                 && ((signifFloor & 1) != 0 || abs().getLowestSetBit() < shift);
3808         long signifRounded = increment ? signifFloor + 1 : signifFloor;
3809         long bits = (long) ((exponent + DoubleConsts.EXP_BIAS))
3810                 << (DoubleConsts.SIGNIFICAND_WIDTH - 1);
3811         bits += signifRounded;
3812         /*
3813          * If signifRounded == 2^53, we'd need to set all of the significand
3814          * bits to zero and add 1 to the exponent. This is exactly the behavior
3815          * we get from just adding signifRounded to bits directly. If the
3816          * exponent is Double.MAX_EXPONENT, we round up (correctly) to
3817          * Double.POSITIVE_INFINITY.
3818          */
3819         bits |= signum & DoubleConsts.SIGN_BIT_MASK;
3820         return Double.longBitsToDouble(bits);
3821     }
3822 
3823     /**
3824      * Returns a copy of the input array stripped of any leading zero bytes.
3825      */
3826     private static int[] stripLeadingZeroInts(int val[]) {
3827         int vlen = val.length;
3828         int keep;
3829 
3830         // Find first nonzero byte
3831         for (keep = 0; keep < vlen && val[keep] == 0; keep++)
3832             ;
3833         return java.util.Arrays.copyOfRange(val, keep, vlen);
3834     }
3835 
3836     /**
3837      * Returns the input array stripped of any leading zero bytes.
3838      * Since the source is trusted the copying may be skipped.
3839      */
3840     private static int[] trustedStripLeadingZeroInts(int val[]) {
3841         int vlen = val.length;
3842         int keep;
3843 
3844         // Find first nonzero byte
3845         for (keep = 0; keep < vlen && val[keep] == 0; keep++)
3846             ;
3847         return keep == 0 ? val : java.util.Arrays.copyOfRange(val, keep, vlen);
3848     }
3849 
3850     /**
3851      * Returns a copy of the input array stripped of any leading zero bytes.
3852      */
3853     private static int[] stripLeadingZeroBytes(byte a[]) {
3854         int byteLength = a.length;
3855         int keep;
3856 
3857         // Find first nonzero byte
3858         for (keep = 0; keep < byteLength && a[keep]==0; keep++)
3859             ;
3860 
3861         // Allocate new array and copy relevant part of input array
3862         int intLength = ((byteLength - keep) + 3) >>> 2;
3863         int[] result = new int[intLength];
3864         int b = byteLength - 1;
3865         for (int i = intLength-1; i >= 0; i--) {
3866             result[i] = a[b--] & 0xff;
3867             int bytesRemaining = b - keep + 1;
3868             int bytesToTransfer = Math.min(3, bytesRemaining);
3869             for (int j=8; j <= (bytesToTransfer << 3); j += 8)
3870                 result[i] |= ((a[b--] & 0xff) << j);
3871         }
3872         return result;
3873     }
3874 
3875     /**
3876      * Takes an array a representing a negative 2's-complement number and
3877      * returns the minimal (no leading zero bytes) unsigned whose value is -a.
3878      */
3879     private static int[] makePositive(byte a[]) {
3880         int keep, k;
3881         int byteLength = a.length;
3882 
3883         // Find first non-sign (0xff) byte of input
3884         for (keep=0; keep<byteLength && a[keep]==-1; keep++)
3885             ;
3886 
3887 
3888         /* Allocate output array.  If all non-sign bytes are 0x00, we must
3889          * allocate space for one extra output byte. */
3890         for (k=keep; k<byteLength && a[k]==0; k++)
3891             ;
3892 
3893         int extraByte = (k==byteLength) ? 1 : 0;
3894         int intLength = ((byteLength - keep + extraByte) + 3)/4;
3895         int result[] = new int[intLength];
3896 
3897         /* Copy one's complement of input into output, leaving extra
3898          * byte (if it exists) == 0x00 */
3899         int b = byteLength - 1;
3900         for (int i = intLength-1; i >= 0; i--) {
3901             result[i] = a[b--] & 0xff;
3902             int numBytesToTransfer = Math.min(3, b-keep+1);
3903             if (numBytesToTransfer < 0)
3904                 numBytesToTransfer = 0;
3905             for (int j=8; j <= 8*numBytesToTransfer; j += 8)
3906                 result[i] |= ((a[b--] & 0xff) << j);
3907 
3908             // Mask indicates which bits must be complemented
3909             int mask = -1 >>> (8*(3-numBytesToTransfer));
3910             result[i] = ~result[i] & mask;
3911         }
3912 
3913         // Add one to one's complement to generate two's complement
3914         for (int i=result.length-1; i>=0; i--) {
3915             result[i] = (int)((result[i] & LONG_MASK) + 1);
3916             if (result[i] != 0)
3917                 break;
3918         }
3919 
3920         return result;
3921     }
3922 
3923     /**
3924      * Takes an array a representing a negative 2's-complement number and
3925      * returns the minimal (no leading zero ints) unsigned whose value is -a.
3926      */
3927     private static int[] makePositive(int a[]) {
3928         int keep, j;
3929 
3930         // Find first non-sign (0xffffffff) int of input
3931         for (keep=0; keep<a.length && a[keep]==-1; keep++)
3932             ;
3933 
3934         /* Allocate output array.  If all non-sign ints are 0x00, we must
3935          * allocate space for one extra output int. */
3936         for (j=keep; j<a.length && a[j]==0; j++)
3937             ;
3938         int extraInt = (j==a.length ? 1 : 0);
3939         int result[] = new int[a.length - keep + extraInt];
3940 
3941         /* Copy one's complement of input into output, leaving extra
3942          * int (if it exists) == 0x00 */
3943         for (int i = keep; i<a.length; i++)
3944             result[i - keep + extraInt] = ~a[i];
3945 
3946         // Add one to one's complement to generate two's complement
3947         for (int i=result.length-1; ++result[i]==0; i--)
3948             ;
3949 
3950         return result;
3951     }
3952 
3953     /*
3954      * The following two arrays are used for fast String conversions.  Both
3955      * are indexed by radix.  The first is the number of digits of the given
3956      * radix that can fit in a Java long without "going negative", i.e., the
3957      * highest integer n such that radix**n < 2**63.  The second is the
3958      * "long radix" that tears each number into "long digits", each of which
3959      * consists of the number of digits in the corresponding element in
3960      * digitsPerLong (longRadix[i] = i**digitPerLong[i]).  Both arrays have
3961      * nonsense values in their 0 and 1 elements, as radixes 0 and 1 are not
3962      * used.
3963      */
3964     private static int digitsPerLong[] = {0, 0,
3965         62, 39, 31, 27, 24, 22, 20, 19, 18, 18, 17, 17, 16, 16, 15, 15, 15, 14,
3966         14, 14, 14, 13, 13, 13, 13, 13, 13, 12, 12, 12, 12, 12, 12, 12, 12};
3967 
3968     private static BigInteger longRadix[] = {null, null,
3969         valueOf(0x4000000000000000L), valueOf(0x383d9170b85ff80bL),
3970         valueOf(0x4000000000000000L), valueOf(0x6765c793fa10079dL),
3971         valueOf(0x41c21cb8e1000000L), valueOf(0x3642798750226111L),
3972         valueOf(0x1000000000000000L), valueOf(0x12bf307ae81ffd59L),
3973         valueOf( 0xde0b6b3a7640000L), valueOf(0x4d28cb56c33fa539L),
3974         valueOf(0x1eca170c00000000L), valueOf(0x780c7372621bd74dL),
3975         valueOf(0x1e39a5057d810000L), valueOf(0x5b27ac993df97701L),
3976         valueOf(0x1000000000000000L), valueOf(0x27b95e997e21d9f1L),
3977         valueOf(0x5da0e1e53c5c8000L), valueOf( 0xb16a458ef403f19L),
3978         valueOf(0x16bcc41e90000000L), valueOf(0x2d04b7fdd9c0ef49L),
3979         valueOf(0x5658597bcaa24000L), valueOf( 0x6feb266931a75b7L),
3980         valueOf( 0xc29e98000000000L), valueOf(0x14adf4b7320334b9L),
3981         valueOf(0x226ed36478bfa000L), valueOf(0x383d9170b85ff80bL),
3982         valueOf(0x5a3c23e39c000000L), valueOf( 0x4e900abb53e6b71L),
3983         valueOf( 0x7600ec618141000L), valueOf( 0xaee5720ee830681L),
3984         valueOf(0x1000000000000000L), valueOf(0x172588ad4f5f0981L),
3985         valueOf(0x211e44f7d02c1000L), valueOf(0x2ee56725f06e5c71L),
3986         valueOf(0x41c21cb8e1000000L)};
3987 
3988     /*
3989      * These two arrays are the integer analogue of above.
3990      */
3991     private static int digitsPerInt[] = {0, 0, 30, 19, 15, 13, 11,
3992         11, 10, 9, 9, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6,
3993         6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5};
3994 
3995     private static int intRadix[] = {0, 0,
3996         0x40000000, 0x4546b3db, 0x40000000, 0x48c27395, 0x159fd800,
3997         0x75db9c97, 0x40000000, 0x17179149, 0x3b9aca00, 0xcc6db61,
3998         0x19a10000, 0x309f1021, 0x57f6c100, 0xa2f1b6f,  0x10000000,
3999         0x18754571, 0x247dbc80, 0x3547667b, 0x4c4b4000, 0x6b5a6e1d,
4000         0x6c20a40,  0x8d2d931,  0xb640000,  0xe8d4a51,  0x1269ae40,
4001         0x17179149, 0x1cb91000, 0x23744899, 0x2b73a840, 0x34e63b41,
4002         0x40000000, 0x4cfa3cc1, 0x5c13d840, 0x6d91b519, 0x39aa400
4003     };
4004 
4005     /**
4006      * These routines provide access to the two's complement representation
4007      * of BigIntegers.
4008      */
4009 
4010     /**
4011      * Returns the length of the two's complement representation in ints,
4012      * including space for at least one sign bit.
4013      */
4014     private int intLength() {
4015         return (bitLength() >>> 5) + 1;
4016     }
4017 
4018     /* Returns sign bit */
4019     private int signBit() {
4020         return signum < 0 ? 1 : 0;
4021     }
4022 
4023     /* Returns an int of sign bits */
4024     private int signInt() {
4025         return signum < 0 ? -1 : 0;
4026     }
4027 
4028     /**
4029      * Returns the specified int of the little-endian two's complement
4030      * representation (int 0 is the least significant).  The int number can
4031      * be arbitrarily high (values are logically preceded by infinitely many
4032      * sign ints).
4033      */
4034     private int getInt(int n) {
4035         if (n < 0)
4036             return 0;
4037         if (n >= mag.length)
4038             return signInt();
4039 
4040         int magInt = mag[mag.length-n-1];
4041 
4042         return (signum >= 0 ? magInt :
4043                 (n <= firstNonzeroIntNum() ? -magInt : ~magInt));
4044     }
4045 
4046     /**
4047      * Returns the index of the int that contains the first nonzero int in the
4048      * little-endian binary representation of the magnitude (int 0 is the
4049      * least significant). If the magnitude is zero, return value is undefined.
4050      */
4051     private int firstNonzeroIntNum() {
4052         int fn = firstNonzeroIntNum - 2;
4053         if (fn == -2) { // firstNonzeroIntNum not initialized yet
4054             fn = 0;
4055 
4056             // Search for the first nonzero int
4057             int i;
4058             int mlen = mag.length;
4059             for (i = mlen - 1; i >= 0 && mag[i] == 0; i--)
4060                 ;
4061             fn = mlen - i - 1;
4062             firstNonzeroIntNum = fn + 2; // offset by two to initialize
4063         }
4064         return fn;
4065     }
4066 
4067     /** use serialVersionUID from JDK 1.1. for interoperability */
4068     private static final long serialVersionUID = -8287574255936472291L;
4069 
4070     /**
4071      * Serializable fields for BigInteger.
4072      *
4073      * @serialField signum  int
4074      *              signum of this BigInteger.
4075      * @serialField magnitude int[]
4076      *              magnitude array of this BigInteger.
4077      * @serialField bitCount  int
4078      *              number of bits in this BigInteger
4079      * @serialField bitLength int
4080      *              the number of bits in the minimal two's-complement
4081      *              representation of this BigInteger
4082      * @serialField lowestSetBit int
4083      *              lowest set bit in the twos complement representation
4084      */
4085     private static final ObjectStreamField[] serialPersistentFields = {
4086         new ObjectStreamField("signum", Integer.TYPE),
4087         new ObjectStreamField("magnitude", byte[].class),
4088         new ObjectStreamField("bitCount", Integer.TYPE),
4089         new ObjectStreamField("bitLength", Integer.TYPE),
4090         new ObjectStreamField("firstNonzeroByteNum", Integer.TYPE),
4091         new ObjectStreamField("lowestSetBit", Integer.TYPE)
4092         };
4093 
4094     /**
4095      * Reconstitute the {@code BigInteger} instance from a stream (that is,
4096      * deserialize it). The magnitude is read in as an array of bytes
4097      * for historical reasons, but it is converted to an array of ints
4098      * and the byte array is discarded.
4099      * Note:
4100      * The current convention is to initialize the cache fields, bitCount,
4101      * bitLength and lowestSetBit, to 0 rather than some other marker value.
4102      * Therefore, no explicit action to set these fields needs to be taken in
4103      * readObject because those fields already have a 0 value be default since
4104      * defaultReadObject is not being used.
4105      */
4106     private void readObject(java.io.ObjectInputStream s)
4107         throws java.io.IOException, ClassNotFoundException {
4108         /*
4109          * In order to maintain compatibility with previous serialized forms,
4110          * the magnitude of a BigInteger is serialized as an array of bytes.
4111          * The magnitude field is used as a temporary store for the byte array
4112          * that is deserialized. The cached computation fields should be
4113          * transient but are serialized for compatibility reasons.
4114          */
4115 
4116         // prepare to read the alternate persistent fields
4117         ObjectInputStream.GetField fields = s.readFields();
4118 
4119         // Read the alternate persistent fields that we care about
4120         int sign = fields.get("signum", -2);
4121         byte[] magnitude = (byte[])fields.get("magnitude", null);
4122 
4123         // Validate signum
4124         if (sign < -1 || sign > 1) {
4125             String message = "BigInteger: Invalid signum value";
4126             if (fields.defaulted("signum"))
4127                 message = "BigInteger: Signum not present in stream";
4128             throw new java.io.StreamCorruptedException(message);
4129         }
4130         if ((magnitude.length == 0) != (sign == 0)) {
4131             String message = "BigInteger: signum-magnitude mismatch";
4132             if (fields.defaulted("magnitude"))
4133                 message = "BigInteger: Magnitude not present in stream";
4134             throw new java.io.StreamCorruptedException(message);
4135         }
4136 
4137         // Commit final fields via Unsafe
4138         UnsafeHolder.putSign(this, sign);
4139 
4140         // Calculate mag field from magnitude and discard magnitude
4141         UnsafeHolder.putMag(this, stripLeadingZeroBytes(magnitude));
4142     }
4143 
4144     // Support for resetting final fields while deserializing
4145     private static class UnsafeHolder {
4146         private static final sun.misc.Unsafe unsafe;
4147         private static final long signumOffset;
4148         private static final long magOffset;
4149         static {
4150             try {
4151                 unsafe = sun.misc.Unsafe.getUnsafe();
4152                 signumOffset = unsafe.objectFieldOffset
4153                     (BigInteger.class.getDeclaredField("signum"));
4154                 magOffset = unsafe.objectFieldOffset
4155                     (BigInteger.class.getDeclaredField("mag"));
4156             } catch (Exception ex) {
4157                 throw new ExceptionInInitializerError(ex);
4158             }
4159         }
4160 
4161         static void putSign(BigInteger bi, int sign) {
4162             unsafe.putIntVolatile(bi, signumOffset, sign);
4163         }
4164 
4165         static void putMag(BigInteger bi, int[] magnitude) {
4166             unsafe.putObjectVolatile(bi, magOffset, magnitude);
4167         }
4168     }
4169 
4170     /**
4171      * Save the {@code BigInteger} instance to a stream.
4172      * The magnitude of a BigInteger is serialized as a byte array for
4173      * historical reasons.
4174      *
4175      * @serialData two necessary fields are written as well as obsolete
4176      *             fields for compatibility with older versions.
4177      */
4178     private void writeObject(ObjectOutputStream s) throws IOException {
4179         // set the values of the Serializable fields
4180         ObjectOutputStream.PutField fields = s.putFields();
4181         fields.put("signum", signum);
4182         fields.put("magnitude", magSerializedForm());
4183         // The values written for cached fields are compatible with older
4184         // versions, but are ignored in readObject so don't otherwise matter.
4185         fields.put("bitCount", -1);
4186         fields.put("bitLength", -1);
4187         fields.put("lowestSetBit", -2);
4188         fields.put("firstNonzeroByteNum", -2);
4189 
4190         // save them
4191         s.writeFields();
4192 }
4193 
4194     /**
4195      * Returns the mag array as an array of bytes.
4196      */
4197     private byte[] magSerializedForm() {
4198         int len = mag.length;
4199 
4200         int bitLen = (len == 0 ? 0 : ((len - 1) << 5) + bitLengthForInt(mag[0]));
4201         int byteLen = (bitLen + 7) >>> 3;
4202         byte[] result = new byte[byteLen];
4203 
4204         for (int i = byteLen - 1, bytesCopied = 4, intIndex = len - 1, nextInt = 0;
4205              i>=0; i--) {
4206             if (bytesCopied == 4) {
4207                 nextInt = mag[intIndex--];
4208                 bytesCopied = 1;
4209             } else {
4210                 nextInt >>>= 8;
4211                 bytesCopied++;
4212             }
4213             result[i] = (byte)nextInt;
4214         }
4215         return result;
4216     }
4217 
4218     /**
4219      * Converts this {@code BigInteger} to a {@code long}, checking
4220      * for lost information.  If the value of this {@code BigInteger}
4221      * is out of the range of the {@code long} type, then an
4222      * {@code ArithmeticException} is thrown.
4223      *
4224      * @return this {@code BigInteger} converted to a {@code long}.
4225      * @throws ArithmeticException if the value of {@code this} will
4226      * not exactly fit in a {@code long}.
4227      * @see BigInteger#longValue
4228      * @since  1.8
4229      */
4230     public long longValueExact() {
4231         if (mag.length <= 2 && bitLength() <= 63)
4232             return longValue();
4233         else
4234             throw new ArithmeticException("BigInteger out of long range");
4235     }
4236 
4237     /**
4238      * Converts this {@code BigInteger} to an {@code int}, checking
4239      * for lost information.  If the value of this {@code BigInteger}
4240      * is out of the range of the {@code int} type, then an
4241      * {@code ArithmeticException} is thrown.
4242      *
4243      * @return this {@code BigInteger} converted to an {@code int}.
4244      * @throws ArithmeticException if the value of {@code this} will
4245      * not exactly fit in a {@code int}.
4246      * @see BigInteger#intValue
4247      * @since  1.8
4248      */
4249     public int intValueExact() {
4250         if (mag.length <= 1 && bitLength() <= 31)
4251             return intValue();
4252         else
4253             throw new ArithmeticException("BigInteger out of int range");
4254     }
4255 
4256     /**
4257      * Converts this {@code BigInteger} to a {@code short}, checking
4258      * for lost information.  If the value of this {@code BigInteger}
4259      * is out of the range of the {@code short} type, then an
4260      * {@code ArithmeticException} is thrown.
4261      *
4262      * @return this {@code BigInteger} converted to a {@code short}.
4263      * @throws ArithmeticException if the value of {@code this} will
4264      * not exactly fit in a {@code short}.
4265      * @see BigInteger#shortValue
4266      * @since  1.8
4267      */
4268     public short shortValueExact() {
4269         if (mag.length <= 1 && bitLength() <= 31) {
4270             int value = intValue();
4271             if (value >= Short.MIN_VALUE && value <= Short.MAX_VALUE)
4272                 return shortValue();
4273         }
4274         throw new ArithmeticException("BigInteger out of short range");
4275     }
4276 
4277     /**
4278      * Converts this {@code BigInteger} to a {@code byte}, checking
4279      * for lost information.  If the value of this {@code BigInteger}
4280      * is out of the range of the {@code byte} type, then an
4281      * {@code ArithmeticException} is thrown.
4282      *
4283      * @return this {@code BigInteger} converted to a {@code byte}.
4284      * @throws ArithmeticException if the value of {@code this} will
4285      * not exactly fit in a {@code byte}.
4286      * @see BigInteger#byteValue
4287      * @since  1.8
4288      */
4289     public byte byteValueExact() {
4290         if (mag.length <= 1 && bitLength() <= 31) {
4291             int value = intValue();
4292             if (value >= Byte.MIN_VALUE && value <= Byte.MAX_VALUE)
4293                 return byteValue();
4294         }
4295         throw new ArithmeticException("BigInteger out of byte range");
4296     }
4297 }