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