1 /*
   2  * Portions Copyright 1996-2007 Sun Microsystems, Inc.  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.  Sun designates this
   8  * particular file as subject to the "Classpath" exception as provided
   9  * by Sun 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 Sun Microsystems, Inc., 4150 Network Circle, Santa Clara,
  22  * CA 95054 USA or visit www.sun.com if you need additional information or
  23  * have any questions.
  24  */
  25 
  26 /*
  27  * Portions Copyright (c) 1995  Colin Plumb.  All rights reserved.
  28  */
  29 
  30 package java.math;
  31 
  32 import java.util.Random;
  33 import java.io.*;
  34 
  35 /**
  36  * Immutable arbitrary-precision integers.  All operations behave as if
  37  * BigIntegers were represented in two's-complement notation (like Java's
  38  * primitive integer types).  BigInteger provides analogues to all of Java's
  39  * primitive integer operators, and all relevant methods from java.lang.Math.
  40  * Additionally, BigInteger provides operations for modular arithmetic, GCD
  41  * calculation, primality testing, prime generation, bit manipulation,
  42  * and a few other miscellaneous operations.
  43  *
  44  * <p>Semantics of arithmetic operations exactly mimic those of Java's integer
  45  * arithmetic operators, as defined in <i>The Java Language Specification</i>.
  46  * For example, division by zero throws an {@code ArithmeticException}, and
  47  * division of a negative by a positive yields a negative (or zero) remainder.
  48  * All of the details in the Spec concerning overflow are ignored, as
  49  * BigIntegers are made as large as necessary to accommodate the results of an
  50  * operation.
  51  *
  52  * <p>Semantics of shift operations extend those of Java's shift operators
  53  * to allow for negative shift distances.  A right-shift with a negative
  54  * shift distance results in a left shift, and vice-versa.  The unsigned
  55  * right shift operator ({@code >>>}) is omitted, as this operation makes
  56  * little sense in combination with the "infinite word size" abstraction
  57  * provided by this class.
  58  *
  59  * <p>Semantics of bitwise logical operations exactly mimic those of Java's
  60  * bitwise integer operators.  The binary operators ({@code and},
  61  * {@code or}, {@code xor}) implicitly perform sign extension on the shorter
  62  * of the two operands prior to performing the operation.
  63  *
  64  * <p>Comparison operations perform signed integer comparisons, analogous to
  65  * those performed by Java's relational and equality operators.
  66  *
  67  * <p>Modular arithmetic operations are provided to compute residues, perform
  68  * exponentiation, and compute multiplicative inverses.  These methods always
  69  * return a non-negative result, between {@code 0} and {@code (modulus - 1)},
  70  * inclusive.
  71  *
  72  * <p>Bit operations operate on a single bit of the two's-complement
  73  * representation of their operand.  If necessary, the operand is sign-
  74  * extended so that it contains the designated bit.  None of the single-bit
  75  * operations can produce a BigInteger with a different sign from the
  76  * BigInteger being operated on, as they affect only a single bit, and the
  77  * "infinite word size" abstraction provided by this class ensures that there
  78  * are infinitely many "virtual sign bits" preceding each BigInteger.
  79  *
  80  * <p>For the sake of brevity and clarity, pseudo-code is used throughout the
  81  * descriptions of BigInteger methods.  The pseudo-code expression
  82  * {@code (i + j)} is shorthand for "a BigInteger whose value is
  83  * that of the BigInteger {@code i} plus that of the BigInteger {@code j}."
  84  * The pseudo-code expression {@code (i == j)} is shorthand for
  85  * "{@code true} if and only if the BigInteger {@code i} represents the same
  86  * value as the BigInteger {@code j}."  Other pseudo-code expressions are
  87  * interpreted similarly.
  88  *
  89  * <p>All methods and constructors in this class throw
  90  * {@code NullPointerException} when passed
  91  * a null object reference for any input parameter.
  92  *
  93  * @see     BigDecimal
  94  * @author  Josh Bloch
  95  * @author  Michael McCloskey
  96  * @since JDK1.1
  97  */
  98 
  99 public class BigInteger extends Number implements Comparable<BigInteger> {
 100     /**
 101      * The signum of this BigInteger: -1 for negative, 0 for zero, or
 102      * 1 for positive.  Note that the BigInteger zero <i>must</i> have
 103      * a signum of 0.  This is necessary to ensures that there is exactly one
 104      * representation for each BigInteger value.
 105      *
 106      * @serial
 107      */
 108     int signum;
 109 
 110     /**
 111      * The magnitude of this BigInteger, in <i>big-endian</i> order: the
 112      * zeroth element of this array is the most-significant int of the
 113      * magnitude.  The magnitude must be "minimal" in that the most-significant
 114      * int ({@code mag[0]}) must be non-zero.  This is necessary to
 115      * ensure that there is exactly one representation for each BigInteger
 116      * value.  Note that this implies that the BigInteger zero has a
 117      * zero-length mag array.
 118      */
 119     int[] mag;
 120 
 121     // These "redundant fields" are initialized with recognizable nonsense
 122     // values, and cached the first time they are needed (or never, if they
 123     // aren't needed).
 124 
 125     /**
 126      * The bitCount of this BigInteger, as returned by bitCount(), or -1
 127      * (either value is acceptable).
 128      *
 129      * @serial
 130      * @see #bitCount
 131      */
 132     private int bitCount =  -1;
 133 
 134     /**
 135      * The bitLength of this BigInteger, as returned by bitLength(), or -1
 136      * (either value is acceptable).
 137      *
 138      * @serial
 139      * @see #bitLength()
 140      */
 141     private int bitLength = -1;
 142 
 143     /**
 144      * The lowest set bit of this BigInteger, as returned by getLowestSetBit(),
 145      * or -2 (either value is acceptable).
 146      *
 147      * @serial
 148      * @see #getLowestSetBit
 149      */
 150     private int lowestSetBit = -2;
 151 
 152     /**
 153      * The index of the lowest-order byte in the magnitude of this BigInteger
 154      * that contains a nonzero byte, or -2 (either value is acceptable).  The
 155      * least significant byte has int-number 0, the next byte in order of
 156      * increasing significance has byte-number 1, and so forth.
 157      *
 158      * @serial
 159      */
 160     private int firstNonzeroByteNum = -2;
 161 
 162     /**
 163      * The index of the lowest-order int in the magnitude of this BigInteger
 164      * that contains a nonzero int, or -2 (either value is acceptable).  The
 165      * least significant int has int-number 0, the next int in order of
 166      * increasing significance has int-number 1, and so forth.
 167      */
 168     private int firstNonzeroIntNum = -2;
 169 
 170     /**
 171      * This mask is used to obtain the value of an int as if it were unsigned.
 172      */
 173     private final static long LONG_MASK = 0xffffffffL;
 174 
 175     //Constructors
 176 
 177     /**
 178      * Translates a byte array containing the two's-complement binary
 179      * representation of a BigInteger into a BigInteger.  The input array is
 180      * assumed to be in <i>big-endian</i> byte-order: the most significant
 181      * byte is in the zeroth element.
 182      *
 183      * @param  val big-endian two's-complement binary representation of
 184      *         BigInteger.
 185      * @throws NumberFormatException {@code val} is zero bytes long.
 186      */
 187     public BigInteger(byte[] val) {
 188         if (val.length == 0)
 189             throw new NumberFormatException("Zero length BigInteger");
 190 
 191         if (val[0] < 0) {
 192             mag = makePositive(val);
 193             signum = -1;
 194         } else {
 195             mag = stripLeadingZeroBytes(val);
 196             signum = (mag.length == 0 ? 0 : 1);
 197         }
 198     }
 199 
 200     /**
 201      * This private constructor translates an int array containing the
 202      * two's-complement binary representation of a BigInteger into a
 203      * BigInteger. The input array is assumed to be in <i>big-endian</i>
 204      * int-order: the most significant int is in the zeroth element.
 205      */
 206     private BigInteger(int[] val) {
 207         if (val.length == 0)
 208             throw new NumberFormatException("Zero length BigInteger");
 209 
 210         if (val[0] < 0) {
 211             mag = makePositive(val);
 212             signum = -1;
 213         } else {
 214             mag = trustedStripLeadingZeroInts(val);
 215             signum = (mag.length == 0 ? 0 : 1);
 216         }
 217     }
 218 
 219     /**
 220      * Translates the sign-magnitude representation of a BigInteger into a
 221      * BigInteger.  The sign is represented as an integer signum value: -1 for
 222      * negative, 0 for zero, or 1 for positive.  The magnitude is a byte array
 223      * in <i>big-endian</i> byte-order: the most significant byte is in the
 224      * zeroth element.  A zero-length magnitude array is permissible, and will
 225      * result in a BigInteger value of 0, whether signum is -1, 0 or 1.
 226      *
 227      * @param  signum signum of the number (-1 for negative, 0 for zero, 1
 228      *         for positive).
 229      * @param  magnitude big-endian binary representation of the magnitude of
 230      *         the number.
 231      * @throws NumberFormatException {@code signum} is not one of the three
 232      *         legal values (-1, 0, and 1), or {@code signum} is 0 and
 233      *         {@code magnitude} contains one or more non-zero bytes.
 234      */
 235     public BigInteger(int signum, byte[] magnitude) {
 236         this.mag = stripLeadingZeroBytes(magnitude);
 237 
 238         if (signum < -1 || signum > 1)
 239             throw(new NumberFormatException("Invalid signum value"));
 240 
 241         if (this.mag.length==0) {
 242             this.signum = 0;
 243         } else {
 244             if (signum == 0)
 245                 throw(new NumberFormatException("signum-magnitude mismatch"));
 246             this.signum = signum;
 247         }
 248     }
 249 
 250     /**
 251      * A constructor for internal use that translates the sign-magnitude
 252      * representation of a BigInteger into a BigInteger. It checks the
 253      * arguments and copies the magnitude so this constructor would be
 254      * safe for external use.
 255      */
 256     private BigInteger(int signum, int[] magnitude) {
 257         this.mag = stripLeadingZeroInts(magnitude);
 258 
 259         if (signum < -1 || signum > 1)
 260             throw(new NumberFormatException("Invalid signum value"));
 261 
 262         if (this.mag.length==0) {
 263             this.signum = 0;
 264         } else {
 265             if (signum == 0)
 266                 throw(new NumberFormatException("signum-magnitude mismatch"));
 267             this.signum = signum;
 268         }
 269     }
 270 
 271     /**
 272      * Translates the String representation of a BigInteger in the
 273      * specified radix into a BigInteger.  The String representation
 274      * consists of an optional minus or plus sign followed by a
 275      * sequence of one or more digits in the specified radix.  The
 276      * character-to-digit mapping is provided by {@code
 277      * Character.digit}.  The String may not contain any extraneous
 278      * characters (whitespace, for example).
 279      *
 280      * @param val String representation of BigInteger.
 281      * @param radix radix to be used in interpreting {@code val}.
 282      * @throws NumberFormatException {@code val} is not a valid representation
 283      *         of a BigInteger in the specified radix, or {@code radix} is
 284      *         outside the range from {@link Character#MIN_RADIX} to
 285      *         {@link Character#MAX_RADIX}, inclusive.
 286      * @see    Character#digit
 287      */
 288     public BigInteger(String val, int radix) {
 289         int cursor = 0, numDigits;
 290         int len = val.length();
 291 
 292         if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
 293             throw new NumberFormatException("Radix out of range");
 294         if (val.length() == 0)
 295             throw new NumberFormatException("Zero length BigInteger");
 296 
 297         // Check for at most one leading sign
 298         signum = 1;
 299         int index1 = val.lastIndexOf('-');
 300         int index2 = val.lastIndexOf('+');
 301         if ((index1 + index2) <= -1) {
 302             // No leading sign character or at most one leading sign character
 303             if (index1 == 0 || index2 == 0) {
 304                 cursor = 1;
 305                 if (val.length() == 1)
 306                     throw new NumberFormatException("Zero length BigInteger");
 307             }
 308             if (index1 == 0)
 309                 signum = -1;
 310         } else
 311             throw new NumberFormatException("Illegal embedded sign character");
 312 
 313         // Skip leading zeros and compute number of digits in magnitude
 314         while (cursor < len &&
 315                Character.digit(val.charAt(cursor), radix) == 0)
 316             cursor++;
 317         if (cursor == len) {
 318             signum = 0;
 319             mag = ZERO.mag;
 320             return;
 321         } else {
 322             numDigits = len - cursor;
 323         }
 324 
 325         // Pre-allocate array of expected size. May be too large but can
 326         // never be too small. Typically exact.
 327         int numBits = (int)(((numDigits * bitsPerDigit[radix]) >>> 10) + 1);
 328         int numWords = (numBits + 31) /32;
 329         mag = new int[numWords];
 330 
 331         // Process first (potentially short) digit group
 332         int firstGroupLen = numDigits % digitsPerInt[radix];
 333         if (firstGroupLen == 0)
 334             firstGroupLen = digitsPerInt[radix];
 335         String group = val.substring(cursor, cursor += firstGroupLen);
 336         mag[mag.length - 1] = Integer.parseInt(group, radix);
 337         if (mag[mag.length - 1] < 0)
 338             throw new NumberFormatException("Illegal digit");
 339 
 340         // Process remaining digit groups
 341         int superRadix = intRadix[radix];
 342         int groupVal = 0;
 343         while (cursor < val.length()) {
 344             group = val.substring(cursor, cursor += digitsPerInt[radix]);
 345             groupVal = Integer.parseInt(group, radix);
 346             if (groupVal < 0)
 347                 throw new NumberFormatException("Illegal digit");
 348             destructiveMulAdd(mag, superRadix, groupVal);
 349         }
 350         // Required for cases where the array was overallocated.
 351         mag = trustedStripLeadingZeroInts(mag);
 352     }
 353 
 354     // Constructs a new BigInteger using a char array with radix=10
 355     BigInteger(char[] val) {
 356         int cursor = 0, numDigits;
 357         int len = val.length;
 358 
 359         // Check for leading minus sign
 360         signum = 1;
 361         if (val[0] == '-') {
 362             if (len == 1)
 363                 throw new NumberFormatException("Zero length BigInteger");
 364             signum = -1;
 365             cursor = 1;
 366         } else if (val[0] == '+') {
 367             if (len == 1)
 368                 throw new NumberFormatException("Zero length BigInteger");
 369             cursor = 1;
 370         }
 371 
 372         // Skip leading zeros and compute number of digits in magnitude
 373         while (cursor < len && Character.digit(val[cursor], 10) == 0)
 374             cursor++;
 375         if (cursor == len) {
 376             signum = 0;
 377             mag = ZERO.mag;
 378             return;
 379         } else {
 380             numDigits = len - cursor;
 381         }
 382 
 383         // Pre-allocate array of expected size
 384         int numWords;
 385         if (len < 10) {
 386             numWords = 1;
 387         } else {
 388             int numBits = (int)(((numDigits * bitsPerDigit[10]) >>> 10) + 1);
 389             numWords = (numBits + 31) /32;
 390         }
 391         mag = new int[numWords];
 392 
 393         // Process first (potentially short) digit group
 394         int firstGroupLen = numDigits % digitsPerInt[10];
 395         if (firstGroupLen == 0)
 396             firstGroupLen = digitsPerInt[10];
 397         mag[mag.length-1] = parseInt(val, cursor,  cursor += firstGroupLen);
 398 
 399         // Process remaining digit groups
 400         while (cursor < len) {
 401             int groupVal = parseInt(val, cursor, cursor += digitsPerInt[10]);
 402             destructiveMulAdd(mag, intRadix[10], groupVal);
 403         }
 404         mag = trustedStripLeadingZeroInts(mag);
 405     }
 406 
 407     // Create an integer with the digits between the two indexes
 408     // Assumes start < end. The result may be negative, but it
 409     // is to be treated as an unsigned value.
 410     private int parseInt(char[] source, int start, int end) {
 411         int result = Character.digit(source[start++], 10);
 412         if (result == -1)
 413             throw new NumberFormatException(new String(source));
 414 
 415         for (int index = start; index<end; index++) {
 416             int nextVal = Character.digit(source[index], 10);
 417             if (nextVal == -1)
 418                 throw new NumberFormatException(new String(source));
 419             result = 10*result + nextVal;
 420         }
 421 
 422         return result;
 423     }
 424 
 425     // bitsPerDigit in the given radix times 1024
 426     // Rounded up to avoid underallocation.
 427     private static long bitsPerDigit[] = { 0, 0,
 428         1024, 1624, 2048, 2378, 2648, 2875, 3072, 3247, 3402, 3543, 3672,
 429         3790, 3899, 4001, 4096, 4186, 4271, 4350, 4426, 4498, 4567, 4633,
 430         4696, 4756, 4814, 4870, 4923, 4975, 5025, 5074, 5120, 5166, 5210,
 431                                            5253, 5295};
 432 
 433     // Multiply x array times word y in place, and add word z
 434     private static void destructiveMulAdd(int[] x, int y, int z) {
 435         // Perform the multiplication word by word
 436         long ylong = y & LONG_MASK;
 437         long zlong = z & LONG_MASK;
 438         int len = x.length;
 439 
 440         long product = 0;
 441         long carry = 0;
 442         for (int i = len-1; i >= 0; i--) {
 443             product = ylong * (x[i] & LONG_MASK) + carry;
 444             x[i] = (int)product;
 445             carry = product >>> 32;
 446         }
 447 
 448         // Perform the addition
 449         long sum = (x[len-1] & LONG_MASK) + zlong;
 450         x[len-1] = (int)sum;
 451         carry = sum >>> 32;
 452         for (int i = len-2; i >= 0; i--) {
 453             sum = (x[i] & LONG_MASK) + carry;
 454             x[i] = (int)sum;
 455             carry = sum >>> 32;
 456         }
 457     }
 458 
 459     /**
 460      * Translates the decimal String representation of a BigInteger into a
 461      * BigInteger.  The String representation consists of an optional minus
 462      * sign followed by a sequence of one or more decimal digits.  The
 463      * character-to-digit mapping is provided by {@code Character.digit}.
 464      * The String may not contain any extraneous characters (whitespace, for
 465      * example).
 466      *
 467      * @param val decimal String representation of BigInteger.
 468      * @throws NumberFormatException {@code val} is not a valid representation
 469      *         of a BigInteger.
 470      * @see    Character#digit
 471      */
 472     public BigInteger(String val) {
 473         this(val, 10);
 474     }
 475 
 476     /**
 477      * Constructs a randomly generated BigInteger, uniformly distributed over
 478      * the range {@code 0} to (2<sup>{@code numBits}</sup> - 1), inclusive.
 479      * The uniformity of the distribution assumes that a fair source of random
 480      * bits is provided in {@code rnd}.  Note that this constructor always
 481      * constructs a non-negative BigInteger.
 482      *
 483      * @param  numBits maximum bitLength of the new BigInteger.
 484      * @param  rnd source of randomness to be used in computing the new
 485      *         BigInteger.
 486      * @throws IllegalArgumentException {@code numBits} is negative.
 487      * @see #bitLength()
 488      */
 489     public BigInteger(int numBits, Random rnd) {
 490         this(1, randomBits(numBits, rnd));
 491     }
 492 
 493     private static byte[] randomBits(int numBits, Random rnd) {
 494         if (numBits < 0)
 495             throw new IllegalArgumentException("numBits must be non-negative");
 496         int numBytes = (int)(((long)numBits+7)/8); // avoid overflow
 497         byte[] randomBits = new byte[numBytes];
 498 
 499         // Generate random bytes and mask out any excess bits
 500         if (numBytes > 0) {
 501             rnd.nextBytes(randomBits);
 502             int excessBits = 8*numBytes - numBits;
 503             randomBits[0] &= (1 << (8-excessBits)) - 1;
 504         }
 505         return randomBits;
 506     }
 507 
 508     /**
 509      * Constructs a randomly generated positive BigInteger that is probably
 510      * prime, with the specified bitLength.
 511      *
 512      * <p>It is recommended that the {@link #probablePrime probablePrime}
 513      * method be used in preference to this constructor unless there
 514      * is a compelling need to specify a certainty.
 515      *
 516      * @param  bitLength bitLength of the returned BigInteger.
 517      * @param  certainty a measure of the uncertainty that the caller is
 518      *         willing to tolerate.  The probability that the new BigInteger
 519      *         represents a prime number will exceed
 520      *         (1 - 1/2<sup>{@code certainty}</sup>).  The execution time of
 521      *         this constructor is proportional to the value of this parameter.
 522      * @param  rnd source of random bits used to select candidates to be
 523      *         tested for primality.
 524      * @throws ArithmeticException {@code bitLength < 2}.
 525      * @see    #bitLength()
 526      */
 527     public BigInteger(int bitLength, int certainty, Random rnd) {
 528         BigInteger prime;
 529 
 530         if (bitLength < 2)
 531             throw new ArithmeticException("bitLength < 2");
 532         // The cutoff of 95 was chosen empirically for best performance
 533         prime = (bitLength < 95 ? smallPrime(bitLength, certainty, rnd)
 534                                 : largePrime(bitLength, certainty, rnd));
 535         signum = 1;
 536         mag = prime.mag;
 537     }
 538 
 539     // Minimum size in bits that the requested prime number has
 540     // before we use the large prime number generating algorithms
 541     private static final int SMALL_PRIME_THRESHOLD = 95;
 542 
 543     // Certainty required to meet the spec of probablePrime
 544     private static final int DEFAULT_PRIME_CERTAINTY = 100;
 545 
 546     /**
 547      * Returns a positive BigInteger that is probably prime, with the
 548      * specified bitLength. The probability that a BigInteger returned
 549      * by this method is composite does not exceed 2<sup>-100</sup>.
 550      *
 551      * @param  bitLength bitLength of the returned BigInteger.
 552      * @param  rnd source of random bits used to select candidates to be
 553      *         tested for primality.
 554      * @return a BigInteger of {@code bitLength} bits that is probably prime
 555      * @throws ArithmeticException {@code bitLength < 2}.
 556      * @see    #bitLength()
 557      * @since 1.4
 558      */
 559     public static BigInteger probablePrime(int bitLength, Random rnd) {
 560         if (bitLength < 2)
 561             throw new ArithmeticException("bitLength < 2");
 562 
 563         // The cutoff of 95 was chosen empirically for best performance
 564         return (bitLength < SMALL_PRIME_THRESHOLD ?
 565                 smallPrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd) :
 566                 largePrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd));
 567     }
 568 
 569     /**
 570      * Find a random number of the specified bitLength that is probably prime.
 571      * This method is used for smaller primes, its performance degrades on
 572      * larger bitlengths.
 573      *
 574      * This method assumes bitLength > 1.
 575      */
 576     private static BigInteger smallPrime(int bitLength, int certainty, Random rnd) {
 577         int magLen = (bitLength + 31) >>> 5;
 578         int temp[] = new int[magLen];
 579         int highBit = 1 << ((bitLength+31) & 0x1f);  // High bit of high int
 580         int highMask = (highBit << 1) - 1;  // Bits to keep in high int
 581 
 582         while(true) {
 583             // Construct a candidate
 584             for (int i=0; i<magLen; i++)
 585                 temp[i] = rnd.nextInt();
 586             temp[0] = (temp[0] & highMask) | highBit;  // Ensure exact length
 587             if (bitLength > 2)
 588                 temp[magLen-1] |= 1;  // Make odd if bitlen > 2
 589 
 590             BigInteger p = new BigInteger(temp, 1);
 591 
 592             // Do cheap "pre-test" if applicable
 593             if (bitLength > 6) {
 594                 long r = p.remainder(SMALL_PRIME_PRODUCT).longValue();
 595                 if ((r%3==0)  || (r%5==0)  || (r%7==0)  || (r%11==0) ||
 596                     (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) ||
 597                     (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0))
 598                     continue; // Candidate is composite; try another
 599             }
 600 
 601             // All candidates of bitLength 2 and 3 are prime by this point
 602             if (bitLength < 4)
 603                 return p;
 604 
 605             // Do expensive test if we survive pre-test (or it's inapplicable)
 606             if (p.primeToCertainty(certainty, rnd))
 607                 return p;
 608         }
 609     }
 610 
 611     private static final BigInteger SMALL_PRIME_PRODUCT
 612                        = valueOf(3L*5*7*11*13*17*19*23*29*31*37*41);
 613 
 614     /**
 615      * Find a random number of the specified bitLength that is probably prime.
 616      * This method is more appropriate for larger bitlengths since it uses
 617      * a sieve to eliminate most composites before using a more expensive
 618      * test.
 619      */
 620     private static BigInteger largePrime(int bitLength, int certainty, Random rnd) {
 621         BigInteger p;
 622         p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
 623         p.mag[p.mag.length-1] &= 0xfffffffe;
 624 
 625         // Use a sieve length likely to contain the next prime number
 626         int searchLen = (bitLength / 20) * 64;
 627         BitSieve searchSieve = new BitSieve(p, searchLen);
 628         BigInteger candidate = searchSieve.retrieve(p, certainty, rnd);
 629 
 630         while ((candidate == null) || (candidate.bitLength() != bitLength)) {
 631             p = p.add(BigInteger.valueOf(2*searchLen));
 632             if (p.bitLength() != bitLength)
 633                 p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
 634             p.mag[p.mag.length-1] &= 0xfffffffe;
 635             searchSieve = new BitSieve(p, searchLen);
 636             candidate = searchSieve.retrieve(p, certainty, rnd);
 637         }
 638         return candidate;
 639     }
 640 
 641    /**
 642     * Returns the first integer greater than this {@code BigInteger} that
 643     * is probably prime.  The probability that the number returned by this
 644     * method is composite does not exceed 2<sup>-100</sup>. This method will
 645     * never skip over a prime when searching: if it returns {@code p}, there
 646     * is no prime {@code q} such that {@code this < q < p}.
 647     *
 648     * @return the first integer greater than this {@code BigInteger} that
 649     *         is probably prime.
 650     * @throws ArithmeticException {@code this < 0}.
 651     * @since 1.5
 652     */
 653     public BigInteger nextProbablePrime() {
 654         if (this.signum < 0)
 655             throw new ArithmeticException("start < 0: " + this);
 656 
 657         // Handle trivial cases
 658         if ((this.signum == 0) || this.equals(ONE))
 659             return TWO;
 660 
 661         BigInteger result = this.add(ONE);
 662 
 663         // Fastpath for small numbers
 664         if (result.bitLength() < SMALL_PRIME_THRESHOLD) {
 665 
 666             // Ensure an odd number
 667             if (!result.testBit(0))
 668                 result = result.add(ONE);
 669 
 670             while(true) {
 671                 // Do cheap "pre-test" if applicable
 672                 if (result.bitLength() > 6) {
 673                     long r = result.remainder(SMALL_PRIME_PRODUCT).longValue();
 674                     if ((r%3==0)  || (r%5==0)  || (r%7==0)  || (r%11==0) ||
 675                         (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) ||
 676                         (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0)) {
 677                         result = result.add(TWO);
 678                         continue; // Candidate is composite; try another
 679                     }
 680                 }
 681 
 682                 // All candidates of bitLength 2 and 3 are prime by this point
 683                 if (result.bitLength() < 4)
 684                     return result;
 685 
 686                 // The expensive test
 687                 if (result.primeToCertainty(DEFAULT_PRIME_CERTAINTY, null))
 688                     return result;
 689 
 690                 result = result.add(TWO);
 691             }
 692         }
 693 
 694         // Start at previous even number
 695         if (result.testBit(0))
 696             result = result.subtract(ONE);
 697 
 698         // Looking for the next large prime
 699         int searchLen = (result.bitLength() / 20) * 64;
 700 
 701         while(true) {
 702            BitSieve searchSieve = new BitSieve(result, searchLen);
 703            BigInteger candidate = searchSieve.retrieve(result,
 704                                                  DEFAULT_PRIME_CERTAINTY, null);
 705            if (candidate != null)
 706                return candidate;
 707            result = result.add(BigInteger.valueOf(2 * searchLen));
 708         }
 709     }
 710 
 711     /**
 712      * Returns {@code true} if this BigInteger is probably prime,
 713      * {@code false} if it's definitely composite.
 714      *
 715      * This method assumes bitLength > 2.
 716      *
 717      * @param  certainty a measure of the uncertainty that the caller is
 718      *         willing to tolerate: if the call returns {@code true}
 719      *         the probability that this BigInteger is prime exceeds
 720      *         {@code (1 - 1/2<sup>certainty</sup>)}.  The execution time of
 721      *         this method is proportional to the value of this parameter.
 722      * @return {@code true} if this BigInteger is probably prime,
 723      *         {@code false} if it's definitely composite.
 724      */
 725     boolean primeToCertainty(int certainty, Random random) {
 726         int rounds = 0;
 727         int n = (Math.min(certainty, Integer.MAX_VALUE-1)+1)/2;
 728 
 729         // The relationship between the certainty and the number of rounds
 730         // we perform is given in the draft standard ANSI X9.80, "PRIME
 731         // NUMBER GENERATION, PRIMALITY TESTING, AND PRIMALITY CERTIFICATES".
 732         int sizeInBits = this.bitLength();
 733         if (sizeInBits < 100) {
 734             rounds = 50;
 735             rounds = n < rounds ? n : rounds;
 736             return passesMillerRabin(rounds, random);
 737         }
 738 
 739         if (sizeInBits < 256) {
 740             rounds = 27;
 741         } else if (sizeInBits < 512) {
 742             rounds = 15;
 743         } else if (sizeInBits < 768) {
 744             rounds = 8;
 745         } else if (sizeInBits < 1024) {
 746             rounds = 4;
 747         } else {
 748             rounds = 2;
 749         }
 750         rounds = n < rounds ? n : rounds;
 751 
 752         return passesMillerRabin(rounds, random) && passesLucasLehmer();
 753     }
 754 
 755     /**
 756      * Returns true iff this BigInteger is a Lucas-Lehmer probable prime.
 757      *
 758      * The following assumptions are made:
 759      * This BigInteger is a positive, odd number.
 760      */
 761     private boolean passesLucasLehmer() {
 762         BigInteger thisPlusOne = this.add(ONE);
 763 
 764         // Step 1
 765         int d = 5;
 766         while (jacobiSymbol(d, this) != -1) {
 767             // 5, -7, 9, -11, ...
 768             d = (d<0) ? Math.abs(d)+2 : -(d+2);
 769         }
 770 
 771         // Step 2
 772         BigInteger u = lucasLehmerSequence(d, thisPlusOne, this);
 773 
 774         // Step 3
 775         return u.mod(this).equals(ZERO);
 776     }
 777 
 778     /**
 779      * Computes Jacobi(p,n).
 780      * Assumes n positive, odd, n>=3.
 781      */
 782     private static int jacobiSymbol(int p, BigInteger n) {
 783         if (p == 0)
 784             return 0;
 785 
 786         // Algorithm and comments adapted from Colin Plumb's C library.
 787         int j = 1;
 788         int u = n.mag[n.mag.length-1];
 789 
 790         // Make p positive
 791         if (p < 0) {
 792             p = -p;
 793             int n8 = u & 7;
 794             if ((n8 == 3) || (n8 == 7))
 795                 j = -j; // 3 (011) or 7 (111) mod 8
 796         }
 797 
 798         // Get rid of factors of 2 in p
 799         while ((p & 3) == 0)
 800             p >>= 2;
 801         if ((p & 1) == 0) {
 802             p >>= 1;
 803             if (((u ^ (u>>1)) & 2) != 0)
 804                 j = -j; // 3 (011) or 5 (101) mod 8
 805         }
 806         if (p == 1)
 807             return j;
 808         // Then, apply quadratic reciprocity
 809         if ((p & u & 2) != 0)   // p = u = 3 (mod 4)?
 810             j = -j;
 811         // And reduce u mod p
 812         u = n.mod(BigInteger.valueOf(p)).intValue();
 813 
 814         // Now compute Jacobi(u,p), u < p
 815         while (u != 0) {
 816             while ((u & 3) == 0)
 817                 u >>= 2;
 818             if ((u & 1) == 0) {
 819                 u >>= 1;
 820                 if (((p ^ (p>>1)) & 2) != 0)
 821                     j = -j;     // 3 (011) or 5 (101) mod 8
 822             }
 823             if (u == 1)
 824                 return j;
 825             // Now both u and p are odd, so use quadratic reciprocity
 826             assert (u < p);
 827             int t = u; u = p; p = t;
 828             if ((u & p & 2) != 0) // u = p = 3 (mod 4)?
 829                 j = -j;
 830             // Now u >= p, so it can be reduced
 831             u %= p;
 832         }
 833         return 0;
 834     }
 835 
 836     private static BigInteger lucasLehmerSequence(int z, BigInteger k, BigInteger n) {
 837         BigInteger d = BigInteger.valueOf(z);
 838         BigInteger u = ONE; BigInteger u2;
 839         BigInteger v = ONE; BigInteger v2;
 840 
 841         for (int i=k.bitLength()-2; i>=0; i--) {
 842             u2 = u.multiply(v).mod(n);
 843 
 844             v2 = v.square().add(d.multiply(u.square())).mod(n);
 845             if (v2.testBit(0)) {
 846                 v2 = n.subtract(v2);
 847                 v2.signum = - v2.signum;
 848             }
 849             v2 = v2.shiftRight(1);
 850 
 851             u = u2; v = v2;
 852             if (k.testBit(i)) {
 853                 u2 = u.add(v).mod(n);
 854                 if (u2.testBit(0)) {
 855                     u2 = n.subtract(u2);
 856                     u2.signum = - u2.signum;
 857                 }
 858                 u2 = u2.shiftRight(1);
 859 
 860                 v2 = v.add(d.multiply(u)).mod(n);
 861                 if (v2.testBit(0)) {
 862                     v2 = n.subtract(v2);
 863                     v2.signum = - v2.signum;
 864                 }
 865                 v2 = v2.shiftRight(1);
 866 
 867                 u = u2; v = v2;
 868             }
 869         }
 870         return u;
 871     }
 872 
 873     private static volatile Random staticRandom;
 874 
 875     private static Random getSecureRandom() {
 876         if (staticRandom == null) {
 877             staticRandom = new java.security.SecureRandom();
 878         }
 879         return staticRandom;
 880     }
 881 
 882     /**
 883      * Returns true iff this BigInteger passes the specified number of
 884      * Miller-Rabin tests. This test is taken from the DSA spec (NIST FIPS
 885      * 186-2).
 886      *
 887      * The following assumptions are made:
 888      * This BigInteger is a positive, odd number greater than 2.
 889      * iterations<=50.
 890      */
 891     private boolean passesMillerRabin(int iterations, Random rnd) {
 892         // Find a and m such that m is odd and this == 1 + 2**a * m
 893         BigInteger thisMinusOne = this.subtract(ONE);
 894         BigInteger m = thisMinusOne;
 895         int a = m.getLowestSetBit();
 896         m = m.shiftRight(a);
 897 
 898         // Do the tests
 899         if (rnd == null) {
 900             rnd = getSecureRandom();
 901         }
 902         for (int i=0; i<iterations; i++) {
 903             // Generate a uniform random on (1, this)
 904             BigInteger b;
 905             do {
 906                 b = new BigInteger(this.bitLength(), rnd);
 907             } while (b.compareTo(ONE) <= 0 || b.compareTo(this) >= 0);
 908 
 909             int j = 0;
 910             BigInteger z = b.modPow(m, this);
 911             while(!((j==0 && z.equals(ONE)) || z.equals(thisMinusOne))) {
 912                 if (j>0 && z.equals(ONE) || ++j==a)
 913                     return false;
 914                 z = z.modPow(TWO, this);
 915             }
 916         }
 917         return true;
 918     }
 919 
 920     /**
 921      * This private constructor differs from its public cousin
 922      * with the arguments reversed in two ways: it assumes that its
 923      * arguments are correct, and it doesn't copy the magnitude array.
 924      */
 925     private BigInteger(int[] magnitude, int signum) {
 926         this.signum = (magnitude.length==0 ? 0 : signum);
 927         this.mag = magnitude;
 928     }
 929 
 930     /**
 931      * This private constructor is for internal use and assumes that its
 932      * arguments are correct.
 933      */
 934     private BigInteger(byte[] magnitude, int signum) {
 935         this.signum = (magnitude.length==0 ? 0 : signum);
 936         this.mag = stripLeadingZeroBytes(magnitude);
 937     }
 938 
 939     /**
 940      * This private constructor is for internal use in converting
 941      * from a MutableBigInteger object into a BigInteger.
 942      */
 943     BigInteger(MutableBigInteger val, int sign) {
 944         if (val.offset > 0 || val.value.length != val.intLen) {
 945             mag = new int[val.intLen];
 946             for(int i=0; i<val.intLen; i++)
 947                 mag[i] = val.value[val.offset+i];
 948         } else {
 949             mag = val.value;
 950         }
 951 
 952         this.signum = (val.intLen == 0) ? 0 : sign;
 953     }
 954 
 955     //Static Factory Methods
 956 
 957     /**
 958      * Returns a BigInteger whose value is equal to that of the
 959      * specified {@code long}.  This "static factory method" is
 960      * provided in preference to a ({@code long}) constructor
 961      * because it allows for reuse of frequently used BigIntegers.
 962      *
 963      * @param  val value of the BigInteger to return.
 964      * @return a BigInteger with the specified value.
 965      */
 966     public static BigInteger valueOf(long val) {
 967         // If -MAX_CONSTANT < val < MAX_CONSTANT, return stashed constant
 968         if (val == 0)
 969             return ZERO;
 970         if (val > 0 && val <= MAX_CONSTANT)
 971             return posConst[(int) val];
 972         else if (val < 0 && val >= -MAX_CONSTANT)
 973             return negConst[(int) -val];
 974 
 975         return new BigInteger(val);
 976     }
 977 
 978     /**
 979      * Constructs a BigInteger with the specified value, which may not be zero.
 980      */
 981     private BigInteger(long val) {
 982         if (val < 0) {
 983             signum = -1;
 984             val = -val;
 985         } else {
 986             signum = 1;
 987         }
 988 
 989         int highWord = (int)(val >>> 32);
 990         if (highWord==0) {
 991             mag = new int[1];
 992             mag[0] = (int)val;
 993         } else {
 994             mag = new int[2];
 995             mag[0] = highWord;
 996             mag[1] = (int)val;
 997         }
 998     }
 999 
1000     /**
1001      * Returns a BigInteger with the given two's complement representation.
1002      * Assumes that the input array will not be modified (the returned
1003      * BigInteger will reference the input array if feasible).
1004      */
1005     private static BigInteger valueOf(int val[]) {
1006         return (val[0]>0 ? new BigInteger(val, 1) : new BigInteger(val));
1007     }
1008 
1009     // Constants
1010 
1011     /**
1012      * Initialize static constant array when class is loaded.
1013      */
1014     private final static int MAX_CONSTANT = 16;
1015     private static BigInteger posConst[] = new BigInteger[MAX_CONSTANT+1];
1016     private static BigInteger negConst[] = new BigInteger[MAX_CONSTANT+1];
1017     static {
1018         for (int i = 1; i <= MAX_CONSTANT; i++) {
1019             int[] magnitude = new int[1];
1020             magnitude[0] = i;
1021             posConst[i] = new BigInteger(magnitude,  1);
1022             negConst[i] = new BigInteger(magnitude, -1);
1023         }
1024     }
1025 
1026     /**
1027      * The BigInteger constant zero.
1028      *
1029      * @since   1.2
1030      */
1031     public static final BigInteger ZERO = new BigInteger(new int[0], 0);
1032 
1033     /**
1034      * The BigInteger constant one.
1035      *
1036      * @since   1.2
1037      */
1038     public static final BigInteger ONE = valueOf(1);
1039 
1040     /**
1041      * The BigInteger constant two.  (Not exported.)
1042      */
1043     private static final BigInteger TWO = valueOf(2);
1044 
1045     /**
1046      * The BigInteger constant ten.
1047      *
1048      * @since   1.5
1049      */
1050     public static final BigInteger TEN = valueOf(10);
1051 
1052     // Arithmetic Operations
1053 
1054     /**
1055      * Returns a BigInteger whose value is {@code (this + val)}.
1056      *
1057      * @param  val value to be added to this BigInteger.
1058      * @return {@code this + val}
1059      */
1060     public BigInteger add(BigInteger val) {
1061         int[] resultMag;
1062         if (val.signum == 0)
1063             return this;
1064         if (signum == 0)
1065             return val;
1066         if (val.signum == signum)
1067             return new BigInteger(add(mag, val.mag), signum);
1068 
1069         int cmp = intArrayCmp(mag, val.mag);
1070         if (cmp==0)
1071             return ZERO;
1072         resultMag = (cmp>0 ? subtract(mag, val.mag)
1073                            : subtract(val.mag, mag));
1074         resultMag = trustedStripLeadingZeroInts(resultMag);
1075 
1076         return new BigInteger(resultMag, cmp*signum);
1077     }
1078 
1079     /**
1080      * Adds the contents of the int arrays x and y. This method allocates
1081      * a new int array to hold the answer and returns a reference to that
1082      * array.
1083      */
1084     private static int[] add(int[] x, int[] y) {
1085         // If x is shorter, swap the two arrays
1086         if (x.length < y.length) {
1087             int[] tmp = x;
1088             x = y;
1089             y = tmp;
1090         }
1091 
1092         int xIndex = x.length;
1093         int yIndex = y.length;
1094         int result[] = new int[xIndex];
1095         long sum = 0;
1096 
1097         // Add common parts of both numbers
1098         while(yIndex > 0) {
1099             sum = (x[--xIndex] & LONG_MASK) +
1100                   (y[--yIndex] & LONG_MASK) + (sum >>> 32);
1101             result[xIndex] = (int)sum;
1102         }
1103 
1104         // Copy remainder of longer number while carry propagation is required
1105         boolean carry = (sum >>> 32 != 0);
1106         while (xIndex > 0 && carry)
1107             carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
1108 
1109         // Copy remainder of longer number
1110         while (xIndex > 0)
1111             result[--xIndex] = x[xIndex];
1112 
1113         // Grow result if necessary
1114         if (carry) {
1115             int newLen = result.length + 1;
1116             int temp[] = new int[newLen];
1117             for (int i = 1; i<newLen; i++)
1118                 temp[i] = result[i-1];
1119             temp[0] = 0x01;
1120             result = temp;
1121         }
1122         return result;
1123     }
1124 
1125     /**
1126      * Returns a BigInteger whose value is {@code (this - val)}.
1127      *
1128      * @param  val value to be subtracted from this BigInteger.
1129      * @return {@code this - val}
1130      */
1131     public BigInteger subtract(BigInteger val) {
1132         int[] resultMag;
1133         if (val.signum == 0)
1134             return this;
1135         if (signum == 0)
1136             return val.negate();
1137         if (val.signum != signum)
1138             return new BigInteger(add(mag, val.mag), signum);
1139 
1140         int cmp = intArrayCmp(mag, val.mag);
1141         if (cmp==0)
1142             return ZERO;
1143         resultMag = (cmp>0 ? subtract(mag, val.mag)
1144                            : subtract(val.mag, mag));
1145         resultMag = trustedStripLeadingZeroInts(resultMag);
1146         return new BigInteger(resultMag, cmp*signum);
1147     }
1148 
1149     /**
1150      * Subtracts the contents of the second int arrays (little) from the
1151      * first (big).  The first int array (big) must represent a larger number
1152      * than the second.  This method allocates the space necessary to hold the
1153      * answer.
1154      */
1155     private static int[] subtract(int[] big, int[] little) {
1156         int bigIndex = big.length;
1157         int result[] = new int[bigIndex];
1158         int littleIndex = little.length;
1159         long difference = 0;
1160 
1161         // Subtract common parts of both numbers
1162         while(littleIndex > 0) {
1163             difference = (big[--bigIndex] & LONG_MASK) -
1164                          (little[--littleIndex] & LONG_MASK) +
1165                          (difference >> 32);
1166             result[bigIndex] = (int)difference;
1167         }
1168 
1169         // Subtract remainder of longer number while borrow propagates
1170         boolean borrow = (difference >> 32 != 0);
1171         while (bigIndex > 0 && borrow)
1172             borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
1173 
1174         // Copy remainder of longer number
1175         while (bigIndex > 0)
1176             result[--bigIndex] = big[bigIndex];
1177 
1178         return result;
1179     }
1180 
1181     /**
1182      * Returns a BigInteger whose value is {@code (this * val)}.
1183      *
1184      * @param  val value to be multiplied by this BigInteger.
1185      * @return {@code this * val}
1186      */
1187     public BigInteger multiply(BigInteger val) {
1188         if (val.signum == 0 || signum == 0)
1189             return ZERO;
1190 
1191         int[] result = multiplyToLen(mag, mag.length,
1192                                      val.mag, val.mag.length, null);
1193         result = trustedStripLeadingZeroInts(result);
1194         return new BigInteger(result, signum*val.signum);
1195     }
1196 
1197     /**
1198      * Multiplies int arrays x and y to the specified lengths and places
1199      * the result into z.
1200      */
1201     private int[] multiplyToLen(int[] x, int xlen, int[] y, int ylen, int[] z) {
1202         int xstart = xlen - 1;
1203         int ystart = ylen - 1;
1204 
1205         if (z == null || z.length < (xlen+ ylen))
1206             z = new int[xlen+ylen];
1207 
1208         long carry = 0;
1209         for (int j=ystart, k=ystart+1+xstart; j>=0; j--, k--) {
1210             long product = (y[j] & LONG_MASK) *
1211                            (x[xstart] & LONG_MASK) + carry;
1212             z[k] = (int)product;
1213             carry = product >>> 32;
1214         }
1215         z[xstart] = (int)carry;
1216 
1217         for (int i = xstart-1; i >= 0; i--) {
1218             carry = 0;
1219             for (int j=ystart, k=ystart+1+i; j>=0; j--, k--) {
1220                 long product = (y[j] & LONG_MASK) *
1221                                (x[i] & LONG_MASK) +
1222                                (z[k] & LONG_MASK) + carry;
1223                 z[k] = (int)product;
1224                 carry = product >>> 32;
1225             }
1226             z[i] = (int)carry;
1227         }
1228         return z;
1229     }
1230 
1231     /**
1232      * Returns a BigInteger whose value is {@code (this<sup>2</sup>)}.
1233      *
1234      * @return {@code this<sup>2</sup>}
1235      */
1236     private BigInteger square() {
1237         if (signum == 0)
1238             return ZERO;
1239         int[] z = squareToLen(mag, mag.length, null);
1240         return new BigInteger(trustedStripLeadingZeroInts(z), 1);
1241     }
1242 
1243     /**
1244      * Squares the contents of the int array x. The result is placed into the
1245      * int array z.  The contents of x are not changed.
1246      */
1247     private static final int[] squareToLen(int[] x, int len, int[] z) {
1248         /*
1249          * The algorithm used here is adapted from Colin Plumb's C library.
1250          * Technique: Consider the partial products in the multiplication
1251          * of "abcde" by itself:
1252          *
1253          *               a  b  c  d  e
1254          *            *  a  b  c  d  e
1255          *          ==================
1256          *              ae be ce de ee
1257          *           ad bd cd dd de
1258          *        ac bc cc cd ce
1259          *     ab bb bc bd be
1260          *  aa ab ac ad ae
1261          *
1262          * Note that everything above the main diagonal:
1263          *              ae be ce de = (abcd) * e
1264          *           ad bd cd       = (abc) * d
1265          *        ac bc             = (ab) * c
1266          *     ab                   = (a) * b
1267          *
1268          * is a copy of everything below the main diagonal:
1269          *                       de
1270          *                 cd ce
1271          *           bc bd be
1272          *     ab ac ad ae
1273          *
1274          * Thus, the sum is 2 * (off the diagonal) + diagonal.
1275          *
1276          * This is accumulated beginning with the diagonal (which
1277          * consist of the squares of the digits of the input), which is then
1278          * divided by two, the off-diagonal added, and multiplied by two
1279          * again.  The low bit is simply a copy of the low bit of the
1280          * input, so it doesn't need special care.
1281          */
1282         int zlen = len << 1;
1283         if (z == null || z.length < zlen)
1284             z = new int[zlen];
1285 
1286         // Store the squares, right shifted one bit (i.e., divided by 2)
1287         int lastProductLowWord = 0;
1288         for (int j=0, i=0; j<len; j++) {
1289             long piece = (x[j] & LONG_MASK);
1290             long product = piece * piece;
1291             z[i++] = (lastProductLowWord << 31) | (int)(product >>> 33);
1292             z[i++] = (int)(product >>> 1);
1293             lastProductLowWord = (int)product;
1294         }
1295 
1296         // Add in off-diagonal sums
1297         for (int i=len, offset=1; i>0; i--, offset+=2) {
1298             int t = x[i-1];
1299             t = mulAdd(z, x, offset, i-1, t);
1300             addOne(z, offset-1, i, t);
1301         }
1302 
1303         // Shift back up and set low bit
1304         primitiveLeftShift(z, zlen, 1);
1305         z[zlen-1] |= x[len-1] & 1;
1306 
1307         return z;
1308     }
1309 
1310     /**
1311      * Returns a BigInteger whose value is {@code (this / val)}.
1312      *
1313      * @param  val value by which this BigInteger is to be divided.
1314      * @return {@code this / val}
1315      * @throws ArithmeticException {@code val==0}
1316      */
1317     public BigInteger divide(BigInteger val) {
1318         MutableBigInteger q = new MutableBigInteger(),
1319                           r = new MutableBigInteger(),
1320                           a = new MutableBigInteger(this.mag),
1321                           b = new MutableBigInteger(val.mag);
1322 
1323         a.divide(b, q, r);
1324         return new BigInteger(q, this.signum * val.signum);
1325     }
1326 
1327     /**
1328      * Returns an array of two BigIntegers containing {@code (this / val)}
1329      * followed by {@code (this % val)}.
1330      *
1331      * @param  val value by which this BigInteger is to be divided, and the
1332      *         remainder computed.
1333      * @return an array of two BigIntegers: the quotient {@code (this / val)}
1334      *         is the initial element, and the remainder {@code (this % val)}
1335      *         is the final element.
1336      * @throws ArithmeticException {@code val==0}
1337      */
1338     public BigInteger[] divideAndRemainder(BigInteger val) {
1339         BigInteger[] result = new BigInteger[2];
1340         MutableBigInteger q = new MutableBigInteger(),
1341                           r = new MutableBigInteger(),
1342                           a = new MutableBigInteger(this.mag),
1343                           b = new MutableBigInteger(val.mag);
1344         a.divide(b, q, r);
1345         result[0] = new BigInteger(q, this.signum * val.signum);
1346         result[1] = new BigInteger(r, this.signum);
1347         return result;
1348     }
1349 
1350     /**
1351      * Returns a BigInteger whose value is {@code (this % val)}.
1352      *
1353      * @param  val value by which this BigInteger is to be divided, and the
1354      *         remainder computed.
1355      * @return {@code this % val}
1356      * @throws ArithmeticException {@code val==0}
1357      */
1358     public BigInteger remainder(BigInteger val) {
1359         MutableBigInteger q = new MutableBigInteger(),
1360                           r = new MutableBigInteger(),
1361                           a = new MutableBigInteger(this.mag),
1362                           b = new MutableBigInteger(val.mag);
1363 
1364         a.divide(b, q, r);
1365         return new BigInteger(r, this.signum);
1366     }
1367 
1368     /**
1369      * Returns a BigInteger whose value is <tt>(this<sup>exponent</sup>)</tt>.
1370      * Note that {@code exponent} is an integer rather than a BigInteger.
1371      *
1372      * @param  exponent exponent to which this BigInteger is to be raised.
1373      * @return <tt>this<sup>exponent</sup></tt>
1374      * @throws ArithmeticException {@code exponent} is negative.  (This would
1375      *         cause the operation to yield a non-integer value.)
1376      */
1377     public BigInteger pow(int exponent) {
1378         if (exponent < 0)
1379             throw new ArithmeticException("Negative exponent");
1380         if (signum==0)
1381             return (exponent==0 ? ONE : this);
1382 
1383         // Perform exponentiation using repeated squaring trick
1384         int newSign = (signum<0 && (exponent&1)==1 ? -1 : 1);
1385         int[] baseToPow2 = this.mag;
1386         int[] result = {1};
1387 
1388         while (exponent != 0) {
1389             if ((exponent & 1)==1) {
1390                 result = multiplyToLen(result, result.length,
1391                                        baseToPow2, baseToPow2.length, null);
1392                 result = trustedStripLeadingZeroInts(result);
1393             }
1394             if ((exponent >>>= 1) != 0) {
1395                 baseToPow2 = squareToLen(baseToPow2, baseToPow2.length, null);
1396                 baseToPow2 = trustedStripLeadingZeroInts(baseToPow2);
1397             }
1398         }
1399         return new BigInteger(result, newSign);
1400     }
1401 
1402     /**
1403      * Returns a BigInteger whose value is the greatest common divisor of
1404      * {@code abs(this)} and {@code abs(val)}.  Returns 0 if
1405      * {@code this==0 && val==0}.
1406      *
1407      * @param  val value with which the GCD is to be computed.
1408      * @return {@code GCD(abs(this), abs(val))}
1409      */
1410     public BigInteger gcd(BigInteger val) {
1411         if (val.signum == 0)
1412             return this.abs();
1413         else if (this.signum == 0)
1414             return val.abs();
1415 
1416         MutableBigInteger a = new MutableBigInteger(this);
1417         MutableBigInteger b = new MutableBigInteger(val);
1418 
1419         MutableBigInteger result = a.hybridGCD(b);
1420 
1421         return new BigInteger(result, 1);
1422     }
1423 
1424     /**
1425      * Left shift int array a up to len by n bits. Returns the array that
1426      * results from the shift since space may have to be reallocated.
1427      */
1428     private static int[] leftShift(int[] a, int len, int n) {
1429         int nInts = n >>> 5;
1430         int nBits = n&0x1F;
1431         int bitsInHighWord = bitLen(a[0]);
1432 
1433         // If shift can be done without recopy, do so
1434         if (n <= (32-bitsInHighWord)) {
1435             primitiveLeftShift(a, len, nBits);
1436             return a;
1437         } else { // Array must be resized
1438             if (nBits <= (32-bitsInHighWord)) {
1439                 int result[] = new int[nInts+len];
1440                 for (int i=0; i<len; i++)
1441                     result[i] = a[i];
1442                 primitiveLeftShift(result, result.length, nBits);
1443                 return result;
1444             } else {
1445                 int result[] = new int[nInts+len+1];
1446                 for (int i=0; i<len; i++)
1447                     result[i] = a[i];
1448                 primitiveRightShift(result, result.length, 32 - nBits);
1449                 return result;
1450             }
1451         }
1452     }
1453 
1454     // shifts a up to len right n bits assumes no leading zeros, 0<n<32
1455     static void primitiveRightShift(int[] a, int len, int n) {
1456         int n2 = 32 - n;
1457         for (int i=len-1, c=a[i]; i>0; i--) {
1458             int b = c;
1459             c = a[i-1];
1460             a[i] = (c << n2) | (b >>> n);
1461         }
1462         a[0] >>>= n;
1463     }
1464 
1465     // shifts a up to len left n bits assumes no leading zeros, 0<=n<32
1466     static void primitiveLeftShift(int[] a, int len, int n) {
1467         if (len == 0 || n == 0)
1468             return;
1469 
1470         int n2 = 32 - n;
1471         for (int i=0, c=a[i], m=i+len-1; i<m; i++) {
1472             int b = c;
1473             c = a[i+1];
1474             a[i] = (b << n) | (c >>> n2);
1475         }
1476         a[len-1] <<= n;
1477     }
1478 
1479     /**
1480      * Calculate bitlength of contents of the first len elements an int array,
1481      * assuming there are no leading zero ints.
1482      */
1483     private static int bitLength(int[] val, int len) {
1484         if (len==0)
1485             return 0;
1486         return ((len-1)<<5) + bitLen(val[0]);
1487     }
1488 
1489     /**
1490      * Returns a BigInteger whose value is the absolute value of this
1491      * BigInteger.
1492      *
1493      * @return {@code abs(this)}
1494      */
1495     public BigInteger abs() {
1496         return (signum >= 0 ? this : this.negate());
1497     }
1498 
1499     /**
1500      * Returns a BigInteger whose value is {@code (-this)}.
1501      *
1502      * @return {@code -this}
1503      */
1504     public BigInteger negate() {
1505         return new BigInteger(this.mag, -this.signum);
1506     }
1507 
1508     /**
1509      * Returns the signum function of this BigInteger.
1510      *
1511      * @return -1, 0 or 1 as the value of this BigInteger is negative, zero or
1512      *         positive.
1513      */
1514     public int signum() {
1515         return this.signum;
1516     }
1517 
1518     // Modular Arithmetic Operations
1519 
1520     /**
1521      * Returns a BigInteger whose value is {@code (this mod m}).  This method
1522      * differs from {@code remainder} in that it always returns a
1523      * <i>non-negative</i> BigInteger.
1524      *
1525      * @param  m the modulus.
1526      * @return {@code this mod m}
1527      * @throws ArithmeticException {@code m <= 0}
1528      * @see    #remainder
1529      */
1530     public BigInteger mod(BigInteger m) {
1531         if (m.signum <= 0)
1532             throw new ArithmeticException("BigInteger: modulus not positive");
1533 
1534         BigInteger result = this.remainder(m);
1535         return (result.signum >= 0 ? result : result.add(m));
1536     }
1537 
1538     /**
1539      * Returns a BigInteger whose value is
1540      * <tt>(this<sup>exponent</sup> mod m)</tt>.  (Unlike {@code pow}, this
1541      * method permits negative exponents.)
1542      *
1543      * @param  exponent the exponent.
1544      * @param  m the modulus.
1545      * @return <tt>this<sup>exponent</sup> mod m</tt>
1546      * @throws ArithmeticException {@code m <= 0}
1547      * @see    #modInverse
1548      */
1549     public BigInteger modPow(BigInteger exponent, BigInteger m) {
1550         if (m.signum <= 0)
1551             throw new ArithmeticException("BigInteger: modulus not positive");
1552 
1553         // Trivial cases
1554         if (exponent.signum == 0)
1555             return (m.equals(ONE) ? ZERO : ONE);
1556 
1557         if (this.equals(ONE))
1558             return (m.equals(ONE) ? ZERO : ONE);
1559 
1560         if (this.equals(ZERO) && exponent.signum >= 0)
1561             return ZERO;
1562 
1563         if (this.equals(negConst[1]) && (!exponent.testBit(0)))
1564             return (m.equals(ONE) ? ZERO : ONE);
1565 
1566         boolean invertResult;
1567         if ((invertResult = (exponent.signum < 0)))
1568             exponent = exponent.negate();
1569 
1570         BigInteger base = (this.signum < 0 || this.compareTo(m) >= 0
1571                            ? this.mod(m) : this);
1572         BigInteger result;
1573         if (m.testBit(0)) { // odd modulus
1574             result = base.oddModPow(exponent, m);
1575         } else {
1576             /*
1577              * Even modulus.  Tear it into an "odd part" (m1) and power of two
1578              * (m2), exponentiate mod m1, manually exponentiate mod m2, and
1579              * use Chinese Remainder Theorem to combine results.
1580              */
1581 
1582             // Tear m apart into odd part (m1) and power of 2 (m2)
1583             int p = m.getLowestSetBit();   // Max pow of 2 that divides m
1584 
1585             BigInteger m1 = m.shiftRight(p);  // m/2**p
1586             BigInteger m2 = ONE.shiftLeft(p); // 2**p
1587 
1588             // Calculate new base from m1
1589             BigInteger base2 = (this.signum < 0 || this.compareTo(m1) >= 0
1590                                 ? this.mod(m1) : this);
1591 
1592             // Caculate (base ** exponent) mod m1.
1593             BigInteger a1 = (m1.equals(ONE) ? ZERO :
1594                              base2.oddModPow(exponent, m1));
1595 
1596             // Calculate (this ** exponent) mod m2
1597             BigInteger a2 = base.modPow2(exponent, p);
1598 
1599             // Combine results using Chinese Remainder Theorem
1600             BigInteger y1 = m2.modInverse(m1);
1601             BigInteger y2 = m1.modInverse(m2);
1602 
1603             result = a1.multiply(m2).multiply(y1).add
1604                      (a2.multiply(m1).multiply(y2)).mod(m);
1605         }
1606 
1607         return (invertResult ? result.modInverse(m) : result);
1608     }
1609 
1610     static int[] bnExpModThreshTable = {7, 25, 81, 241, 673, 1793,
1611                                                 Integer.MAX_VALUE}; // Sentinel
1612 
1613     /**
1614      * Returns a BigInteger whose value is x to the power of y mod z.
1615      * Assumes: z is odd && x < z.
1616      */
1617     private BigInteger oddModPow(BigInteger y, BigInteger z) {
1618     /*
1619      * The algorithm is adapted from Colin Plumb's C library.
1620      *
1621      * The window algorithm:
1622      * The idea is to keep a running product of b1 = n^(high-order bits of exp)
1623      * and then keep appending exponent bits to it.  The following patterns
1624      * apply to a 3-bit window (k = 3):
1625      * To append   0: square
1626      * To append   1: square, multiply by n^1
1627      * To append  10: square, multiply by n^1, square
1628      * To append  11: square, square, multiply by n^3
1629      * To append 100: square, multiply by n^1, square, square
1630      * To append 101: square, square, square, multiply by n^5
1631      * To append 110: square, square, multiply by n^3, square
1632      * To append 111: square, square, square, multiply by n^7
1633      *
1634      * Since each pattern involves only one multiply, the longer the pattern
1635      * the better, except that a 0 (no multiplies) can be appended directly.
1636      * We precompute a table of odd powers of n, up to 2^k, and can then
1637      * multiply k bits of exponent at a time.  Actually, assuming random
1638      * exponents, there is on average one zero bit between needs to
1639      * multiply (1/2 of the time there's none, 1/4 of the time there's 1,
1640      * 1/8 of the time, there's 2, 1/32 of the time, there's 3, etc.), so
1641      * you have to do one multiply per k+1 bits of exponent.
1642      *
1643      * The loop walks down the exponent, squaring the result buffer as
1644      * it goes.  There is a wbits+1 bit lookahead buffer, buf, that is
1645      * filled with the upcoming exponent bits.  (What is read after the
1646      * end of the exponent is unimportant, but it is filled with zero here.)
1647      * When the most-significant bit of this buffer becomes set, i.e.
1648      * (buf & tblmask) != 0, we have to decide what pattern to multiply
1649      * by, and when to do it.  We decide, remember to do it in future
1650      * after a suitable number of squarings have passed (e.g. a pattern
1651      * of "100" in the buffer requires that we multiply by n^1 immediately;
1652      * a pattern of "110" calls for multiplying by n^3 after one more
1653      * squaring), clear the buffer, and continue.
1654      *
1655      * When we start, there is one more optimization: the result buffer
1656      * is implcitly one, so squaring it or multiplying by it can be
1657      * optimized away.  Further, if we start with a pattern like "100"
1658      * in the lookahead window, rather than placing n into the buffer
1659      * and then starting to square it, we have already computed n^2
1660      * to compute the odd-powers table, so we can place that into
1661      * the buffer and save a squaring.
1662      *
1663      * This means that if you have a k-bit window, to compute n^z,
1664      * where z is the high k bits of the exponent, 1/2 of the time
1665      * it requires no squarings.  1/4 of the time, it requires 1
1666      * squaring, ... 1/2^(k-1) of the time, it reqires k-2 squarings.
1667      * And the remaining 1/2^(k-1) of the time, the top k bits are a
1668      * 1 followed by k-1 0 bits, so it again only requires k-2
1669      * squarings, not k-1.  The average of these is 1.  Add that
1670      * to the one squaring we have to do to compute the table,
1671      * and you'll see that a k-bit window saves k-2 squarings
1672      * as well as reducing the multiplies.  (It actually doesn't
1673      * hurt in the case k = 1, either.)
1674      */
1675         // Special case for exponent of one
1676         if (y.equals(ONE))
1677             return this;
1678 
1679         // Special case for base of zero
1680         if (signum==0)
1681             return ZERO;
1682 
1683         int[] base = mag.clone();
1684         int[] exp = y.mag;
1685         int[] mod = z.mag;
1686         int modLen = mod.length;
1687 
1688         // Select an appropriate window size
1689         int wbits = 0;
1690         int ebits = bitLength(exp, exp.length);
1691         // if exponent is 65537 (0x10001), use minimum window size
1692         if ((ebits != 17) || (exp[0] != 65537)) {
1693             while (ebits > bnExpModThreshTable[wbits]) {
1694                 wbits++;
1695             }
1696         }
1697 
1698         // Calculate appropriate table size
1699         int tblmask = 1 << wbits;
1700 
1701         // Allocate table for precomputed odd powers of base in Montgomery form
1702         int[][] table = new int[tblmask][];
1703         for (int i=0; i<tblmask; i++)
1704             table[i] = new int[modLen];
1705 
1706         // Compute the modular inverse
1707         int inv = -MutableBigInteger.inverseMod32(mod[modLen-1]);
1708 
1709         // Convert base to Montgomery form
1710         int[] a = leftShift(base, base.length, modLen << 5);
1711 
1712         MutableBigInteger q = new MutableBigInteger(),
1713                           r = new MutableBigInteger(),
1714                           a2 = new MutableBigInteger(a),
1715                           b2 = new MutableBigInteger(mod);
1716 
1717         a2.divide(b2, q, r);
1718         table[0] = r.toIntArray();
1719 
1720         // Pad table[0] with leading zeros so its length is at least modLen
1721         if (table[0].length < modLen) {
1722            int offset = modLen - table[0].length;
1723            int[] t2 = new int[modLen];
1724            for (int i=0; i<table[0].length; i++)
1725                t2[i+offset] = table[0][i];
1726            table[0] = t2;
1727         }
1728 
1729         // Set b to the square of the base
1730         int[] b = squareToLen(table[0], modLen, null);
1731         b = montReduce(b, mod, modLen, inv);
1732 
1733         // Set t to high half of b
1734         int[] t = new int[modLen];
1735         for(int i=0; i<modLen; i++)
1736             t[i] = b[i];
1737 
1738         // Fill in the table with odd powers of the base
1739         for (int i=1; i<tblmask; i++) {
1740             int[] prod = multiplyToLen(t, modLen, table[i-1], modLen, null);
1741             table[i] = montReduce(prod, mod, modLen, inv);
1742         }
1743 
1744         // Pre load the window that slides over the exponent
1745         int bitpos = 1 << ((ebits-1) & (32-1));
1746 
1747         int buf = 0;
1748         int elen = exp.length;
1749         int eIndex = 0;
1750         for (int i = 0; i <= wbits; i++) {
1751             buf = (buf << 1) | (((exp[eIndex] & bitpos) != 0)?1:0);
1752             bitpos >>>= 1;
1753             if (bitpos == 0) {
1754                 eIndex++;
1755                 bitpos = 1 << (32-1);
1756                 elen--;
1757             }
1758         }
1759 
1760         int multpos = ebits;
1761 
1762         // The first iteration, which is hoisted out of the main loop
1763         ebits--;
1764         boolean isone = true;
1765 
1766         multpos = ebits - wbits;
1767         while ((buf & 1) == 0) {
1768             buf >>>= 1;
1769             multpos++;
1770         }
1771 
1772         int[] mult = table[buf >>> 1];
1773 
1774         buf = 0;
1775         if (multpos == ebits)
1776             isone = false;
1777 
1778         // The main loop
1779         while(true) {
1780             ebits--;
1781             // Advance the window
1782             buf <<= 1;
1783 
1784             if (elen != 0) {
1785                 buf |= ((exp[eIndex] & bitpos) != 0) ? 1 : 0;
1786                 bitpos >>>= 1;
1787                 if (bitpos == 0) {
1788                     eIndex++;
1789                     bitpos = 1 << (32-1);
1790                     elen--;
1791                 }
1792             }
1793 
1794             // Examine the window for pending multiplies
1795             if ((buf & tblmask) != 0) {
1796                 multpos = ebits - wbits;
1797                 while ((buf & 1) == 0) {
1798                     buf >>>= 1;
1799                     multpos++;
1800                 }
1801                 mult = table[buf >>> 1];
1802                 buf = 0;
1803             }
1804 
1805             // Perform multiply
1806             if (ebits == multpos) {
1807                 if (isone) {
1808                     b = mult.clone();
1809                     isone = false;
1810                 } else {
1811                     t = b;
1812                     a = multiplyToLen(t, modLen, mult, modLen, a);
1813                     a = montReduce(a, mod, modLen, inv);
1814                     t = a; a = b; b = t;
1815                 }
1816             }
1817 
1818             // Check if done
1819             if (ebits == 0)
1820                 break;
1821 
1822             // Square the input
1823             if (!isone) {
1824                 t = b;
1825                 a = squareToLen(t, modLen, a);
1826                 a = montReduce(a, mod, modLen, inv);
1827                 t = a; a = b; b = t;
1828             }
1829         }
1830 
1831         // Convert result out of Montgomery form and return
1832         int[] t2 = new int[2*modLen];
1833         for(int i=0; i<modLen; i++)
1834             t2[i+modLen] = b[i];
1835 
1836         b = montReduce(t2, mod, modLen, inv);
1837 
1838         t2 = new int[modLen];
1839         for(int i=0; i<modLen; i++)
1840             t2[i] = b[i];
1841 
1842         return new BigInteger(1, t2);
1843     }
1844 
1845     /**
1846      * Montgomery reduce n, modulo mod.  This reduces modulo mod and divides
1847      * by 2^(32*mlen). Adapted from Colin Plumb's C library.
1848      */
1849     private static int[] montReduce(int[] n, int[] mod, int mlen, int inv) {
1850         int c=0;
1851         int len = mlen;
1852         int offset=0;
1853 
1854         do {
1855             int nEnd = n[n.length-1-offset];
1856             int carry = mulAdd(n, mod, offset, mlen, inv * nEnd);
1857             c += addOne(n, offset, mlen, carry);
1858             offset++;
1859         } while(--len > 0);
1860 
1861         while(c>0)
1862             c += subN(n, mod, mlen);
1863 
1864         while (intArrayCmpToLen(n, mod, mlen) >= 0)
1865             subN(n, mod, mlen);
1866 
1867         return n;
1868     }
1869 
1870 
1871     /*
1872      * Returns -1, 0 or +1 as big-endian unsigned int array arg1 is less than,
1873      * equal to, or greater than arg2 up to length len.
1874      */
1875     private static int intArrayCmpToLen(int[] arg1, int[] arg2, int len) {
1876         for (int i=0; i<len; i++) {
1877             long b1 = arg1[i] & LONG_MASK;
1878             long b2 = arg2[i] & LONG_MASK;
1879             if (b1 < b2)
1880                 return -1;
1881             if (b1 > b2)
1882                 return 1;
1883         }
1884         return 0;
1885     }
1886 
1887     /**
1888      * Subtracts two numbers of same length, returning borrow.
1889      */
1890     private static int subN(int[] a, int[] b, int len) {
1891         long sum = 0;
1892 
1893         while(--len >= 0) {
1894             sum = (a[len] & LONG_MASK) -
1895                  (b[len] & LONG_MASK) + (sum >> 32);
1896             a[len] = (int)sum;
1897         }
1898 
1899         return (int)(sum >> 32);
1900     }
1901 
1902     /**
1903      * Multiply an array by one word k and add to result, return the carry
1904      */
1905     static int mulAdd(int[] out, int[] in, int offset, int len, int k) {
1906         long kLong = k & LONG_MASK;
1907         long carry = 0;
1908 
1909         offset = out.length-offset - 1;
1910         for (int j=len-1; j >= 0; j--) {
1911             long product = (in[j] & LONG_MASK) * kLong +
1912                            (out[offset] & LONG_MASK) + carry;
1913             out[offset--] = (int)product;
1914             carry = product >>> 32;
1915         }
1916         return (int)carry;
1917     }
1918 
1919     /**
1920      * Add one word to the number a mlen words into a. Return the resulting
1921      * carry.
1922      */
1923     static int addOne(int[] a, int offset, int mlen, int carry) {
1924         offset = a.length-1-mlen-offset;
1925         long t = (a[offset] & LONG_MASK) + (carry & LONG_MASK);
1926 
1927         a[offset] = (int)t;
1928         if ((t >>> 32) == 0)
1929             return 0;
1930         while (--mlen >= 0) {
1931             if (--offset < 0) { // Carry out of number
1932                 return 1;
1933             } else {
1934                 a[offset]++;
1935                 if (a[offset] != 0)
1936                     return 0;
1937             }
1938         }
1939         return 1;
1940     }
1941 
1942     /**
1943      * Returns a BigInteger whose value is (this ** exponent) mod (2**p)
1944      */
1945     private BigInteger modPow2(BigInteger exponent, int p) {
1946         /*
1947          * Perform exponentiation using repeated squaring trick, chopping off
1948          * high order bits as indicated by modulus.
1949          */
1950         BigInteger result = valueOf(1);
1951         BigInteger baseToPow2 = this.mod2(p);
1952         int expOffset = 0;
1953 
1954         int limit = exponent.bitLength();
1955 
1956         if (this.testBit(0))
1957            limit = (p-1) < limit ? (p-1) : limit;
1958 
1959         while (expOffset < limit) {
1960             if (exponent.testBit(expOffset))
1961                 result = result.multiply(baseToPow2).mod2(p);
1962             expOffset++;
1963             if (expOffset < limit)
1964                 baseToPow2 = baseToPow2.square().mod2(p);
1965         }
1966 
1967         return result;
1968     }
1969 
1970     /**
1971      * Returns a BigInteger whose value is this mod(2**p).
1972      * Assumes that this {@code BigInteger >= 0} and {@code p > 0}.
1973      */
1974     private BigInteger mod2(int p) {
1975         if (bitLength() <= p)
1976             return this;
1977 
1978         // Copy remaining ints of mag
1979         int numInts = (p+31)/32;
1980         int[] mag = new int[numInts];
1981         for (int i=0; i<numInts; i++)
1982             mag[i] = this.mag[i + (this.mag.length - numInts)];
1983 
1984         // Mask out any excess bits
1985         int excessBits = (numInts << 5) - p;
1986         mag[0] &= (1L << (32-excessBits)) - 1;
1987 
1988         return (mag[0]==0 ? new BigInteger(1, mag) : new BigInteger(mag, 1));
1989     }
1990 
1991     /**
1992      * Returns a BigInteger whose value is {@code (this}<sup>-1</sup> {@code mod m)}.
1993      *
1994      * @param  m the modulus.
1995      * @return {@code this}<sup>-1</sup> {@code mod m}.
1996      * @throws ArithmeticException {@code  m <= 0}, or this BigInteger
1997      *         has no multiplicative inverse mod m (that is, this BigInteger
1998      *         is not <i>relatively prime</i> to m).
1999      */
2000     public BigInteger modInverse(BigInteger m) {
2001         if (m.signum != 1)
2002             throw new ArithmeticException("BigInteger: modulus not positive");
2003 
2004         if (m.equals(ONE))
2005             return ZERO;
2006 
2007         // Calculate (this mod m)
2008         BigInteger modVal = this;
2009         if (signum < 0 || (intArrayCmp(mag, m.mag) >= 0))
2010             modVal = this.mod(m);
2011 
2012         if (modVal.equals(ONE))
2013             return ONE;
2014 
2015         MutableBigInteger a = new MutableBigInteger(modVal);
2016         MutableBigInteger b = new MutableBigInteger(m);
2017 
2018         MutableBigInteger result = a.mutableModInverse(b);
2019         return new BigInteger(result, 1);
2020     }
2021 
2022     // Shift Operations
2023 
2024     /**
2025      * Returns a BigInteger whose value is {@code (this << n)}.
2026      * The shift distance, {@code n}, may be negative, in which case
2027      * this method performs a right shift.
2028      * (Computes <tt>floor(this * 2<sup>n</sup>)</tt>.)
2029      *
2030      * @param  n shift distance, in bits.
2031      * @return {@code this << n}
2032      * @see #shiftRight
2033      */
2034     public BigInteger shiftLeft(int n) {
2035         if (signum == 0)
2036             return ZERO;
2037         if (n==0)
2038             return this;
2039         if (n<0)
2040             return shiftRight(-n);
2041 
2042         int nInts = n >>> 5;
2043         int nBits = n & 0x1f;
2044         int magLen = mag.length;
2045         int newMag[] = null;
2046 
2047         if (nBits == 0) {
2048             newMag = new int[magLen + nInts];
2049             for (int i=0; i<magLen; i++)
2050                 newMag[i] = mag[i];
2051         } else {
2052             int i = 0;
2053             int nBits2 = 32 - nBits;
2054             int highBits = mag[0] >>> nBits2;
2055             if (highBits != 0) {
2056                 newMag = new int[magLen + nInts + 1];
2057                 newMag[i++] = highBits;
2058             } else {
2059                 newMag = new int[magLen + nInts];
2060             }
2061             int j=0;
2062             while (j < magLen-1)
2063                 newMag[i++] = mag[j++] << nBits | mag[j] >>> nBits2;
2064             newMag[i] = mag[j] << nBits;
2065         }
2066 
2067         return new BigInteger(newMag, signum);
2068     }
2069 
2070     /**
2071      * Returns a BigInteger whose value is {@code (this >> n)}.  Sign
2072      * extension is performed.  The shift distance, {@code n}, may be
2073      * negative, in which case this method performs a left shift.
2074      * (Computes <tt>floor(this / 2<sup>n</sup>)</tt>.)
2075      *
2076      * @param  n shift distance, in bits.
2077      * @return {@code this >> n}
2078      * @see #shiftLeft
2079      */
2080     public BigInteger shiftRight(int n) {
2081         if (n==0)
2082             return this;
2083         if (n<0)
2084             return shiftLeft(-n);
2085 
2086         int nInts = n >>> 5;
2087         int nBits = n & 0x1f;
2088         int magLen = mag.length;
2089         int newMag[] = null;
2090 
2091         // Special case: entire contents shifted off the end
2092         if (nInts >= magLen)
2093             return (signum >= 0 ? ZERO : negConst[1]);
2094 
2095         if (nBits == 0) {
2096             int newMagLen = magLen - nInts;
2097             newMag = new int[newMagLen];
2098             for (int i=0; i<newMagLen; i++)
2099                 newMag[i] = mag[i];
2100         } else {
2101             int i = 0;
2102             int highBits = mag[0] >>> nBits;
2103             if (highBits != 0) {
2104                 newMag = new int[magLen - nInts];
2105                 newMag[i++] = highBits;
2106             } else {
2107                 newMag = new int[magLen - nInts -1];
2108             }
2109 
2110             int nBits2 = 32 - nBits;
2111             int j=0;
2112             while (j < magLen - nInts - 1)
2113                 newMag[i++] = (mag[j++] << nBits2) | (mag[j] >>> nBits);
2114         }
2115 
2116         if (signum < 0) {
2117             // Find out whether any one-bits were shifted off the end.
2118             boolean onesLost = false;
2119             for (int i=magLen-1, j=magLen-nInts; i>=j && !onesLost; i--)
2120                 onesLost = (mag[i] != 0);
2121             if (!onesLost && nBits != 0)
2122                 onesLost = (mag[magLen - nInts - 1] << (32 - nBits) != 0);
2123 
2124             if (onesLost)
2125                 newMag = javaIncrement(newMag);
2126         }
2127 
2128         return new BigInteger(newMag, signum);
2129     }
2130 
2131     int[] javaIncrement(int[] val) {
2132         int lastSum = 0;
2133         for (int i=val.length-1;  i >= 0 && lastSum == 0; i--)
2134             lastSum = (val[i] += 1);
2135         if (lastSum == 0) {
2136             val = new int[val.length+1];
2137             val[0] = 1;
2138         }
2139         return val;
2140     }
2141 
2142     // Bitwise Operations
2143 
2144     /**
2145      * Returns a BigInteger whose value is {@code (this & val)}.  (This
2146      * method returns a negative BigInteger if and only if this and val are
2147      * both negative.)
2148      *
2149      * @param val value to be AND'ed with this BigInteger.
2150      * @return {@code this & val}
2151      */
2152     public BigInteger and(BigInteger val) {
2153         int[] result = new int[Math.max(intLength(), val.intLength())];
2154         for (int i=0; i<result.length; i++)
2155             result[i] = (getInt(result.length-i-1)
2156                          & val.getInt(result.length-i-1));
2157 
2158         return valueOf(result);
2159     }
2160 
2161     /**
2162      * Returns a BigInteger whose value is {@code (this | val)}.  (This method
2163      * returns a negative BigInteger if and only if either this or val is
2164      * negative.)
2165      *
2166      * @param val value to be OR'ed with this BigInteger.
2167      * @return {@code this | val}
2168      */
2169     public BigInteger or(BigInteger val) {
2170         int[] result = new int[Math.max(intLength(), val.intLength())];
2171         for (int i=0; i<result.length; i++)
2172             result[i] = (getInt(result.length-i-1)
2173                          | val.getInt(result.length-i-1));
2174 
2175         return valueOf(result);
2176     }
2177 
2178     /**
2179      * Returns a BigInteger whose value is {@code (this ^ val)}.  (This method
2180      * returns a negative BigInteger if and only if exactly one of this and
2181      * val are negative.)
2182      *
2183      * @param val value to be XOR'ed with this BigInteger.
2184      * @return {@code this ^ val}
2185      */
2186     public BigInteger xor(BigInteger val) {
2187         int[] result = new int[Math.max(intLength(), val.intLength())];
2188         for (int i=0; i<result.length; i++)
2189             result[i] = (getInt(result.length-i-1)
2190                          ^ val.getInt(result.length-i-1));
2191 
2192         return valueOf(result);
2193     }
2194 
2195     /**
2196      * Returns a BigInteger whose value is {@code (~this)}.  (This method
2197      * returns a negative value if and only if this BigInteger is
2198      * non-negative.)
2199      *
2200      * @return {@code ~this}
2201      */
2202     public BigInteger not() {
2203         int[] result = new int[intLength()];
2204         for (int i=0; i<result.length; i++)
2205             result[i] = ~getInt(result.length-i-1);
2206 
2207         return valueOf(result);
2208     }
2209 
2210     /**
2211      * Returns a BigInteger whose value is {@code (this & ~val)}.  This
2212      * method, which is equivalent to {@code and(val.not())}, is provided as
2213      * a convenience for masking operations.  (This method returns a negative
2214      * BigInteger if and only if {@code this} is negative and {@code val} is
2215      * positive.)
2216      *
2217      * @param val value to be complemented and AND'ed with this BigInteger.
2218      * @return {@code this & ~val}
2219      */
2220     public BigInteger andNot(BigInteger val) {
2221         int[] result = new int[Math.max(intLength(), val.intLength())];
2222         for (int i=0; i<result.length; i++)
2223             result[i] = (getInt(result.length-i-1)
2224                          & ~val.getInt(result.length-i-1));
2225 
2226         return valueOf(result);
2227     }
2228 
2229 
2230     // Single Bit Operations
2231 
2232     /**
2233      * Returns {@code true} if and only if the designated bit is set.
2234      * (Computes {@code ((this & (1<<n)) != 0)}.)
2235      *
2236      * @param  n index of bit to test.
2237      * @return {@code true} if and only if the designated bit is set.
2238      * @throws ArithmeticException {@code n} is negative.
2239      */
2240     public boolean testBit(int n) {
2241         if (n<0)
2242             throw new ArithmeticException("Negative bit address");
2243 
2244         return (getInt(n/32) & (1 << (n%32))) != 0;
2245     }
2246 
2247     /**
2248      * Returns a BigInteger whose value is equivalent to this BigInteger
2249      * with the designated bit set.  (Computes {@code (this | (1<<n))}.)
2250      *
2251      * @param  n index of bit to set.
2252      * @return {@code this | (1<<n)}
2253      * @throws ArithmeticException {@code n} is negative.
2254      */
2255     public BigInteger setBit(int n) {
2256         if (n<0)
2257             throw new ArithmeticException("Negative bit address");
2258 
2259         int intNum = n/32;
2260         int[] result = new int[Math.max(intLength(), intNum+2)];
2261 
2262         for (int i=0; i<result.length; i++)
2263             result[result.length-i-1] = getInt(i);
2264 
2265         result[result.length-intNum-1] |= (1 << (n%32));
2266 
2267         return valueOf(result);
2268     }
2269 
2270     /**
2271      * Returns a BigInteger whose value is equivalent to this BigInteger
2272      * with the designated bit cleared.
2273      * (Computes {@code (this & ~(1<<n))}.)
2274      *
2275      * @param  n index of bit to clear.
2276      * @return {@code this & ~(1<<n)}
2277      * @throws ArithmeticException {@code n} is negative.
2278      */
2279     public BigInteger clearBit(int n) {
2280         if (n<0)
2281             throw new ArithmeticException("Negative bit address");
2282 
2283         int intNum = n/32;
2284         int[] result = new int[Math.max(intLength(), (n+1)/32+1)];
2285 
2286         for (int i=0; i<result.length; i++)
2287             result[result.length-i-1] = getInt(i);
2288 
2289         result[result.length-intNum-1] &= ~(1 << (n%32));
2290 
2291         return valueOf(result);
2292     }
2293 
2294     /**
2295      * Returns a BigInteger whose value is equivalent to this BigInteger
2296      * with the designated bit flipped.
2297      * (Computes {@code (this ^ (1<<n))}.)
2298      *
2299      * @param  n index of bit to flip.
2300      * @return {@code this ^ (1<<n)}
2301      * @throws ArithmeticException {@code n} is negative.
2302      */
2303     public BigInteger flipBit(int n) {
2304         if (n<0)
2305             throw new ArithmeticException("Negative bit address");
2306 
2307         int intNum = n/32;
2308         int[] result = new int[Math.max(intLength(), intNum+2)];
2309 
2310         for (int i=0; i<result.length; i++)
2311             result[result.length-i-1] = getInt(i);
2312 
2313         result[result.length-intNum-1] ^= (1 << (n%32));
2314 
2315         return valueOf(result);
2316     }
2317 
2318     /**
2319      * Returns the index of the rightmost (lowest-order) one bit in this
2320      * BigInteger (the number of zero bits to the right of the rightmost
2321      * one bit).  Returns -1 if this BigInteger contains no one bits.
2322      * (Computes {@code (this==0? -1 : log2(this & -this))}.)
2323      *
2324      * @return index of the rightmost one bit in this BigInteger.
2325      */
2326     public int getLowestSetBit() {
2327         /*
2328          * Initialize lowestSetBit field the first time this method is
2329          * executed. This method depends on the atomicity of int modifies;
2330          * without this guarantee, it would have to be synchronized.
2331          */
2332         if (lowestSetBit == -2) {
2333             if (signum == 0) {
2334                 lowestSetBit = -1;
2335             } else {
2336                 // Search for lowest order nonzero int
2337                 int i,b;
2338                 for (i=0; (b = getInt(i))==0; i++)
2339                     ;
2340                 lowestSetBit = (i << 5) + trailingZeroCnt(b);
2341             }
2342         }
2343         return lowestSetBit;
2344     }
2345 
2346 
2347     // Miscellaneous Bit Operations
2348 
2349     /**
2350      * Returns the number of bits in the minimal two's-complement
2351      * representation of this BigInteger, <i>excluding</i> a sign bit.
2352      * For positive BigIntegers, this is equivalent to the number of bits in
2353      * the ordinary binary representation.  (Computes
2354      * {@code (ceil(log2(this < 0 ? -this : this+1)))}.)
2355      *
2356      * @return number of bits in the minimal two's-complement
2357      *         representation of this BigInteger, <i>excluding</i> a sign bit.
2358      */
2359     public int bitLength() {
2360         /*
2361          * Initialize bitLength field the first time this method is executed.
2362          * This method depends on the atomicity of int modifies; without
2363          * this guarantee, it would have to be synchronized.
2364          */
2365         if (bitLength == -1) {
2366             if (signum == 0) {
2367                 bitLength = 0;
2368             } else {
2369                 // Calculate the bit length of the magnitude
2370                 int magBitLength = ((mag.length-1) << 5) + bitLen(mag[0]);
2371 
2372                 if (signum < 0) {
2373                     // Check if magnitude is a power of two
2374                     boolean pow2 = (bitCnt(mag[0]) == 1);
2375                     for(int i=1; i<mag.length && pow2; i++)
2376                         pow2 = (mag[i]==0);
2377 
2378                     bitLength = (pow2 ? magBitLength-1 : magBitLength);
2379                 } else {
2380                     bitLength = magBitLength;
2381                 }
2382             }
2383         }
2384         return bitLength;
2385     }
2386 
2387     /**
2388      * bitLen(val) is the number of bits in val.
2389      */
2390     static int bitLen(int w) {
2391         // Binary search - decision tree (5 tests, rarely 6)
2392         return
2393          (w < 1<<15 ?
2394           (w < 1<<7 ?
2395            (w < 1<<3 ?
2396             (w < 1<<1 ? (w < 1<<0 ? (w<0 ? 32 : 0) : 1) : (w < 1<<2 ? 2 : 3)) :
2397             (w < 1<<5 ? (w < 1<<4 ? 4 : 5) : (w < 1<<6 ? 6 : 7))) :
2398            (w < 1<<11 ?
2399             (w < 1<<9 ? (w < 1<<8 ? 8 : 9) : (w < 1<<10 ? 10 : 11)) :
2400             (w < 1<<13 ? (w < 1<<12 ? 12 : 13) : (w < 1<<14 ? 14 : 15)))) :
2401           (w < 1<<23 ?
2402            (w < 1<<19 ?
2403             (w < 1<<17 ? (w < 1<<16 ? 16 : 17) : (w < 1<<18 ? 18 : 19)) :
2404             (w < 1<<21 ? (w < 1<<20 ? 20 : 21) : (w < 1<<22 ? 22 : 23))) :
2405            (w < 1<<27 ?
2406             (w < 1<<25 ? (w < 1<<24 ? 24 : 25) : (w < 1<<26 ? 26 : 27)) :
2407             (w < 1<<29 ? (w < 1<<28 ? 28 : 29) : (w < 1<<30 ? 30 : 31)))));
2408     }
2409 
2410     /*
2411      * trailingZeroTable[i] is the number of trailing zero bits in the binary
2412      * representation of i.
2413      */
2414     final static byte trailingZeroTable[] = {
2415       -25, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2416         4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2417         5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2418         4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2419         6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2420         4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2421         5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2422         4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2423         7, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2424         4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2425         5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2426         4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2427         6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2428         4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2429         5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
2430         4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0};
2431 
2432     /**
2433      * Returns the number of bits in the two's complement representation
2434      * of this BigInteger that differ from its sign bit.  This method is
2435      * useful when implementing bit-vector style sets atop BigIntegers.
2436      *
2437      * @return number of bits in the two's complement representation
2438      *         of this BigInteger that differ from its sign bit.
2439      */
2440     public int bitCount() {
2441         /*
2442          * Initialize bitCount field the first time this method is executed.
2443          * This method depends on the atomicity of int modifies; without
2444          * this guarantee, it would have to be synchronized.
2445          */
2446         if (bitCount == -1) {
2447             // Count the bits in the magnitude
2448             int magBitCount = 0;
2449             for (int i=0; i<mag.length; i++)
2450                 magBitCount += bitCnt(mag[i]);
2451 
2452             if (signum < 0) {
2453                 // Count the trailing zeros in the magnitude
2454                 int magTrailingZeroCount = 0, j;
2455                 for (j=mag.length-1; mag[j]==0; j--)
2456                     magTrailingZeroCount += 32;
2457                 magTrailingZeroCount +=
2458                             trailingZeroCnt(mag[j]);
2459 
2460                 bitCount = magBitCount + magTrailingZeroCount - 1;
2461             } else {
2462                 bitCount = magBitCount;
2463             }
2464         }
2465         return bitCount;
2466     }
2467 
2468     static int bitCnt(int val) {
2469         val -= (0xaaaaaaaa & val) >>> 1;
2470         val = (val & 0x33333333) + ((val >>> 2) & 0x33333333);
2471         val = val + (val >>> 4) & 0x0f0f0f0f;
2472         val += val >>> 8;
2473         val += val >>> 16;
2474         return val & 0xff;
2475     }
2476 
2477     static int trailingZeroCnt(int val) {
2478         // Loop unrolled for performance
2479         int byteVal = val & 0xff;
2480         if (byteVal != 0)
2481             return trailingZeroTable[byteVal];
2482 
2483         byteVal = (val >>> 8) & 0xff;
2484         if (byteVal != 0)
2485             return trailingZeroTable[byteVal] + 8;
2486 
2487         byteVal = (val >>> 16) & 0xff;
2488         if (byteVal != 0)
2489             return trailingZeroTable[byteVal] + 16;
2490 
2491         byteVal = (val >>> 24) & 0xff;
2492         return trailingZeroTable[byteVal] + 24;
2493     }
2494 
2495     // Primality Testing
2496 
2497     /**
2498      * Returns {@code true} if this BigInteger is probably prime,
2499      * {@code false} if it's definitely composite.  If
2500      * {@code certainty} is {@code  <= 0}, {@code true} is
2501      * returned.
2502      *
2503      * @param  certainty a measure of the uncertainty that the caller is
2504      *         willing to tolerate: if the call returns {@code true}
2505      *         the probability that this BigInteger is prime exceeds
2506      *         (1 - 1/2<sup>{@code certainty}</sup>).  The execution time of
2507      *         this method is proportional to the value of this parameter.
2508      * @return {@code true} if this BigInteger is probably prime,
2509      *         {@code false} if it's definitely composite.
2510      */
2511     public boolean isProbablePrime(int certainty) {
2512         if (certainty <= 0)
2513             return true;
2514         BigInteger w = this.abs();
2515         if (w.equals(TWO))
2516             return true;
2517         if (!w.testBit(0) || w.equals(ONE))
2518             return false;
2519 
2520         return w.primeToCertainty(certainty, null);
2521     }
2522 
2523     // Comparison Operations
2524 
2525     /**
2526      * Compares this BigInteger with the specified BigInteger.  This
2527      * method is provided in preference to individual methods for each
2528      * of the six boolean comparison operators ({@literal <}, ==,
2529      * {@literal >}, {@literal >=}, !=, {@literal <=}).  The suggested
2530      * idiom for performing these comparisons is: {@code
2531      * (x.compareTo(y)} &lt;<i>op</i>&gt; {@code 0)}, where
2532      * &lt;<i>op</i>&gt; is one of the six comparison operators.
2533      *
2534      * @param  val BigInteger to which this BigInteger is to be compared.
2535      * @return -1, 0 or 1 as this BigInteger is numerically less than, equal
2536      *         to, or greater than {@code val}.
2537      */
2538     public int compareTo(BigInteger val) {
2539         return (signum==val.signum
2540                 ? signum*intArrayCmp(mag, val.mag)
2541                 : (signum>val.signum ? 1 : -1));
2542     }
2543 
2544     /*
2545      * Returns -1, 0 or +1 as big-endian unsigned int array arg1 is
2546      * less than, equal to, or greater than arg2.
2547      */
2548     private static int intArrayCmp(int[] arg1, int[] arg2) {
2549         if (arg1.length < arg2.length)
2550             return -1;
2551         if (arg1.length > arg2.length)
2552             return 1;
2553 
2554         // Argument lengths are equal; compare the values
2555         for (int i=0; i<arg1.length; i++) {
2556             long b1 = arg1[i] & LONG_MASK;
2557             long b2 = arg2[i] & LONG_MASK;
2558             if (b1 < b2)
2559                 return -1;
2560             if (b1 > b2)
2561                 return 1;
2562         }
2563         return 0;
2564     }
2565 
2566     /**
2567      * Compares this BigInteger with the specified Object for equality.
2568      *
2569      * @param  x Object to which this BigInteger is to be compared.
2570      * @return {@code true} if and only if the specified Object is a
2571      *         BigInteger whose value is numerically equal to this BigInteger.
2572      */
2573     public boolean equals(Object x) {
2574         // This test is just an optimization, which may or may not help
2575         if (x == this)
2576             return true;
2577 
2578         if (!(x instanceof BigInteger))
2579             return false;
2580         BigInteger xInt = (BigInteger) x;
2581 
2582         if (xInt.signum != signum || xInt.mag.length != mag.length)
2583             return false;
2584 
2585         for (int i=0; i<mag.length; i++)
2586             if (xInt.mag[i] != mag[i])
2587                 return false;
2588 
2589         return true;
2590     }
2591 
2592     /**
2593      * Returns the minimum of this BigInteger and {@code val}.
2594      *
2595      * @param  val value with which the minimum is to be computed.
2596      * @return the BigInteger whose value is the lesser of this BigInteger and
2597      *         {@code val}.  If they are equal, either may be returned.
2598      */
2599     public BigInteger min(BigInteger val) {
2600         return (compareTo(val)<0 ? this : val);
2601     }
2602 
2603     /**
2604      * Returns the maximum of this BigInteger and {@code val}.
2605      *
2606      * @param  val value with which the maximum is to be computed.
2607      * @return the BigInteger whose value is the greater of this and
2608      *         {@code val}.  If they are equal, either may be returned.
2609      */
2610     public BigInteger max(BigInteger val) {
2611         return (compareTo(val)>0 ? this : val);
2612     }
2613 
2614 
2615     // Hash Function
2616 
2617     /**
2618      * Returns the hash code for this BigInteger.
2619      *
2620      * @return hash code for this BigInteger.
2621      */
2622     public int hashCode() {
2623         int hashCode = 0;
2624 
2625         for (int i=0; i<mag.length; i++)
2626             hashCode = (int)(31*hashCode + (mag[i] & LONG_MASK));
2627 
2628         return hashCode * signum;
2629     }
2630 
2631     /**
2632      * Returns the String representation of this BigInteger in the
2633      * given radix.  If the radix is outside the range from {@link
2634      * Character#MIN_RADIX} to {@link Character#MAX_RADIX} inclusive,
2635      * it will default to 10 (as is the case for
2636      * {@code Integer.toString}).  The digit-to-character mapping
2637      * provided by {@code Character.forDigit} is used, and a minus
2638      * sign is prepended if appropriate.  (This representation is
2639      * compatible with the {@link #BigInteger(String, int) (String,
2640      * int)} constructor.)
2641      *
2642      * @param  radix  radix of the String representation.
2643      * @return String representation of this BigInteger in the given radix.
2644      * @see    Integer#toString
2645      * @see    Character#forDigit
2646      * @see    #BigInteger(java.lang.String, int)
2647      */
2648     public String toString(int radix) {
2649         if (signum == 0)
2650             return "0";
2651         if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
2652             radix = 10;
2653 
2654         // Compute upper bound on number of digit groups and allocate space
2655         int maxNumDigitGroups = (4*mag.length + 6)/7;
2656         String digitGroup[] = new String[maxNumDigitGroups];
2657 
2658         // Translate number to string, a digit group at a time
2659         BigInteger tmp = this.abs();
2660         int numGroups = 0;
2661         while (tmp.signum != 0) {
2662             BigInteger d = longRadix[radix];
2663 
2664             MutableBigInteger q = new MutableBigInteger(),
2665                               r = new MutableBigInteger(),
2666                               a = new MutableBigInteger(tmp.mag),
2667                               b = new MutableBigInteger(d.mag);
2668             a.divide(b, q, r);
2669             BigInteger q2 = new BigInteger(q, tmp.signum * d.signum);
2670             BigInteger r2 = new BigInteger(r, tmp.signum * d.signum);
2671 
2672             digitGroup[numGroups++] = Long.toString(r2.longValue(), radix);
2673             tmp = q2;
2674         }
2675 
2676         // Put sign (if any) and first digit group into result buffer
2677         StringBuilder buf = new StringBuilder(numGroups*digitsPerLong[radix]+1);
2678         if (signum<0)
2679             buf.append('-');
2680         buf.append(digitGroup[numGroups-1]);
2681 
2682         // Append remaining digit groups padded with leading zeros
2683         for (int i=numGroups-2; i>=0; i--) {
2684             // Prepend (any) leading zeros for this digit group
2685             int numLeadingZeros = digitsPerLong[radix]-digitGroup[i].length();
2686             if (numLeadingZeros != 0)
2687                 buf.append(zeros[numLeadingZeros]);
2688             buf.append(digitGroup[i]);
2689         }
2690         return buf.toString();
2691     }
2692 
2693     /* zero[i] is a string of i consecutive zeros. */
2694     private static String zeros[] = new String[64];
2695     static {
2696         zeros[63] =
2697             "000000000000000000000000000000000000000000000000000000000000000";
2698         for (int i=0; i<63; i++)
2699             zeros[i] = zeros[63].substring(0, i);
2700     }
2701 
2702     /**
2703      * Returns the decimal String representation of this BigInteger.
2704      * The digit-to-character mapping provided by
2705      * {@code Character.forDigit} is used, and a minus sign is
2706      * prepended if appropriate.  (This representation is compatible
2707      * with the {@link #BigInteger(String) (String)} constructor, and
2708      * allows for String concatenation with Java's + operator.)
2709      *
2710      * @return decimal String representation of this BigInteger.
2711      * @see    Character#forDigit
2712      * @see    #BigInteger(java.lang.String)
2713      */
2714     public String toString() {
2715         return toString(10);
2716     }
2717 
2718     /**
2719      * Returns a byte array containing the two's-complement
2720      * representation of this BigInteger.  The byte array will be in
2721      * <i>big-endian</i> byte-order: the most significant byte is in
2722      * the zeroth element.  The array will contain the minimum number
2723      * of bytes required to represent this BigInteger, including at
2724      * least one sign bit, which is {@code (ceil((this.bitLength() +
2725      * 1)/8))}.  (This representation is compatible with the
2726      * {@link #BigInteger(byte[]) (byte[])} constructor.)
2727      *
2728      * @return a byte array containing the two's-complement representation of
2729      *         this BigInteger.
2730      * @see    #BigInteger(byte[])
2731      */
2732     public byte[] toByteArray() {
2733         int byteLen = bitLength()/8 + 1;
2734         byte[] byteArray = new byte[byteLen];
2735 
2736         for (int i=byteLen-1, bytesCopied=4, nextInt=0, intIndex=0; i>=0; i--) {
2737             if (bytesCopied == 4) {
2738                 nextInt = getInt(intIndex++);
2739                 bytesCopied = 1;
2740             } else {
2741                 nextInt >>>= 8;
2742                 bytesCopied++;
2743             }
2744             byteArray[i] = (byte)nextInt;
2745         }
2746         return byteArray;
2747     }
2748 
2749     /**
2750      * Converts this BigInteger to an {@code int}.  This
2751      * conversion is analogous to a <a
2752      * href="http://java.sun.com/docs/books/jls/second_edition/html/conversions.doc.html#25363"><i>narrowing
2753      * primitive conversion</i></a> from {@code long} to
2754      * {@code int} as defined in the <a
2755      * href="http://java.sun.com/docs/books/jls/html/">Java Language
2756      * Specification</a>: if this BigInteger is too big to fit in an
2757      * {@code int}, only the low-order 32 bits are returned.
2758      * Note that this conversion can lose information about the
2759      * overall magnitude of the BigInteger value as well as return a
2760      * result with the opposite sign.
2761      *
2762      * @return this BigInteger converted to an {@code int}.
2763      */
2764     public int intValue() {
2765         int result = 0;
2766         result = getInt(0);
2767         return result;
2768     }
2769 
2770     /**
2771      * Converts this BigInteger to a {@code long}.  This
2772      * conversion is analogous to a <a
2773      * href="http://java.sun.com/docs/books/jls/second_edition/html/conversions.doc.html#25363"><i>narrowing
2774      * primitive conversion</i></a> from {@code long} to
2775      * {@code int} as defined in the <a
2776      * href="http://java.sun.com/docs/books/jls/html/">Java Language
2777      * Specification</a>: if this BigInteger is too big to fit in a
2778      * {@code long}, only the low-order 64 bits are returned.
2779      * Note that this conversion can lose information about the
2780      * overall magnitude of the BigInteger value as well as return a
2781      * result with the opposite sign.
2782      *
2783      * @return this BigInteger converted to a {@code long}.
2784      */
2785     public long longValue() {
2786         long result = 0;
2787 
2788         for (int i=1; i>=0; i--)
2789             result = (result << 32) + (getInt(i) & LONG_MASK);
2790         return result;
2791     }
2792 
2793     /**
2794      * Converts this BigInteger to a {@code float}.  This
2795      * conversion is similar to the <a
2796      * href="http://java.sun.com/docs/books/jls/second_edition/html/conversions.doc.html#25363"><i>narrowing
2797      * primitive conversion</i></a> from {@code double} to
2798      * {@code float} defined in the <a
2799      * href="http://java.sun.com/docs/books/jls/html/">Java Language
2800      * Specification</a>: if this BigInteger has too great a magnitude
2801      * to represent as a {@code float}, it will be converted to
2802      * {@link Float#NEGATIVE_INFINITY} or {@link
2803      * Float#POSITIVE_INFINITY} as appropriate.  Note that even when
2804      * the return value is finite, this conversion can lose
2805      * information about the precision of the BigInteger value.
2806      *
2807      * @return this BigInteger converted to a {@code float}.
2808      */
2809     public float floatValue() {
2810         // Somewhat inefficient, but guaranteed to work.
2811         return Float.parseFloat(this.toString());
2812     }
2813 
2814     /**
2815      * Converts this BigInteger to a {@code double}.  This
2816      * conversion is similar to the <a
2817      * href="http://java.sun.com/docs/books/jls/second_edition/html/conversions.doc.html#25363"><i>narrowing
2818      * primitive conversion</i></a> from {@code double} to
2819      * {@code float} defined in the <a
2820      * href="http://java.sun.com/docs/books/jls/html/">Java Language
2821      * Specification</a>: if this BigInteger has too great a magnitude
2822      * to represent as a {@code double}, it will be converted to
2823      * {@link Double#NEGATIVE_INFINITY} or {@link
2824      * Double#POSITIVE_INFINITY} as appropriate.  Note that even when
2825      * the return value is finite, this conversion can lose
2826      * information about the precision of the BigInteger value.
2827      *
2828      * @return this BigInteger converted to a {@code double}.
2829      */
2830     public double doubleValue() {
2831         // Somewhat inefficient, but guaranteed to work.
2832         return Double.parseDouble(this.toString());
2833     }
2834 
2835     /**
2836      * Returns a copy of the input array stripped of any leading zero bytes.
2837      */
2838     private static int[] stripLeadingZeroInts(int val[]) {
2839         int byteLength = val.length;
2840         int keep;
2841 
2842         // Find first nonzero byte
2843         for (keep=0; keep<val.length && val[keep]==0; keep++)
2844             ;
2845 
2846         int result[] = new int[val.length - keep];
2847         for(int i=0; i<val.length - keep; i++)
2848             result[i] = val[keep+i];
2849 
2850         return result;
2851     }
2852 
2853     /**
2854      * Returns the input array stripped of any leading zero bytes.
2855      * Since the source is trusted the copying may be skipped.
2856      */
2857     private static int[] trustedStripLeadingZeroInts(int val[]) {
2858         int byteLength = val.length;
2859         int keep;
2860 
2861         // Find first nonzero byte
2862         for (keep=0; keep<val.length && val[keep]==0; keep++)
2863             ;
2864 
2865         // Only perform copy if necessary
2866         if (keep > 0) {
2867             int result[] = new int[val.length - keep];
2868             for(int i=0; i<val.length - keep; i++)
2869                result[i] = val[keep+i];
2870             return result;
2871         }
2872         return val;
2873     }
2874 
2875     /**
2876      * Returns a copy of the input array stripped of any leading zero bytes.
2877      */
2878     private static int[] stripLeadingZeroBytes(byte a[]) {
2879         int byteLength = a.length;
2880         int keep;
2881 
2882         // Find first nonzero byte
2883         for (keep=0; keep<a.length && a[keep]==0; keep++)
2884             ;
2885 
2886         // Allocate new array and copy relevant part of input array
2887         int intLength = ((byteLength - keep) + 3)/4;
2888         int[] result = new int[intLength];
2889         int b = byteLength - 1;
2890         for (int i = intLength-1; i >= 0; i--) {
2891             result[i] = a[b--] & 0xff;
2892             int bytesRemaining = b - keep + 1;
2893             int bytesToTransfer = Math.min(3, bytesRemaining);
2894             for (int j=8; j <= 8*bytesToTransfer; j += 8)
2895                 result[i] |= ((a[b--] & 0xff) << j);
2896         }
2897         return result;
2898     }
2899 
2900     /**
2901      * Takes an array a representing a negative 2's-complement number and
2902      * returns the minimal (no leading zero bytes) unsigned whose value is -a.
2903      */
2904     private static int[] makePositive(byte a[]) {
2905         int keep, k;
2906         int byteLength = a.length;
2907 
2908         // Find first non-sign (0xff) byte of input
2909         for (keep=0; keep<byteLength && a[keep]==-1; keep++)
2910             ;
2911 
2912 
2913         /* Allocate output array.  If all non-sign bytes are 0x00, we must
2914          * allocate space for one extra output byte. */
2915         for (k=keep; k<byteLength && a[k]==0; k++)
2916             ;
2917 
2918         int extraByte = (k==byteLength) ? 1 : 0;
2919         int intLength = ((byteLength - keep + extraByte) + 3)/4;
2920         int result[] = new int[intLength];
2921 
2922         /* Copy one's complement of input into output, leaving extra
2923          * byte (if it exists) == 0x00 */
2924         int b = byteLength - 1;
2925         for (int i = intLength-1; i >= 0; i--) {
2926             result[i] = a[b--] & 0xff;
2927             int numBytesToTransfer = Math.min(3, b-keep+1);
2928             if (numBytesToTransfer < 0)
2929                 numBytesToTransfer = 0;
2930             for (int j=8; j <= 8*numBytesToTransfer; j += 8)
2931                 result[i] |= ((a[b--] & 0xff) << j);
2932 
2933             // Mask indicates which bits must be complemented
2934             int mask = -1 >>> (8*(3-numBytesToTransfer));
2935             result[i] = ~result[i] & mask;
2936         }
2937 
2938         // Add one to one's complement to generate two's complement
2939         for (int i=result.length-1; i>=0; i--) {
2940             result[i] = (int)((result[i] & LONG_MASK) + 1);
2941             if (result[i] != 0)
2942                 break;
2943         }
2944 
2945         return result;
2946     }
2947 
2948     /**
2949      * Takes an array a representing a negative 2's-complement number and
2950      * returns the minimal (no leading zero ints) unsigned whose value is -a.
2951      */
2952     private static int[] makePositive(int a[]) {
2953         int keep, j;
2954 
2955         // Find first non-sign (0xffffffff) int of input
2956         for (keep=0; keep<a.length && a[keep]==-1; keep++)
2957             ;
2958 
2959         /* Allocate output array.  If all non-sign ints are 0x00, we must
2960          * allocate space for one extra output int. */
2961         for (j=keep; j<a.length && a[j]==0; j++)
2962             ;
2963         int extraInt = (j==a.length ? 1 : 0);
2964         int result[] = new int[a.length - keep + extraInt];
2965 
2966         /* Copy one's complement of input into output, leaving extra
2967          * int (if it exists) == 0x00 */
2968         for (int i = keep; i<a.length; i++)
2969             result[i - keep + extraInt] = ~a[i];
2970 
2971         // Add one to one's complement to generate two's complement
2972         for (int i=result.length-1; ++result[i]==0; i--)
2973             ;
2974 
2975         return result;
2976     }
2977 
2978     /*
2979      * The following two arrays are used for fast String conversions.  Both
2980      * are indexed by radix.  The first is the number of digits of the given
2981      * radix that can fit in a Java long without "going negative", i.e., the
2982      * highest integer n such that radix**n < 2**63.  The second is the
2983      * "long radix" that tears each number into "long digits", each of which
2984      * consists of the number of digits in the corresponding element in
2985      * digitsPerLong (longRadix[i] = i**digitPerLong[i]).  Both arrays have
2986      * nonsense values in their 0 and 1 elements, as radixes 0 and 1 are not
2987      * used.
2988      */
2989     private static int digitsPerLong[] = {0, 0,
2990         62, 39, 31, 27, 24, 22, 20, 19, 18, 18, 17, 17, 16, 16, 15, 15, 15, 14,
2991         14, 14, 14, 13, 13, 13, 13, 13, 13, 12, 12, 12, 12, 12, 12, 12, 12};
2992 
2993     private static BigInteger longRadix[] = {null, null,
2994         valueOf(0x4000000000000000L), valueOf(0x383d9170b85ff80bL),
2995         valueOf(0x4000000000000000L), valueOf(0x6765c793fa10079dL),
2996         valueOf(0x41c21cb8e1000000L), valueOf(0x3642798750226111L),
2997         valueOf(0x1000000000000000L), valueOf(0x12bf307ae81ffd59L),
2998         valueOf( 0xde0b6b3a7640000L), valueOf(0x4d28cb56c33fa539L),
2999         valueOf(0x1eca170c00000000L), valueOf(0x780c7372621bd74dL),
3000         valueOf(0x1e39a5057d810000L), valueOf(0x5b27ac993df97701L),
3001         valueOf(0x1000000000000000L), valueOf(0x27b95e997e21d9f1L),
3002         valueOf(0x5da0e1e53c5c8000L), valueOf( 0xb16a458ef403f19L),
3003         valueOf(0x16bcc41e90000000L), valueOf(0x2d04b7fdd9c0ef49L),
3004         valueOf(0x5658597bcaa24000L), valueOf( 0x6feb266931a75b7L),
3005         valueOf( 0xc29e98000000000L), valueOf(0x14adf4b7320334b9L),
3006         valueOf(0x226ed36478bfa000L), valueOf(0x383d9170b85ff80bL),
3007         valueOf(0x5a3c23e39c000000L), valueOf( 0x4e900abb53e6b71L),
3008         valueOf( 0x7600ec618141000L), valueOf( 0xaee5720ee830681L),
3009         valueOf(0x1000000000000000L), valueOf(0x172588ad4f5f0981L),
3010         valueOf(0x211e44f7d02c1000L), valueOf(0x2ee56725f06e5c71L),
3011         valueOf(0x41c21cb8e1000000L)};
3012 
3013     /*
3014      * These two arrays are the integer analogue of above.
3015      */
3016     private static int digitsPerInt[] = {0, 0, 30, 19, 15, 13, 11,
3017         11, 10, 9, 9, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6,
3018         6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5};
3019 
3020     private static int intRadix[] = {0, 0,
3021         0x40000000, 0x4546b3db, 0x40000000, 0x48c27395, 0x159fd800,
3022         0x75db9c97, 0x40000000, 0x17179149, 0x3b9aca00, 0xcc6db61,
3023         0x19a10000, 0x309f1021, 0x57f6c100, 0xa2f1b6f,  0x10000000,
3024         0x18754571, 0x247dbc80, 0x3547667b, 0x4c4b4000, 0x6b5a6e1d,
3025         0x6c20a40,  0x8d2d931,  0xb640000,  0xe8d4a51,  0x1269ae40,
3026         0x17179149, 0x1cb91000, 0x23744899, 0x2b73a840, 0x34e63b41,
3027         0x40000000, 0x4cfa3cc1, 0x5c13d840, 0x6d91b519, 0x39aa400
3028     };
3029 
3030     /**
3031      * These routines provide access to the two's complement representation
3032      * of BigIntegers.
3033      */
3034 
3035     /**
3036      * Returns the length of the two's complement representation in ints,
3037      * including space for at least one sign bit.
3038      */
3039     private int intLength() {
3040         return bitLength()/32 + 1;
3041     }
3042 
3043     /* Returns sign bit */
3044     private int signBit() {
3045         return signum < 0 ? 1 : 0;
3046     }
3047 
3048     /* Returns an int of sign bits */
3049     private int signInt() {
3050         return signum < 0 ? -1 : 0;
3051     }
3052 
3053     /**
3054      * Returns the specified int of the little-endian two's complement
3055      * representation (int 0 is the least significant).  The int number can
3056      * be arbitrarily high (values are logically preceded by infinitely many
3057      * sign ints).
3058      */
3059     private int getInt(int n) {
3060         if (n < 0)
3061             return 0;
3062         if (n >= mag.length)
3063             return signInt();
3064 
3065         int magInt = mag[mag.length-n-1];
3066 
3067         return (signum >= 0 ? magInt :
3068                 (n <= firstNonzeroIntNum() ? -magInt : ~magInt));
3069     }
3070 
3071     /**
3072      * Returns the index of the int that contains the first nonzero int in the
3073      * little-endian binary representation of the magnitude (int 0 is the
3074      * least significant). If the magnitude is zero, return value is undefined.
3075      */
3076      private int firstNonzeroIntNum() {
3077         /*
3078          * Initialize firstNonzeroIntNum field the first time this method is
3079          * executed. This method depends on the atomicity of int modifies;
3080          * without this guarantee, it would have to be synchronized.
3081          */
3082         if (firstNonzeroIntNum == -2) {
3083             // Search for the first nonzero int
3084             int i;
3085             for (i=mag.length-1; i>=0 && mag[i]==0; i--)
3086                 ;
3087             firstNonzeroIntNum = mag.length-i-1;
3088         }
3089         return firstNonzeroIntNum;
3090     }
3091 
3092     /** use serialVersionUID from JDK 1.1. for interoperability */
3093     private static final long serialVersionUID = -8287574255936472291L;
3094 
3095     /**
3096      * Serializable fields for BigInteger.
3097      *
3098      * @serialField signum  int
3099      *              signum of this BigInteger.
3100      * @serialField magnitude int[]
3101      *              magnitude array of this BigInteger.
3102      * @serialField bitCount  int
3103      *              number of bits in this BigInteger
3104      * @serialField bitLength int
3105      *              the number of bits in the minimal two's-complement
3106      *              representation of this BigInteger
3107      * @serialField lowestSetBit int
3108      *              lowest set bit in the twos complement representation
3109      */
3110     private static final ObjectStreamField[] serialPersistentFields = {
3111         new ObjectStreamField("signum", Integer.TYPE),
3112         new ObjectStreamField("magnitude", byte[].class),
3113         new ObjectStreamField("bitCount", Integer.TYPE),
3114         new ObjectStreamField("bitLength", Integer.TYPE),
3115         new ObjectStreamField("firstNonzeroByteNum", Integer.TYPE),
3116         new ObjectStreamField("lowestSetBit", Integer.TYPE)
3117         };
3118 
3119     /**
3120      * Reconstitute the {@code BigInteger} instance from a stream (that is,
3121      * deserialize it). The magnitude is read in as an array of bytes
3122      * for historical reasons, but it is converted to an array of ints
3123      * and the byte array is discarded.
3124      */
3125     private void readObject(java.io.ObjectInputStream s)
3126         throws java.io.IOException, ClassNotFoundException {
3127         /*
3128          * In order to maintain compatibility with previous serialized forms,
3129          * the magnitude of a BigInteger is serialized as an array of bytes.
3130          * The magnitude field is used as a temporary store for the byte array
3131          * that is deserialized. The cached computation fields should be
3132          * transient but are serialized for compatibility reasons.
3133          */
3134 
3135         // prepare to read the alternate persistent fields
3136         ObjectInputStream.GetField fields = s.readFields();
3137 
3138         // Read the alternate persistent fields that we care about
3139         signum = fields.get("signum", -2);
3140         byte[] magnitude = (byte[])fields.get("magnitude", null);
3141 
3142         // Validate signum
3143         if (signum < -1 || signum > 1) {
3144             String message = "BigInteger: Invalid signum value";
3145             if (fields.defaulted("signum"))
3146                 message = "BigInteger: Signum not present in stream";
3147             throw new java.io.StreamCorruptedException(message);
3148         }
3149         if ((magnitude.length==0) != (signum==0)) {
3150             String message = "BigInteger: signum-magnitude mismatch";
3151             if (fields.defaulted("magnitude"))
3152                 message = "BigInteger: Magnitude not present in stream";
3153             throw new java.io.StreamCorruptedException(message);
3154         }
3155 
3156         // Set "cached computation" fields to their initial values
3157         bitCount = bitLength = -1;
3158         lowestSetBit = firstNonzeroByteNum = firstNonzeroIntNum = -2;
3159 
3160         // Calculate mag field from magnitude and discard magnitude
3161         mag = stripLeadingZeroBytes(magnitude);
3162     }
3163 
3164     /**
3165      * Save the {@code BigInteger} instance to a stream.
3166      * The magnitude of a BigInteger is serialized as a byte array for
3167      * historical reasons.
3168      *
3169      * @serialData two necessary fields are written as well as obsolete
3170      *             fields for compatibility with older versions.
3171      */
3172     private void writeObject(ObjectOutputStream s) throws IOException {
3173         // set the values of the Serializable fields
3174         ObjectOutputStream.PutField fields = s.putFields();
3175         fields.put("signum", signum);
3176         fields.put("magnitude", magSerializedForm());
3177         fields.put("bitCount", -1);
3178         fields.put("bitLength", -1);
3179         fields.put("lowestSetBit", -2);
3180         fields.put("firstNonzeroByteNum", -2);
3181 
3182         // save them
3183         s.writeFields();
3184 }
3185 
3186     /**
3187      * Returns the mag array as an array of bytes.
3188      */
3189     private byte[] magSerializedForm() {
3190         int bitLen = (mag.length == 0 ? 0 :
3191                       ((mag.length - 1) << 5) + bitLen(mag[0]));
3192         int byteLen = (bitLen + 7)/8;
3193         byte[] result = new byte[byteLen];
3194 
3195         for (int i=byteLen-1, bytesCopied=4, intIndex=mag.length-1, nextInt=0;
3196              i>=0; i--) {
3197             if (bytesCopied == 4) {
3198                 nextInt = mag[intIndex--];
3199                 bytesCopied = 1;
3200             } else {
3201                 nextInt >>>= 8;
3202                 bytesCopied++;
3203             }
3204             result[i] = (byte)nextInt;
3205         }
3206         return result;
3207     }
3208 }