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