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