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