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         /*
1967          * The algorithm used here is adapted from Colin Plumb's C library.
1968          * Technique: Consider the partial products in the multiplication
1969          * of "abcde" by itself:
1970          *
1971          *               a  b  c  d  e
1972          *            *  a  b  c  d  e
1973          *          ==================
1974          *              ae be ce de ee
1975          *           ad bd cd dd de
1976          *        ac bc cc cd ce
1977          *     ab bb bc bd be
1978          *  aa ab ac ad ae
1979          *
1980          * Note that everything above the main diagonal:
1981          *              ae be ce de = (abcd) * e
1982          *           ad bd cd       = (abc) * d
1983          *        ac bc             = (ab) * c
1984          *     ab                   = (a) * b
1985          *
1986          * is a copy of everything below the main diagonal:
1987          *                       de
1988          *                 cd ce
1989          *           bc bd be
1990          *     ab ac ad ae
1991          *
1992          * Thus, the sum is 2 * (off the diagonal) + diagonal.
1993          *
1994          * This is accumulated beginning with the diagonal (which
1995          * consist of the squares of the digits of the input), which is then
1996          * divided by two, the off-diagonal added, and multiplied by two
1997          * again.  The low bit is simply a copy of the low bit of the
1998          * input, so it doesn't need special care.
1999          */
2000         int zlen = len << 1;
2001         if (z == null || z.length < zlen)
2002             z = new int[zlen];
2003 
2004         // Store the squares, right shifted one bit (i.e., divided by 2)
2005         int lastProductLowWord = 0;
2006         for (int j=0, i=0; j < len; j++) {
2007             long piece = (x[j] & LONG_MASK);
2008             long product = piece * piece;
2009             z[i++] = (lastProductLowWord << 31) | (int)(product >>> 33);
2010             z[i++] = (int)(product >>> 1);
2011             lastProductLowWord = (int)product;
2012         }
2013 
2014         // Add in off-diagonal sums
2015         for (int i=len, offset=1; i > 0; i--, offset+=2) {
2016             int t = x[i-1];
2017             t = mulAdd(z, x, offset, i-1, t);
2018             addOne(z, offset-1, i, t);
2019         }
2020 
2021         // Shift back up and set low bit
2022         primitiveLeftShift(z, zlen, 1);
2023         z[zlen-1] |= x[len-1] & 1;
2024 
2025         return z;
2026     }
2027 
2028     /**
2029      * Squares a BigInteger using the Karatsuba squaring algorithm.  It should
2030      * be used when both numbers are larger than a certain threshold (found
2031      * experimentally).  It is a recursive divide-and-conquer algorithm that
2032      * has better asymptotic performance than the algorithm used in
2033      * squareToLen.
2034      */
2035     private BigInteger squareKaratsuba() {
2036         int half = (mag.length+1) / 2;
2037 
2038         BigInteger xl = getLower(half);
2039         BigInteger xh = getUpper(half);
2040 
2041         BigInteger xhs = xh.square();  // xhs = xh^2
2042         BigInteger xls = xl.square();  // xls = xl^2
2043 
2044         // xh^2 << 64  +  (((xl+xh)^2 - (xh^2 + xl^2)) << 32) + xl^2
2045         return xhs.shiftLeft(half*32).add(xl.add(xh).square().subtract(xhs.add(xls))).shiftLeft(half*32).add(xls);
2046     }
2047 
2048     /**
2049      * Squares a BigInteger using the 3-way Toom-Cook squaring algorithm.  It
2050      * should be used when both numbers are larger than a certain threshold
2051      * (found experimentally).  It is a recursive divide-and-conquer algorithm
2052      * that has better asymptotic performance than the algorithm used in
2053      * squareToLen or squareKaratsuba.
2054      */
2055     private BigInteger squareToomCook3() {
2056         int len = mag.length;
2057 
2058         // k is the size (in ints) of the lower-order slices.
2059         int k = (len+2)/3;   // Equal to ceil(largest/3)
2060 
2061         // r is the size (in ints) of the highest-order slice.
2062         int r = len - 2*k;
2063 
2064         // Obtain slices of the numbers. a2 is the most significant
2065         // bits of the number, and a0 the least significant.
2066         BigInteger a0, a1, a2;
2067         a2 = getToomSlice(k, r, 0, len);
2068         a1 = getToomSlice(k, r, 1, len);
2069         a0 = getToomSlice(k, r, 2, len);
2070         BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1;
2071 
2072         v0 = a0.square();
2073         da1 = a2.add(a0);
2074         vm1 = da1.subtract(a1).square();
2075         da1 = da1.add(a1);
2076         v1 = da1.square();
2077         vinf = a2.square();
2078         v2 = da1.add(a2).shiftLeft(1).subtract(a0).square();
2079 
2080         // The algorithm requires two divisions by 2 and one by 3.
2081         // All divisions are known to be exact, that is, they do not produce
2082         // remainders, and all results are positive.  The divisions by 2 are
2083         // implemented as right shifts which are relatively efficient, leaving
2084         // only a division by 3.
2085         // The division by 3 is done by an optimized algorithm for this case.
2086         t2 = v2.subtract(vm1).exactDivideBy3();
2087         tm1 = v1.subtract(vm1).shiftRight(1);
2088         t1 = v1.subtract(v0);
2089         t2 = t2.subtract(t1).shiftRight(1);
2090         t1 = t1.subtract(tm1).subtract(vinf);
2091         t2 = t2.subtract(vinf.shiftLeft(1));
2092         tm1 = tm1.subtract(t2);
2093 
2094         // Number of bits to shift left.
2095         int ss = k*32;
2096 
2097         return vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
2098     }
2099 
2100     // Division
2101 
2102     /**
2103      * Returns a BigInteger whose value is {@code (this / val)}.
2104      *
2105      * @param  val value by which this BigInteger is to be divided.
2106      * @return {@code this / val}
2107      * @throws ArithmeticException if {@code val} is zero.
2108      */
2109     public BigInteger divide(BigInteger val) {
2110         if (val.mag.length < BURNIKEL_ZIEGLER_THRESHOLD ||
2111                 mag.length - val.mag.length < BURNIKEL_ZIEGLER_OFFSET) {
2112             return divideKnuth(val);
2113         } else {
2114             return divideBurnikelZiegler(val);
2115         }
2116     }
2117 
2118     /**
2119      * Returns a BigInteger whose value is {@code (this / val)} using an O(n^2) algorithm from Knuth.
2120      *
2121      * @param  val value by which this BigInteger is to be divided.
2122      * @return {@code this / val}
2123      * @throws ArithmeticException if {@code val} is zero.
2124      * @see MutableBigInteger#divideKnuth(MutableBigInteger, MutableBigInteger, boolean)
2125      */
2126     private BigInteger divideKnuth(BigInteger val) {
2127         MutableBigInteger q = new MutableBigInteger(),
2128                           a = new MutableBigInteger(this.mag),
2129                           b = new MutableBigInteger(val.mag);
2130 
2131         a.divideKnuth(b, q, false);
2132         return q.toBigInteger(this.signum * val.signum);
2133     }
2134 
2135     /**
2136      * Returns an array of two BigIntegers containing {@code (this / val)}
2137      * followed by {@code (this % val)}.
2138      *
2139      * @param  val value by which this BigInteger is to be divided, and the
2140      *         remainder computed.
2141      * @return an array of two BigIntegers: the quotient {@code (this / val)}
2142      *         is the initial element, and the remainder {@code (this % val)}
2143      *         is the final element.
2144      * @throws ArithmeticException if {@code val} is zero.
2145      */
2146     public BigInteger[] divideAndRemainder(BigInteger val) {
2147         if (val.mag.length < BURNIKEL_ZIEGLER_THRESHOLD ||
2148                 mag.length - val.mag.length < BURNIKEL_ZIEGLER_OFFSET) {
2149             return divideAndRemainderKnuth(val);
2150         } else {
2151             return divideAndRemainderBurnikelZiegler(val);
2152         }
2153     }
2154 
2155     /** Long division */
2156     private BigInteger[] divideAndRemainderKnuth(BigInteger val) {
2157         BigInteger[] result = new BigInteger[2];
2158         MutableBigInteger q = new MutableBigInteger(),
2159                           a = new MutableBigInteger(this.mag),
2160                           b = new MutableBigInteger(val.mag);
2161         MutableBigInteger r = a.divideKnuth(b, q);
2162         result[0] = q.toBigInteger(this.signum == val.signum ? 1 : -1);
2163         result[1] = r.toBigInteger(this.signum);
2164         return result;
2165     }
2166 
2167     /**
2168      * Returns a BigInteger whose value is {@code (this % val)}.
2169      *
2170      * @param  val value by which this BigInteger is to be divided, and the
2171      *         remainder computed.
2172      * @return {@code this % val}
2173      * @throws ArithmeticException if {@code val} is zero.
2174      */
2175     public BigInteger remainder(BigInteger val) {
2176         if (val.mag.length < BURNIKEL_ZIEGLER_THRESHOLD ||
2177                 mag.length - val.mag.length < BURNIKEL_ZIEGLER_OFFSET) {
2178             return remainderKnuth(val);
2179         } else {
2180             return remainderBurnikelZiegler(val);
2181         }
2182     }
2183 
2184     /** Long division */
2185     private BigInteger remainderKnuth(BigInteger val) {
2186         MutableBigInteger q = new MutableBigInteger(),
2187                           a = new MutableBigInteger(this.mag),
2188                           b = new MutableBigInteger(val.mag);
2189 
2190         return a.divideKnuth(b, q).toBigInteger(this.signum);
2191     }
2192 
2193     /**
2194      * Calculates {@code this / val} using the Burnikel-Ziegler algorithm.
2195      * @param  val the divisor
2196      * @return {@code this / val}
2197      */
2198     private BigInteger divideBurnikelZiegler(BigInteger val) {
2199         return divideAndRemainderBurnikelZiegler(val)[0];
2200     }
2201 
2202     /**
2203      * Calculates {@code this % val} using the Burnikel-Ziegler algorithm.
2204      * @param val the divisor
2205      * @return {@code this % val}
2206      */
2207     private BigInteger remainderBurnikelZiegler(BigInteger val) {
2208         return divideAndRemainderBurnikelZiegler(val)[1];
2209     }
2210 
2211     /**
2212      * Computes {@code this / val} and {@code this % val} using the
2213      * Burnikel-Ziegler algorithm.
2214      * @param val the divisor
2215      * @return an array containing the quotient and remainder
2216      */
2217     private BigInteger[] divideAndRemainderBurnikelZiegler(BigInteger val) {
2218         MutableBigInteger q = new MutableBigInteger();
2219         MutableBigInteger r = new MutableBigInteger(this).divideAndRemainderBurnikelZiegler(new MutableBigInteger(val), q);
2220         BigInteger qBigInt = q.isZero() ? ZERO : q.toBigInteger(signum*val.signum);
2221         BigInteger rBigInt = r.isZero() ? ZERO : r.toBigInteger(signum);
2222         return new BigInteger[] {qBigInt, rBigInt};
2223     }
2224 
2225     /**
2226      * Returns a BigInteger whose value is <tt>(this<sup>exponent</sup>)</tt>.
2227      * Note that {@code exponent} is an integer rather than a BigInteger.
2228      *
2229      * @param  exponent exponent to which this BigInteger is to be raised.
2230      * @return <tt>this<sup>exponent</sup></tt>
2231      * @throws ArithmeticException {@code exponent} is negative.  (This would
2232      *         cause the operation to yield a non-integer value.)
2233      */
2234     public BigInteger pow(int exponent) {
2235         if (exponent < 0) {
2236             throw new ArithmeticException("Negative exponent");
2237         }
2238         if (signum == 0) {
2239             return (exponent == 0 ? ONE : this);
2240         }
2241 
2242         BigInteger partToSquare = this.abs();
2243 
2244         // Factor out powers of two from the base, as the exponentiation of
2245         // these can be done by left shifts only.
2246         // The remaining part can then be exponentiated faster.  The
2247         // powers of two will be multiplied back at the end.
2248         int powersOfTwo = partToSquare.getLowestSetBit();
2249         long bitsToShift = (long)powersOfTwo * exponent;
2250         if (bitsToShift > Integer.MAX_VALUE) {
2251             reportOverflow();
2252         }
2253 
2254         int remainingBits;
2255 
2256         // Factor the powers of two out quickly by shifting right, if needed.
2257         if (powersOfTwo > 0) {
2258             partToSquare = partToSquare.shiftRight(powersOfTwo);
2259             remainingBits = partToSquare.bitLength();
2260             if (remainingBits == 1) {  // Nothing left but +/- 1?
2261                 if (signum < 0 && (exponent&1) == 1) {
2262                     return NEGATIVE_ONE.shiftLeft(powersOfTwo*exponent);
2263                 } else {
2264                     return ONE.shiftLeft(powersOfTwo*exponent);
2265                 }
2266             }
2267         } else {
2268             remainingBits = partToSquare.bitLength();
2269             if (remainingBits == 1) { // Nothing left but +/- 1?
2270                 if (signum < 0  && (exponent&1) == 1) {
2271                     return NEGATIVE_ONE;
2272                 } else {
2273                     return ONE;
2274                 }
2275             }
2276         }
2277 
2278         // This is a quick way to approximate the size of the result,
2279         // similar to doing log2[n] * exponent.  This will give an upper bound
2280         // of how big the result can be, and which algorithm to use.
2281         long scaleFactor = (long)remainingBits * exponent;
2282 
2283         // Use slightly different algorithms for small and large operands.
2284         // See if the result will safely fit into a long. (Largest 2^63-1)
2285         if (partToSquare.mag.length == 1 && scaleFactor <= 62) {
2286             // Small number algorithm.  Everything fits into a long.
2287             int newSign = (signum <0  && (exponent&1) == 1 ? -1 : 1);
2288             long result = 1;
2289             long baseToPow2 = partToSquare.mag[0] & LONG_MASK;
2290 
2291             int workingExponent = exponent;
2292 
2293             // Perform exponentiation using repeated squaring trick
2294             while (workingExponent != 0) {
2295                 if ((workingExponent & 1) == 1) {
2296                     result = result * baseToPow2;
2297                 }
2298 
2299                 if ((workingExponent >>>= 1) != 0) {
2300                     baseToPow2 = baseToPow2 * baseToPow2;
2301                 }
2302             }
2303 
2304             // Multiply back the powers of two (quickly, by shifting left)
2305             if (powersOfTwo > 0) {
2306                 if (bitsToShift + scaleFactor <= 62) { // Fits in long?
2307                     return valueOf((result << bitsToShift) * newSign);
2308                 } else {
2309                     return valueOf(result*newSign).shiftLeft((int) bitsToShift);
2310                 }
2311             }
2312             else {
2313                 return valueOf(result*newSign);
2314             }
2315         } else {
2316             // Large number algorithm.  This is basically identical to
2317             // the algorithm above, but calls multiply() and square()
2318             // which may use more efficient algorithms for large numbers.
2319             BigInteger answer = ONE;
2320 
2321             int workingExponent = exponent;
2322             // Perform exponentiation using repeated squaring trick
2323             while (workingExponent != 0) {
2324                 if ((workingExponent & 1) == 1) {
2325                     answer = answer.multiply(partToSquare);
2326                 }
2327 
2328                 if ((workingExponent >>>= 1) != 0) {
2329                     partToSquare = partToSquare.square();
2330                 }
2331             }
2332             // Multiply back the (exponentiated) powers of two (quickly,
2333             // by shifting left)
2334             if (powersOfTwo > 0) {
2335                 answer = answer.shiftLeft(powersOfTwo*exponent);
2336             }
2337 
2338             if (signum < 0 && (exponent&1) == 1) {
2339                 return answer.negate();
2340             } else {
2341                 return answer;
2342             }
2343         }
2344     }
2345 
2346     /**
2347      * Returns a BigInteger whose value is the greatest common divisor of
2348      * {@code abs(this)} and {@code abs(val)}.  Returns 0 if
2349      * {@code this == 0 && val == 0}.
2350      *
2351      * @param  val value with which the GCD is to be computed.
2352      * @return {@code GCD(abs(this), abs(val))}
2353      */
2354     public BigInteger gcd(BigInteger val) {
2355         if (val.signum == 0)
2356             return this.abs();
2357         else if (this.signum == 0)
2358             return val.abs();
2359 
2360         MutableBigInteger a = new MutableBigInteger(this);
2361         MutableBigInteger b = new MutableBigInteger(val);
2362 
2363         MutableBigInteger result = a.hybridGCD(b);
2364 
2365         return result.toBigInteger(1);
2366     }
2367 
2368     /**
2369      * Package private method to return bit length for an integer.
2370      */
2371     static int bitLengthForInt(int n) {
2372         return 32 - Integer.numberOfLeadingZeros(n);
2373     }
2374 
2375     /**
2376      * Left shift int array a up to len by n bits. Returns the array that
2377      * results from the shift since space may have to be reallocated.
2378      */
2379     private static int[] leftShift(int[] a, int len, int n) {
2380         int nInts = n >>> 5;
2381         int nBits = n&0x1F;
2382         int bitsInHighWord = bitLengthForInt(a[0]);
2383 
2384         // If shift can be done without recopy, do so
2385         if (n <= (32-bitsInHighWord)) {
2386             primitiveLeftShift(a, len, nBits);
2387             return a;
2388         } else { // Array must be resized
2389             if (nBits <= (32-bitsInHighWord)) {
2390                 int result[] = new int[nInts+len];
2391                 System.arraycopy(a, 0, result, 0, len);
2392                 primitiveLeftShift(result, result.length, nBits);
2393                 return result;
2394             } else {
2395                 int result[] = new int[nInts+len+1];
2396                 System.arraycopy(a, 0, result, 0, len);
2397                 primitiveRightShift(result, result.length, 32 - nBits);
2398                 return result;
2399             }
2400         }
2401     }
2402 
2403     // shifts a up to len right n bits assumes no leading zeros, 0<n<32
2404     static void primitiveRightShift(int[] a, int len, int n) {
2405         int n2 = 32 - n;
2406         for (int i=len-1, c=a[i]; i > 0; i--) {
2407             int b = c;
2408             c = a[i-1];
2409             a[i] = (c << n2) | (b >>> n);
2410         }
2411         a[0] >>>= n;
2412     }
2413 
2414     // shifts a up to len left n bits assumes no leading zeros, 0<=n<32
2415     static void primitiveLeftShift(int[] a, int len, int n) {
2416         if (len == 0 || n == 0)
2417             return;
2418 
2419         int n2 = 32 - n;
2420         for (int i=0, c=a[i], m=i+len-1; i < m; i++) {
2421             int b = c;
2422             c = a[i+1];
2423             a[i] = (b << n) | (c >>> n2);
2424         }
2425         a[len-1] <<= n;
2426     }
2427 
2428     /**
2429      * Calculate bitlength of contents of the first len elements an int array,
2430      * assuming there are no leading zero ints.
2431      */
2432     private static int bitLength(int[] val, int len) {
2433         if (len == 0)
2434             return 0;
2435         return ((len - 1) << 5) + bitLengthForInt(val[0]);
2436     }
2437 
2438     /**
2439      * Returns a BigInteger whose value is the absolute value of this
2440      * BigInteger.
2441      *
2442      * @return {@code abs(this)}
2443      */
2444     public BigInteger abs() {
2445         return (signum >= 0 ? this : this.negate());
2446     }
2447 
2448     /**
2449      * Returns a BigInteger whose value is {@code (-this)}.
2450      *
2451      * @return {@code -this}
2452      */
2453     public BigInteger negate() {
2454         return new BigInteger(this.mag, -this.signum);
2455     }
2456 
2457     /**
2458      * Returns the signum function of this BigInteger.
2459      *
2460      * @return -1, 0 or 1 as the value of this BigInteger is negative, zero or
2461      *         positive.
2462      */
2463     public int signum() {
2464         return this.signum;
2465     }
2466 
2467     // Modular Arithmetic Operations
2468 
2469     /**
2470      * Returns a BigInteger whose value is {@code (this mod m}).  This method
2471      * differs from {@code remainder} in that it always returns a
2472      * <i>non-negative</i> BigInteger.
2473      *
2474      * @param  m the modulus.
2475      * @return {@code this mod m}
2476      * @throws ArithmeticException {@code m} &le; 0
2477      * @see    #remainder
2478      */
2479     public BigInteger mod(BigInteger m) {
2480         if (m.signum <= 0)
2481             throw new ArithmeticException("BigInteger: modulus not positive");
2482 
2483         BigInteger result = this.remainder(m);
2484         return (result.signum >= 0 ? result : result.add(m));
2485     }
2486 
2487     /**
2488      * Returns a BigInteger whose value is
2489      * <tt>(this<sup>exponent</sup> mod m)</tt>.  (Unlike {@code pow}, this
2490      * method permits negative exponents.)
2491      *
2492      * @param  exponent the exponent.
2493      * @param  m the modulus.
2494      * @return <tt>this<sup>exponent</sup> mod m</tt>
2495      * @throws ArithmeticException {@code m} &le; 0 or the exponent is
2496      *         negative and this BigInteger is not <i>relatively
2497      *         prime</i> to {@code m}.
2498      * @see    #modInverse
2499      */
2500     public BigInteger modPow(BigInteger exponent, BigInteger m) {
2501         if (m.signum <= 0)
2502             throw new ArithmeticException("BigInteger: modulus not positive");
2503 
2504         // Trivial cases
2505         if (exponent.signum == 0)
2506             return (m.equals(ONE) ? ZERO : ONE);
2507 
2508         if (this.equals(ONE))
2509             return (m.equals(ONE) ? ZERO : ONE);
2510 
2511         if (this.equals(ZERO) && exponent.signum >= 0)
2512             return ZERO;
2513 
2514         if (this.equals(negConst[1]) && (!exponent.testBit(0)))
2515             return (m.equals(ONE) ? ZERO : ONE);
2516 
2517         boolean invertResult;
2518         if ((invertResult = (exponent.signum < 0)))
2519             exponent = exponent.negate();
2520 
2521         BigInteger base = (this.signum < 0 || this.compareTo(m) >= 0
2522                            ? this.mod(m) : this);
2523         BigInteger result;
2524         if (m.testBit(0)) { // odd modulus
2525             result = base.oddModPow(exponent, m);
2526         } else {
2527             /*
2528              * Even modulus.  Tear it into an "odd part" (m1) and power of two
2529              * (m2), exponentiate mod m1, manually exponentiate mod m2, and
2530              * use Chinese Remainder Theorem to combine results.
2531              */
2532 
2533             // Tear m apart into odd part (m1) and power of 2 (m2)
2534             int p = m.getLowestSetBit();   // Max pow of 2 that divides m
2535 
2536             BigInteger m1 = m.shiftRight(p);  // m/2**p
2537             BigInteger m2 = ONE.shiftLeft(p); // 2**p
2538 
2539             // Calculate new base from m1
2540             BigInteger base2 = (this.signum < 0 || this.compareTo(m1) >= 0
2541                                 ? this.mod(m1) : this);
2542 
2543             // Caculate (base ** exponent) mod m1.
2544             BigInteger a1 = (m1.equals(ONE) ? ZERO :
2545                              base2.oddModPow(exponent, m1));
2546 
2547             // Calculate (this ** exponent) mod m2
2548             BigInteger a2 = base.modPow2(exponent, p);
2549 
2550             // Combine results using Chinese Remainder Theorem
2551             BigInteger y1 = m2.modInverse(m1);
2552             BigInteger y2 = m1.modInverse(m2);
2553 
2554             if (m.mag.length < MAX_MAG_LENGTH / 2) {
2555                 result = a1.multiply(m2).multiply(y1).add(a2.multiply(m1).multiply(y2)).mod(m);
2556             } else {
2557                 MutableBigInteger t1 = new MutableBigInteger();
2558                 new MutableBigInteger(a1.multiply(m2)).multiply(new MutableBigInteger(y1), t1);
2559                 MutableBigInteger t2 = new MutableBigInteger();
2560                 new MutableBigInteger(a2.multiply(m1)).multiply(new MutableBigInteger(y2), t2);
2561                 t1.add(t2);
2562                 MutableBigInteger q = new MutableBigInteger();
2563                 result = t1.divide(new MutableBigInteger(m), q).toBigInteger();
2564             }
2565         }
2566 
2567         return (invertResult ? result.modInverse(m) : result);
2568     }
2569 
2570     static int[] bnExpModThreshTable = {7, 25, 81, 241, 673, 1793,
2571                                                 Integer.MAX_VALUE}; // Sentinel
2572 
2573     /**
2574      * Returns a BigInteger whose value is x to the power of y mod z.
2575      * Assumes: z is odd && x < z.
2576      */
2577     private BigInteger oddModPow(BigInteger y, BigInteger z) {
2578     /*
2579      * The algorithm is adapted from Colin Plumb's C library.
2580      *
2581      * The window algorithm:
2582      * The idea is to keep a running product of b1 = n^(high-order bits of exp)
2583      * and then keep appending exponent bits to it.  The following patterns
2584      * apply to a 3-bit window (k = 3):
2585      * To append   0: square
2586      * To append   1: square, multiply by n^1
2587      * To append  10: square, multiply by n^1, square
2588      * To append  11: square, square, multiply by n^3
2589      * To append 100: square, multiply by n^1, square, square
2590      * To append 101: square, square, square, multiply by n^5
2591      * To append 110: square, square, multiply by n^3, square
2592      * To append 111: square, square, square, multiply by n^7
2593      *
2594      * Since each pattern involves only one multiply, the longer the pattern
2595      * the better, except that a 0 (no multiplies) can be appended directly.
2596      * We precompute a table of odd powers of n, up to 2^k, and can then
2597      * multiply k bits of exponent at a time.  Actually, assuming random
2598      * exponents, there is on average one zero bit between needs to
2599      * multiply (1/2 of the time there's none, 1/4 of the time there's 1,
2600      * 1/8 of the time, there's 2, 1/32 of the time, there's 3, etc.), so
2601      * you have to do one multiply per k+1 bits of exponent.
2602      *
2603      * The loop walks down the exponent, squaring the result buffer as
2604      * it goes.  There is a wbits+1 bit lookahead buffer, buf, that is
2605      * filled with the upcoming exponent bits.  (What is read after the
2606      * end of the exponent is unimportant, but it is filled with zero here.)
2607      * When the most-significant bit of this buffer becomes set, i.e.
2608      * (buf & tblmask) != 0, we have to decide what pattern to multiply
2609      * by, and when to do it.  We decide, remember to do it in future
2610      * after a suitable number of squarings have passed (e.g. a pattern
2611      * of "100" in the buffer requires that we multiply by n^1 immediately;
2612      * a pattern of "110" calls for multiplying by n^3 after one more
2613      * squaring), clear the buffer, and continue.
2614      *
2615      * When we start, there is one more optimization: the result buffer
2616      * is implcitly one, so squaring it or multiplying by it can be
2617      * optimized away.  Further, if we start with a pattern like "100"
2618      * in the lookahead window, rather than placing n into the buffer
2619      * and then starting to square it, we have already computed n^2
2620      * to compute the odd-powers table, so we can place that into
2621      * the buffer and save a squaring.
2622      *
2623      * This means that if you have a k-bit window, to compute n^z,
2624      * where z is the high k bits of the exponent, 1/2 of the time
2625      * it requires no squarings.  1/4 of the time, it requires 1
2626      * squaring, ... 1/2^(k-1) of the time, it reqires k-2 squarings.
2627      * And the remaining 1/2^(k-1) of the time, the top k bits are a
2628      * 1 followed by k-1 0 bits, so it again only requires k-2
2629      * squarings, not k-1.  The average of these is 1.  Add that
2630      * to the one squaring we have to do to compute the table,
2631      * and you'll see that a k-bit window saves k-2 squarings
2632      * as well as reducing the multiplies.  (It actually doesn't
2633      * hurt in the case k = 1, either.)
2634      */
2635         // Special case for exponent of one
2636         if (y.equals(ONE))
2637             return this;
2638 
2639         // Special case for base of zero
2640         if (signum == 0)
2641             return ZERO;
2642 
2643         int[] base = mag.clone();
2644         int[] exp = y.mag;
2645         int[] mod = z.mag;
2646         int modLen = mod.length;
2647 
2648         // Select an appropriate window size
2649         int wbits = 0;
2650         int ebits = bitLength(exp, exp.length);
2651         // if exponent is 65537 (0x10001), use minimum window size
2652         if ((ebits != 17) || (exp[0] != 65537)) {
2653             while (ebits > bnExpModThreshTable[wbits]) {
2654                 wbits++;
2655             }
2656         }
2657 
2658         // Calculate appropriate table size
2659         int tblmask = 1 << wbits;
2660 
2661         // Allocate table for precomputed odd powers of base in Montgomery form
2662         int[][] table = new int[tblmask][];
2663         for (int i=0; i < tblmask; i++)
2664             table[i] = new int[modLen];
2665 
2666         // Compute the modular inverse
2667         int inv = -MutableBigInteger.inverseMod32(mod[modLen-1]);
2668 
2669         // Convert base to Montgomery form
2670         int[] a = leftShift(base, base.length, modLen << 5);
2671 
2672         MutableBigInteger q = new MutableBigInteger(),
2673                           a2 = new MutableBigInteger(a),
2674                           b2 = new MutableBigInteger(mod);
2675 
2676         MutableBigInteger r= a2.divide(b2, q);
2677         table[0] = r.toIntArray();
2678 
2679         // Pad table[0] with leading zeros so its length is at least modLen
2680         if (table[0].length < modLen) {
2681            int offset = modLen - table[0].length;
2682            int[] t2 = new int[modLen];
2683            for (int i=0; i < table[0].length; i++)
2684                t2[i+offset] = table[0][i];
2685            table[0] = t2;
2686         }
2687 
2688         // Set b to the square of the base
2689         int[] b = squareToLen(table[0], modLen, null);
2690         b = montReduce(b, mod, modLen, inv);
2691 
2692         // Set t to high half of b
2693         int[] t = Arrays.copyOf(b, modLen);
2694 
2695         // Fill in the table with odd powers of the base
2696         for (int i=1; i < tblmask; i++) {
2697             int[] prod = multiplyToLen(t, modLen, table[i-1], modLen, null);
2698             table[i] = montReduce(prod, mod, modLen, inv);
2699         }
2700 
2701         // Pre load the window that slides over the exponent
2702         int bitpos = 1 << ((ebits-1) & (32-1));
2703 
2704         int buf = 0;
2705         int elen = exp.length;
2706         int eIndex = 0;
2707         for (int i = 0; i <= wbits; i++) {
2708             buf = (buf << 1) | (((exp[eIndex] & bitpos) != 0)?1:0);
2709             bitpos >>>= 1;
2710             if (bitpos == 0) {
2711                 eIndex++;
2712                 bitpos = 1 << (32-1);
2713                 elen--;
2714             }
2715         }
2716 
2717         int multpos = ebits;
2718 
2719         // The first iteration, which is hoisted out of the main loop
2720         ebits--;
2721         boolean isone = true;
2722 
2723         multpos = ebits - wbits;
2724         while ((buf & 1) == 0) {
2725             buf >>>= 1;
2726             multpos++;
2727         }
2728 
2729         int[] mult = table[buf >>> 1];
2730 
2731         buf = 0;
2732         if (multpos == ebits)
2733             isone = false;
2734 
2735         // The main loop
2736         while (true) {
2737             ebits--;
2738             // Advance the window
2739             buf <<= 1;
2740 
2741             if (elen != 0) {
2742                 buf |= ((exp[eIndex] & bitpos) != 0) ? 1 : 0;
2743                 bitpos >>>= 1;
2744                 if (bitpos == 0) {
2745                     eIndex++;
2746                     bitpos = 1 << (32-1);
2747                     elen--;
2748                 }
2749             }
2750 
2751             // Examine the window for pending multiplies
2752             if ((buf & tblmask) != 0) {
2753                 multpos = ebits - wbits;
2754                 while ((buf & 1) == 0) {
2755                     buf >>>= 1;
2756                     multpos++;
2757                 }
2758                 mult = table[buf >>> 1];
2759                 buf = 0;
2760             }
2761 
2762             // Perform multiply
2763             if (ebits == multpos) {
2764                 if (isone) {
2765                     b = mult.clone();
2766                     isone = false;
2767                 } else {
2768                     t = b;
2769                     a = multiplyToLen(t, modLen, mult, modLen, a);
2770                     a = montReduce(a, mod, modLen, inv);
2771                     t = a; a = b; b = t;
2772                 }
2773             }
2774 
2775             // Check if done
2776             if (ebits == 0)
2777                 break;
2778 
2779             // Square the input
2780             if (!isone) {
2781                 t = b;
2782                 a = squareToLen(t, modLen, a);
2783                 a = montReduce(a, mod, modLen, inv);
2784                 t = a; a = b; b = t;
2785             }
2786         }
2787 
2788         // Convert result out of Montgomery form and return
2789         int[] t2 = new int[2*modLen];
2790         System.arraycopy(b, 0, t2, modLen, modLen);
2791 
2792         b = montReduce(t2, mod, modLen, inv);
2793 
2794         t2 = Arrays.copyOf(b, modLen);
2795 
2796         return new BigInteger(1, t2);
2797     }
2798 
2799     /**
2800      * Montgomery reduce n, modulo mod.  This reduces modulo mod and divides
2801      * by 2^(32*mlen). Adapted from Colin Plumb's C library.
2802      */
2803     private static int[] montReduce(int[] n, int[] mod, int mlen, int inv) {
2804         int c=0;
2805         int len = mlen;
2806         int offset=0;
2807 
2808         do {
2809             int nEnd = n[n.length-1-offset];
2810             int carry = mulAdd(n, mod, offset, mlen, inv * nEnd);
2811             c += addOne(n, offset, mlen, carry);
2812             offset++;
2813         } while (--len > 0);
2814 
2815         while (c > 0)
2816             c += subN(n, mod, mlen);
2817 
2818         while (intArrayCmpToLen(n, mod, mlen) >= 0)
2819             subN(n, mod, mlen);
2820 
2821         return n;
2822     }
2823 
2824 
2825     /*
2826      * Returns -1, 0 or +1 as big-endian unsigned int array arg1 is less than,
2827      * equal to, or greater than arg2 up to length len.
2828      */
2829     private static int intArrayCmpToLen(int[] arg1, int[] arg2, int len) {
2830         for (int i=0; i < len; i++) {
2831             long b1 = arg1[i] & LONG_MASK;
2832             long b2 = arg2[i] & LONG_MASK;
2833             if (b1 < b2)
2834                 return -1;
2835             if (b1 > b2)
2836                 return 1;
2837         }
2838         return 0;
2839     }
2840 
2841     /**
2842      * Subtracts two numbers of same length, returning borrow.
2843      */
2844     private static int subN(int[] a, int[] b, int len) {
2845         long sum = 0;
2846 
2847         while (--len >= 0) {
2848             sum = (a[len] & LONG_MASK) -
2849                  (b[len] & LONG_MASK) + (sum >> 32);
2850             a[len] = (int)sum;
2851         }
2852 
2853         return (int)(sum >> 32);
2854     }
2855 
2856     /**
2857      * Multiply an array by one word k and add to result, return the carry
2858      */
2859     static int mulAdd(int[] out, int[] in, int offset, int len, int k) {
2860         long kLong = k & LONG_MASK;
2861         long carry = 0;
2862 
2863         offset = out.length-offset - 1;
2864         for (int j=len-1; j >= 0; j--) {
2865             long product = (in[j] & LONG_MASK) * kLong +
2866                            (out[offset] & LONG_MASK) + carry;
2867             out[offset--] = (int)product;
2868             carry = product >>> 32;
2869         }
2870         return (int)carry;
2871     }
2872 
2873     /**
2874      * Add one word to the number a mlen words into a. Return the resulting
2875      * carry.
2876      */
2877     static int addOne(int[] a, int offset, int mlen, int carry) {
2878         offset = a.length-1-mlen-offset;
2879         long t = (a[offset] & LONG_MASK) + (carry & LONG_MASK);
2880 
2881         a[offset] = (int)t;
2882         if ((t >>> 32) == 0)
2883             return 0;
2884         while (--mlen >= 0) {
2885             if (--offset < 0) { // Carry out of number
2886                 return 1;
2887             } else {
2888                 a[offset]++;
2889                 if (a[offset] != 0)
2890                     return 0;
2891             }
2892         }
2893         return 1;
2894     }
2895 
2896     /**
2897      * Returns a BigInteger whose value is (this ** exponent) mod (2**p)
2898      */
2899     private BigInteger modPow2(BigInteger exponent, int p) {
2900         /*
2901          * Perform exponentiation using repeated squaring trick, chopping off
2902          * high order bits as indicated by modulus.
2903          */
2904         BigInteger result = ONE;
2905         BigInteger baseToPow2 = this.mod2(p);
2906         int expOffset = 0;
2907 
2908         int limit = exponent.bitLength();
2909 
2910         if (this.testBit(0))
2911            limit = (p-1) < limit ? (p-1) : limit;
2912 
2913         while (expOffset < limit) {
2914             if (exponent.testBit(expOffset))
2915                 result = result.multiply(baseToPow2).mod2(p);
2916             expOffset++;
2917             if (expOffset < limit)
2918                 baseToPow2 = baseToPow2.square().mod2(p);
2919         }
2920 
2921         return result;
2922     }
2923 
2924     /**
2925      * Returns a BigInteger whose value is this mod(2**p).
2926      * Assumes that this {@code BigInteger >= 0} and {@code p > 0}.
2927      */
2928     private BigInteger mod2(int p) {
2929         if (bitLength() <= p)
2930             return this;
2931 
2932         // Copy remaining ints of mag
2933         int numInts = (p + 31) >>> 5;
2934         int[] mag = new int[numInts];
2935         System.arraycopy(this.mag, (this.mag.length - numInts), mag, 0, numInts);
2936 
2937         // Mask out any excess bits
2938         int excessBits = (numInts << 5) - p;
2939         mag[0] &= (1L << (32-excessBits)) - 1;
2940 
2941         return (mag[0] == 0 ? new BigInteger(1, mag) : new BigInteger(mag, 1));
2942     }
2943 
2944     /**
2945      * Returns a BigInteger whose value is {@code (this}<sup>-1</sup> {@code mod m)}.
2946      *
2947      * @param  m the modulus.
2948      * @return {@code this}<sup>-1</sup> {@code mod m}.
2949      * @throws ArithmeticException {@code  m} &le; 0, or this BigInteger
2950      *         has no multiplicative inverse mod m (that is, this BigInteger
2951      *         is not <i>relatively prime</i> to m).
2952      */
2953     public BigInteger modInverse(BigInteger m) {
2954         if (m.signum != 1)
2955             throw new ArithmeticException("BigInteger: modulus not positive");
2956 
2957         if (m.equals(ONE))
2958             return ZERO;
2959 
2960         // Calculate (this mod m)
2961         BigInteger modVal = this;
2962         if (signum < 0 || (this.compareMagnitude(m) >= 0))
2963             modVal = this.mod(m);
2964 
2965         if (modVal.equals(ONE))
2966             return ONE;
2967 
2968         MutableBigInteger a = new MutableBigInteger(modVal);
2969         MutableBigInteger b = new MutableBigInteger(m);
2970 
2971         MutableBigInteger result = a.mutableModInverse(b);
2972         return result.toBigInteger(1);
2973     }
2974 
2975     // Shift Operations
2976 
2977     /**
2978      * Returns a BigInteger whose value is {@code (this << n)}.
2979      * The shift distance, {@code n}, may be negative, in which case
2980      * this method performs a right shift.
2981      * (Computes <tt>floor(this * 2<sup>n</sup>)</tt>.)
2982      *
2983      * @param  n shift distance, in bits.
2984      * @return {@code this << n}
2985      * @see #shiftRight
2986      */
2987     public BigInteger shiftLeft(int n) {
2988         if (signum == 0)
2989             return ZERO;
2990         if (n > 0) {
2991             return new BigInteger(shiftLeft(mag, n), signum);
2992         } else if (n == 0) {
2993             return this;
2994         } else {
2995             // Possible int overflow in (-n) is not a trouble,
2996             // because shiftRightImpl considers its argument unsigned
2997             return shiftRightImpl(-n);
2998         }
2999     }
3000 
3001     /**
3002      * Returns a magnitude array whose value is {@code (mag << n)}.
3003      * The shift distance, {@code n}, is considered unnsigned.
3004      * (Computes <tt>this * 2<sup>n</sup></tt>.)
3005      *
3006      * @param mag magnitude, the most-significant int ({@code mag[0]}) must be non-zero.
3007      * @param  n unsigned shift distance, in bits.
3008      * @return {@code mag << n}
3009      */
3010     private static int[] shiftLeft(int[] mag, int n) {
3011         int nInts = n >>> 5;
3012         int nBits = n & 0x1f;
3013         int magLen = mag.length;
3014         int newMag[] = null;
3015 
3016         if (nBits == 0) {
3017             newMag = new int[magLen + nInts];
3018             System.arraycopy(mag, 0, newMag, 0, magLen);
3019         } else {
3020             int i = 0;
3021             int nBits2 = 32 - nBits;
3022             int highBits = mag[0] >>> nBits2;
3023             if (highBits != 0) {
3024                 newMag = new int[magLen + nInts + 1];
3025                 newMag[i++] = highBits;
3026             } else {
3027                 newMag = new int[magLen + nInts];
3028             }
3029             int j=0;
3030             while (j < magLen-1)
3031                 newMag[i++] = mag[j++] << nBits | mag[j] >>> nBits2;
3032             newMag[i] = mag[j] << nBits;
3033         }
3034         return newMag;
3035     }
3036 
3037     /**
3038      * Returns a BigInteger whose value is {@code (this >> n)}.  Sign
3039      * extension is performed.  The shift distance, {@code n}, may be
3040      * negative, in which case this method performs a left shift.
3041      * (Computes <tt>floor(this / 2<sup>n</sup>)</tt>.)
3042      *
3043      * @param  n shift distance, in bits.
3044      * @return {@code this >> n}
3045      * @see #shiftLeft
3046      */
3047     public BigInteger shiftRight(int n) {
3048         if (signum == 0)
3049             return ZERO;
3050         if (n > 0) {
3051             return shiftRightImpl(n);
3052         } else if (n == 0) {
3053             return this;
3054         } else {
3055             // Possible int overflow in {@code -n} is not a trouble,
3056             // because shiftLeft considers its argument unsigned
3057             return new BigInteger(shiftLeft(mag, -n), signum);
3058         }
3059     }
3060 
3061     /**
3062      * Returns a BigInteger whose value is {@code (this >> n)}. The shift
3063      * distance, {@code n}, is considered unsigned.
3064      * (Computes <tt>floor(this * 2<sup>-n</sup>)</tt>.)
3065      *
3066      * @param  n unsigned shift distance, in bits.
3067      * @return {@code this >> n}
3068      */
3069     private BigInteger shiftRightImpl(int n) {
3070         int nInts = n >>> 5;
3071         int nBits = n & 0x1f;
3072         int magLen = mag.length;
3073         int newMag[] = null;
3074 
3075         // Special case: entire contents shifted off the end
3076         if (nInts >= magLen)
3077             return (signum >= 0 ? ZERO : negConst[1]);
3078 
3079         if (nBits == 0) {
3080             int newMagLen = magLen - nInts;
3081             newMag = Arrays.copyOf(mag, newMagLen);
3082         } else {
3083             int i = 0;
3084             int highBits = mag[0] >>> nBits;
3085             if (highBits != 0) {
3086                 newMag = new int[magLen - nInts];
3087                 newMag[i++] = highBits;
3088             } else {
3089                 newMag = new int[magLen - nInts -1];
3090             }
3091 
3092             int nBits2 = 32 - nBits;
3093             int j=0;
3094             while (j < magLen - nInts - 1)
3095                 newMag[i++] = (mag[j++] << nBits2) | (mag[j] >>> nBits);
3096         }
3097 
3098         if (signum < 0) {
3099             // Find out whether any one-bits were shifted off the end.
3100             boolean onesLost = false;
3101             for (int i=magLen-1, j=magLen-nInts; i >= j && !onesLost; i--)
3102                 onesLost = (mag[i] != 0);
3103             if (!onesLost && nBits != 0)
3104                 onesLost = (mag[magLen - nInts - 1] << (32 - nBits) != 0);
3105 
3106             if (onesLost)
3107                 newMag = javaIncrement(newMag);
3108         }
3109 
3110         return new BigInteger(newMag, signum);
3111     }
3112 
3113     int[] javaIncrement(int[] val) {
3114         int lastSum = 0;
3115         for (int i=val.length-1;  i >= 0 && lastSum == 0; i--)
3116             lastSum = (val[i] += 1);
3117         if (lastSum == 0) {
3118             val = new int[val.length+1];
3119             val[0] = 1;
3120         }
3121         return val;
3122     }
3123 
3124     // Bitwise Operations
3125 
3126     /**
3127      * Returns a BigInteger whose value is {@code (this & val)}.  (This
3128      * method returns a negative BigInteger if and only if this and val are
3129      * both negative.)
3130      *
3131      * @param val value to be AND'ed with this BigInteger.
3132      * @return {@code this & val}
3133      */
3134     public BigInteger and(BigInteger val) {
3135         int[] result = new int[Math.max(intLength(), val.intLength())];
3136         for (int i=0; i < result.length; i++)
3137             result[i] = (getInt(result.length-i-1)
3138                          & val.getInt(result.length-i-1));
3139 
3140         return valueOf(result);
3141     }
3142 
3143     /**
3144      * Returns a BigInteger whose value is {@code (this | val)}.  (This method
3145      * returns a negative BigInteger if and only if either this or val is
3146      * negative.)
3147      *
3148      * @param val value to be OR'ed with this BigInteger.
3149      * @return {@code this | val}
3150      */
3151     public BigInteger or(BigInteger val) {
3152         int[] result = new int[Math.max(intLength(), val.intLength())];
3153         for (int i=0; i < result.length; i++)
3154             result[i] = (getInt(result.length-i-1)
3155                          | val.getInt(result.length-i-1));
3156 
3157         return valueOf(result);
3158     }
3159 
3160     /**
3161      * Returns a BigInteger whose value is {@code (this ^ val)}.  (This method
3162      * returns a negative BigInteger if and only if exactly one of this and
3163      * val are negative.)
3164      *
3165      * @param val value to be XOR'ed with this BigInteger.
3166      * @return {@code this ^ val}
3167      */
3168     public BigInteger xor(BigInteger val) {
3169         int[] result = new int[Math.max(intLength(), val.intLength())];
3170         for (int i=0; i < result.length; i++)
3171             result[i] = (getInt(result.length-i-1)
3172                          ^ val.getInt(result.length-i-1));
3173 
3174         return valueOf(result);
3175     }
3176 
3177     /**
3178      * Returns a BigInteger whose value is {@code (~this)}.  (This method
3179      * returns a negative value if and only if this BigInteger is
3180      * non-negative.)
3181      *
3182      * @return {@code ~this}
3183      */
3184     public BigInteger not() {
3185         int[] result = new int[intLength()];
3186         for (int i=0; i < result.length; i++)
3187             result[i] = ~getInt(result.length-i-1);
3188 
3189         return valueOf(result);
3190     }
3191 
3192     /**
3193      * Returns a BigInteger whose value is {@code (this & ~val)}.  This
3194      * method, which is equivalent to {@code and(val.not())}, is provided as
3195      * a convenience for masking operations.  (This method returns a negative
3196      * BigInteger if and only if {@code this} is negative and {@code val} is
3197      * positive.)
3198      *
3199      * @param val value to be complemented and AND'ed with this BigInteger.
3200      * @return {@code this & ~val}
3201      */
3202     public BigInteger andNot(BigInteger val) {
3203         int[] result = new int[Math.max(intLength(), val.intLength())];
3204         for (int i=0; i < result.length; i++)
3205             result[i] = (getInt(result.length-i-1)
3206                          & ~val.getInt(result.length-i-1));
3207 
3208         return valueOf(result);
3209     }
3210 
3211 
3212     // Single Bit Operations
3213 
3214     /**
3215      * Returns {@code true} if and only if the designated bit is set.
3216      * (Computes {@code ((this & (1<<n)) != 0)}.)
3217      *
3218      * @param  n index of bit to test.
3219      * @return {@code true} if and only if the designated bit is set.
3220      * @throws ArithmeticException {@code n} is negative.
3221      */
3222     public boolean testBit(int n) {
3223         if (n < 0)
3224             throw new ArithmeticException("Negative bit address");
3225 
3226         return (getInt(n >>> 5) & (1 << (n & 31))) != 0;
3227     }
3228 
3229     /**
3230      * Returns a BigInteger whose value is equivalent to this BigInteger
3231      * with the designated bit set.  (Computes {@code (this | (1<<n))}.)
3232      *
3233      * @param  n index of bit to set.
3234      * @return {@code this | (1<<n)}
3235      * @throws ArithmeticException {@code n} is negative.
3236      */
3237     public BigInteger setBit(int n) {
3238         if (n < 0)
3239             throw new ArithmeticException("Negative bit address");
3240 
3241         int intNum = n >>> 5;
3242         int[] result = new int[Math.max(intLength(), intNum+2)];
3243 
3244         for (int i=0; i < result.length; i++)
3245             result[result.length-i-1] = getInt(i);
3246 
3247         result[result.length-intNum-1] |= (1 << (n & 31));
3248 
3249         return valueOf(result);
3250     }
3251 
3252     /**
3253      * Returns a BigInteger whose value is equivalent to this BigInteger
3254      * with the designated bit cleared.
3255      * (Computes {@code (this & ~(1<<n))}.)
3256      *
3257      * @param  n index of bit to clear.
3258      * @return {@code this & ~(1<<n)}
3259      * @throws ArithmeticException {@code n} is negative.
3260      */
3261     public BigInteger clearBit(int n) {
3262         if (n < 0)
3263             throw new ArithmeticException("Negative bit address");
3264 
3265         int intNum = n >>> 5;
3266         int[] result = new int[Math.max(intLength(), ((n + 1) >>> 5) + 1)];
3267 
3268         for (int i=0; i < result.length; i++)
3269             result[result.length-i-1] = getInt(i);
3270 
3271         result[result.length-intNum-1] &= ~(1 << (n & 31));
3272 
3273         return valueOf(result);
3274     }
3275 
3276     /**
3277      * Returns a BigInteger whose value is equivalent to this BigInteger
3278      * with the designated bit flipped.
3279      * (Computes {@code (this ^ (1<<n))}.)
3280      *
3281      * @param  n index of bit to flip.
3282      * @return {@code this ^ (1<<n)}
3283      * @throws ArithmeticException {@code n} is negative.
3284      */
3285     public BigInteger flipBit(int n) {
3286         if (n < 0)
3287             throw new ArithmeticException("Negative bit address");
3288 
3289         int intNum = n >>> 5;
3290         int[] result = new int[Math.max(intLength(), intNum+2)];
3291 
3292         for (int i=0; i < result.length; i++)
3293             result[result.length-i-1] = getInt(i);
3294 
3295         result[result.length-intNum-1] ^= (1 << (n & 31));
3296 
3297         return valueOf(result);
3298     }
3299 
3300     /**
3301      * Returns the index of the rightmost (lowest-order) one bit in this
3302      * BigInteger (the number of zero bits to the right of the rightmost
3303      * one bit).  Returns -1 if this BigInteger contains no one bits.
3304      * (Computes {@code (this == 0? -1 : log2(this & -this))}.)
3305      *
3306      * @return index of the rightmost one bit in this BigInteger.
3307      */
3308     public int getLowestSetBit() {
3309         int lsb = lowestSetBitPlusTwo - 2;
3310         if (lsb == -2) {  // lowestSetBit not initialized yet
3311             lsb = 0;
3312             if (signum == 0) {
3313                 lsb -= 1;
3314             } else {
3315                 // Search for lowest order nonzero int
3316                 int i,b;
3317                 for (i=0; (b = getInt(i)) == 0; i++)
3318                     ;
3319                 lsb += (i << 5) + Integer.numberOfTrailingZeros(b);
3320             }
3321             lowestSetBitPlusTwo = lsb + 2;
3322         }
3323         return lsb;
3324     }
3325 
3326 
3327     // Miscellaneous Bit Operations
3328 
3329     /**
3330      * Returns the number of bits in the minimal two's-complement
3331      * representation of this BigInteger, <i>excluding</i> a sign bit.
3332      * For positive BigIntegers, this is equivalent to the number of bits in
3333      * the ordinary binary representation.  (Computes
3334      * {@code (ceil(log2(this < 0 ? -this : this+1)))}.)
3335      *
3336      * @return number of bits in the minimal two's-complement
3337      *         representation of this BigInteger, <i>excluding</i> a sign bit.
3338      */
3339     public int bitLength() {
3340         int n = bitLengthPlusOne - 1;
3341         if (n == -1) { // bitLength not initialized yet
3342             int[] m = mag;
3343             int len = m.length;
3344             if (len == 0) {
3345                 n = 0; // offset by one to initialize
3346             }  else {
3347                 // Calculate the bit length of the magnitude
3348                 int magBitLength = ((len - 1) << 5) + bitLengthForInt(mag[0]);
3349                  if (signum < 0) {
3350                      // Check if magnitude is a power of two
3351                      boolean pow2 = (Integer.bitCount(mag[0]) == 1);
3352                      for (int i=1; i< len && pow2; i++)
3353                          pow2 = (mag[i] == 0);
3354 
3355                      n = (pow2 ? magBitLength -1 : magBitLength);
3356                  } else {
3357                      n = magBitLength;
3358                  }
3359             }
3360             bitLengthPlusOne = n + 1;
3361         }
3362         return n;
3363     }
3364 
3365     /**
3366      * Returns the number of bits in the two's complement representation
3367      * of this BigInteger that differ from its sign bit.  This method is
3368      * useful when implementing bit-vector style sets atop BigIntegers.
3369      *
3370      * @return number of bits in the two's complement representation
3371      *         of this BigInteger that differ from its sign bit.
3372      */
3373     public int bitCount() {
3374         int bc = bitCountPlusOne - 1;
3375         if (bc == -1) {  // bitCount not initialized yet
3376             bc = 0;      // offset by one to initialize
3377             // Count the bits in the magnitude
3378             for (int i=0; i < mag.length; i++)
3379                 bc += Integer.bitCount(mag[i]);
3380             if (signum < 0) {
3381                 // Count the trailing zeros in the magnitude
3382                 int magTrailingZeroCount = 0, j;
3383                 for (j=mag.length-1; mag[j] == 0; j--)
3384                     magTrailingZeroCount += 32;
3385                 magTrailingZeroCount += Integer.numberOfTrailingZeros(mag[j]);
3386                 bc += magTrailingZeroCount - 1;
3387             }
3388             bitCountPlusOne = bc + 1;
3389         }
3390         return bc;
3391     }
3392 
3393     // Primality Testing
3394 
3395     /**
3396      * Returns {@code true} if this BigInteger is probably prime,
3397      * {@code false} if it's definitely composite.  If
3398      * {@code certainty} is &le; 0, {@code true} is
3399      * returned.
3400      *
3401      * @param  certainty a measure of the uncertainty that the caller is
3402      *         willing to tolerate: if the call returns {@code true}
3403      *         the probability that this BigInteger is prime exceeds
3404      *         (1 - 1/2<sup>{@code certainty}</sup>).  The execution time of
3405      *         this method is proportional to the value of this parameter.
3406      * @return {@code true} if this BigInteger is probably prime,
3407      *         {@code false} if it's definitely composite.
3408      */
3409     public boolean isProbablePrime(int certainty) {
3410         if (certainty <= 0)
3411             return true;
3412         BigInteger w = this.abs();
3413         if (w.equals(TWO))
3414             return true;
3415         if (!w.testBit(0) || w.equals(ONE))
3416             return false;
3417 
3418         return w.primeToCertainty(certainty, null);
3419     }
3420 
3421     // Comparison Operations
3422 
3423     /**
3424      * Compares this BigInteger with the specified BigInteger.  This
3425      * method is provided in preference to individual methods for each
3426      * of the six boolean comparison operators ({@literal <}, ==,
3427      * {@literal >}, {@literal >=}, !=, {@literal <=}).  The suggested
3428      * idiom for performing these comparisons is: {@code
3429      * (x.compareTo(y)} &lt;<i>op</i>&gt; {@code 0)}, where
3430      * &lt;<i>op</i>&gt; is one of the six comparison operators.
3431      *
3432      * @param  val BigInteger to which this BigInteger is to be compared.
3433      * @return -1, 0 or 1 as this BigInteger is numerically less than, equal
3434      *         to, or greater than {@code val}.
3435      */
3436     public int compareTo(BigInteger val) {
3437         if (signum == val.signum) {
3438             switch (signum) {
3439             case 1:
3440                 return compareMagnitude(val);
3441             case -1:
3442                 return val.compareMagnitude(this);
3443             default:
3444                 return 0;
3445             }
3446         }
3447         return signum > val.signum ? 1 : -1;
3448     }
3449 
3450     /**
3451      * Compares the magnitude array of this BigInteger with the specified
3452      * BigInteger's. This is the version of compareTo ignoring sign.
3453      *
3454      * @param val BigInteger whose magnitude array to be compared.
3455      * @return -1, 0 or 1 as this magnitude array is less than, equal to or
3456      *         greater than the magnitude aray for the specified BigInteger's.
3457      */
3458     final int compareMagnitude(BigInteger val) {
3459         int[] m1 = mag;
3460         int len1 = m1.length;
3461         int[] m2 = val.mag;
3462         int len2 = m2.length;
3463         if (len1 < len2)
3464             return -1;
3465         if (len1 > len2)
3466             return 1;
3467         for (int i = 0; i < len1; i++) {
3468             int a = m1[i];
3469             int b = m2[i];
3470             if (a != b)
3471                 return ((a & LONG_MASK) < (b & LONG_MASK)) ? -1 : 1;
3472         }
3473         return 0;
3474     }
3475 
3476     /**
3477      * Version of compareMagnitude that compares magnitude with long value.
3478      * val can't be Long.MIN_VALUE.
3479      */
3480     final int compareMagnitude(long val) {
3481         assert val != Long.MIN_VALUE;
3482         int[] m1 = mag;
3483         int len = m1.length;
3484         if (len > 2) {
3485             return 1;
3486         }
3487         if (val < 0) {
3488             val = -val;
3489         }
3490         int highWord = (int)(val >>> 32);
3491         if (highWord == 0) {
3492             if (len < 1)
3493                 return -1;
3494             if (len > 1)
3495                 return 1;
3496             int a = m1[0];
3497             int b = (int)val;
3498             if (a != b) {
3499                 return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
3500             }
3501             return 0;
3502         } else {
3503             if (len < 2)
3504                 return -1;
3505             int a = m1[0];
3506             int b = highWord;
3507             if (a != b) {
3508                 return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
3509             }
3510             a = m1[1];
3511             b = (int)val;
3512             if (a != b) {
3513                 return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
3514             }
3515             return 0;
3516         }
3517     }
3518 
3519     /**
3520      * Compares this BigInteger with the specified Object for equality.
3521      *
3522      * @param  x Object to which this BigInteger is to be compared.
3523      * @return {@code true} if and only if the specified Object is a
3524      *         BigInteger whose value is numerically equal to this BigInteger.
3525      */
3526     public boolean equals(Object x) {
3527         // This test is just an optimization, which may or may not help
3528         if (x == this)
3529             return true;
3530 
3531         if (!(x instanceof BigInteger))
3532             return false;
3533 
3534         BigInteger xInt = (BigInteger) x;
3535         if (xInt.signum != signum)
3536             return false;
3537 
3538         int[] m = mag;
3539         int len = m.length;
3540         int[] xm = xInt.mag;
3541         if (len != xm.length)
3542             return false;
3543 
3544         for (int i = 0; i < len; i++)
3545             if (xm[i] != m[i])
3546                 return false;
3547 
3548         return true;
3549     }
3550 
3551     /**
3552      * Returns the minimum of this BigInteger and {@code val}.
3553      *
3554      * @param  val value with which the minimum is to be computed.
3555      * @return the BigInteger whose value is the lesser of this BigInteger and
3556      *         {@code val}.  If they are equal, either may be returned.
3557      */
3558     public BigInteger min(BigInteger val) {
3559         return (compareTo(val) < 0 ? this : val);
3560     }
3561 
3562     /**
3563      * Returns the maximum of this BigInteger and {@code val}.
3564      *
3565      * @param  val value with which the maximum is to be computed.
3566      * @return the BigInteger whose value is the greater of this and
3567      *         {@code val}.  If they are equal, either may be returned.
3568      */
3569     public BigInteger max(BigInteger val) {
3570         return (compareTo(val) > 0 ? this : val);
3571     }
3572 
3573 
3574     // Hash Function
3575 
3576     /**
3577      * Returns the hash code for this BigInteger.
3578      *
3579      * @return hash code for this BigInteger.
3580      */
3581     public int hashCode() {
3582         int hashCode = 0;
3583 
3584         for (int i=0; i < mag.length; i++)
3585             hashCode = (int)(31*hashCode + (mag[i] & LONG_MASK));
3586 
3587         return hashCode * signum;
3588     }
3589 
3590     /**
3591      * Returns the String representation of this BigInteger in the
3592      * given radix.  If the radix is outside the range from {@link
3593      * Character#MIN_RADIX} to {@link Character#MAX_RADIX} inclusive,
3594      * it will default to 10 (as is the case for
3595      * {@code Integer.toString}).  The digit-to-character mapping
3596      * provided by {@code Character.forDigit} is used, and a minus
3597      * sign is prepended if appropriate.  (This representation is
3598      * compatible with the {@link #BigInteger(String, int) (String,
3599      * int)} constructor.)
3600      *
3601      * @param  radix  radix of the String representation.
3602      * @return String representation of this BigInteger in the given radix.
3603      * @see    Integer#toString
3604      * @see    Character#forDigit
3605      * @see    #BigInteger(java.lang.String, int)
3606      */
3607     public String toString(int radix) {
3608         if (signum == 0)
3609             return "0";
3610         if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
3611             radix = 10;
3612 
3613         // If it's small enough, use smallToString.
3614         if (mag.length <= SCHOENHAGE_BASE_CONVERSION_THRESHOLD)
3615            return smallToString(radix);
3616 
3617         // Otherwise use recursive toString, which requires positive arguments.
3618         // The results will be concatenated into this StringBuilder
3619         StringBuilder sb = new StringBuilder();
3620         if (signum < 0) {
3621             toString(this.negate(), sb, radix, 0);
3622             sb.insert(0, '-');
3623         }
3624         else
3625             toString(this, sb, radix, 0);
3626 
3627         return sb.toString();
3628     }
3629 
3630     /** This method is used to perform toString when arguments are small. */
3631     private String smallToString(int radix) {
3632         if (signum == 0) {
3633             return "0";
3634         }
3635 
3636         // Compute upper bound on number of digit groups and allocate space
3637         int maxNumDigitGroups = (4*mag.length + 6)/7;
3638         String digitGroup[] = new String[maxNumDigitGroups];
3639 
3640         // Translate number to string, a digit group at a time
3641         BigInteger tmp = this.abs();
3642         int numGroups = 0;
3643         while (tmp.signum != 0) {
3644             BigInteger d = longRadix[radix];
3645 
3646             MutableBigInteger q = new MutableBigInteger(),
3647                               a = new MutableBigInteger(tmp.mag),
3648                               b = new MutableBigInteger(d.mag);
3649             MutableBigInteger r = a.divide(b, q);
3650             BigInteger q2 = q.toBigInteger(tmp.signum * d.signum);
3651             BigInteger r2 = r.toBigInteger(tmp.signum * d.signum);
3652 
3653             digitGroup[numGroups++] = Long.toString(r2.longValue(), radix);
3654             tmp = q2;
3655         }
3656 
3657         // Put sign (if any) and first digit group into result buffer
3658         StringBuilder buf = new StringBuilder(numGroups*digitsPerLong[radix]+1);
3659         if (signum < 0) {
3660             buf.append('-');
3661         }
3662         buf.append(digitGroup[numGroups-1]);
3663 
3664         // Append remaining digit groups padded with leading zeros
3665         for (int i=numGroups-2; i >= 0; i--) {
3666             // Prepend (any) leading zeros for this digit group
3667             int numLeadingZeros = digitsPerLong[radix]-digitGroup[i].length();
3668             if (numLeadingZeros != 0) {
3669                 buf.append(zeros[numLeadingZeros]);
3670             }
3671             buf.append(digitGroup[i]);
3672         }
3673         return buf.toString();
3674     }
3675 
3676     /**
3677      * Converts the specified BigInteger to a string and appends to
3678      * {@code sb}.  This implements the recursive Schoenhage algorithm
3679      * for base conversions.
3680      * <p>
3681      * See Knuth, Donald,  _The Art of Computer Programming_, Vol. 2,
3682      * Answers to Exercises (4.4) Question 14.
3683      *
3684      * @param u      The number to convert to a string.
3685      * @param sb     The StringBuilder that will be appended to in place.
3686      * @param radix  The base to convert to.
3687      * @param digits The minimum number of digits to pad to.
3688      */
3689     private static void toString(BigInteger u, StringBuilder sb, int radix,
3690                                  int digits) {
3691         // If we're smaller than a certain threshold, use the smallToString
3692         // method, padding with leading zeroes when necessary.
3693         if (u.mag.length <= SCHOENHAGE_BASE_CONVERSION_THRESHOLD) {
3694             String s = u.smallToString(radix);
3695 
3696             // Pad with internal zeros if necessary.
3697             // Don't pad if we're at the beginning of the string.
3698             if ((s.length() < digits) && (sb.length() > 0)) {
3699                 for (int i=s.length(); i < digits; i++) {
3700                     sb.append('0');
3701                 }
3702             }
3703 
3704             sb.append(s);
3705             return;
3706         }
3707 
3708         int b, n;
3709         b = u.bitLength();
3710 
3711         // Calculate a value for n in the equation radix^(2^n) = u
3712         // and subtract 1 from that value.  This is used to find the
3713         // cache index that contains the best value to divide u.
3714         n = (int) Math.round(Math.log(b * LOG_TWO / logCache[radix]) / LOG_TWO - 1.0);
3715         BigInteger v = getRadixConversionCache(radix, n);
3716         BigInteger[] results;
3717         results = u.divideAndRemainder(v);
3718 
3719         int expectedDigits = 1 << n;
3720 
3721         // Now recursively build the two halves of each number.
3722         toString(results[0], sb, radix, digits-expectedDigits);
3723         toString(results[1], sb, radix, expectedDigits);
3724     }
3725 
3726     /**
3727      * Returns the value radix^(2^exponent) from the cache.
3728      * If this value doesn't already exist in the cache, it is added.
3729      * <p>
3730      * This could be changed to a more complicated caching method using
3731      * {@code Future}.
3732      */
3733     private static BigInteger getRadixConversionCache(int radix, int exponent) {
3734         BigInteger[] cacheLine = powerCache[radix]; // volatile read
3735         if (exponent < cacheLine.length) {
3736             return cacheLine[exponent];
3737         }
3738 
3739         int oldLength = cacheLine.length;
3740         cacheLine = Arrays.copyOf(cacheLine, exponent + 1);
3741         for (int i = oldLength; i <= exponent; i++) {
3742             cacheLine[i] = cacheLine[i - 1].pow(2);
3743         }
3744 
3745         BigInteger[][] pc = powerCache; // volatile read again
3746         if (exponent >= pc[radix].length) {
3747             pc = pc.clone();
3748             pc[radix] = cacheLine;
3749             powerCache = pc; // volatile write, publish
3750         }
3751         return cacheLine[exponent];
3752     }
3753 
3754     /* zero[i] is a string of i consecutive zeros. */
3755     private static String zeros[] = new String[64];
3756     static {
3757         zeros[63] =
3758             "000000000000000000000000000000000000000000000000000000000000000";
3759         for (int i=0; i < 63; i++)
3760             zeros[i] = zeros[63].substring(0, i);
3761     }
3762 
3763     /**
3764      * Returns the decimal String representation of this BigInteger.
3765      * The digit-to-character mapping provided by
3766      * {@code Character.forDigit} is used, and a minus sign is
3767      * prepended if appropriate.  (This representation is compatible
3768      * with the {@link #BigInteger(String) (String)} constructor, and
3769      * allows for String concatenation with Java's + operator.)
3770      *
3771      * @return decimal String representation of this BigInteger.
3772      * @see    Character#forDigit
3773      * @see    #BigInteger(java.lang.String)
3774      */
3775     public String toString() {
3776         return toString(10);
3777     }
3778 
3779     /**
3780      * Returns a byte array containing the two's-complement
3781      * representation of this BigInteger.  The byte array will be in
3782      * <i>big-endian</i> byte-order: the most significant byte is in
3783      * the zeroth element.  The array will contain the minimum number
3784      * of bytes required to represent this BigInteger, including at
3785      * least one sign bit, which is {@code (ceil((this.bitLength() +
3786      * 1)/8))}.  (This representation is compatible with the
3787      * {@link #BigInteger(byte[]) (byte[])} constructor.)
3788      *
3789      * @return a byte array containing the two's-complement representation of
3790      *         this BigInteger.
3791      * @see    #BigInteger(byte[])
3792      */
3793     public byte[] toByteArray() {
3794         int byteLen = bitLength()/8 + 1;
3795         byte[] byteArray = new byte[byteLen];
3796 
3797         for (int i=byteLen-1, bytesCopied=4, nextInt=0, intIndex=0; i >= 0; i--) {
3798             if (bytesCopied == 4) {
3799                 nextInt = getInt(intIndex++);
3800                 bytesCopied = 1;
3801             } else {
3802                 nextInt >>>= 8;
3803                 bytesCopied++;
3804             }
3805             byteArray[i] = (byte)nextInt;
3806         }
3807         return byteArray;
3808     }
3809 
3810     /**
3811      * Converts this BigInteger to an {@code int}.  This
3812      * conversion is analogous to a
3813      * <i>narrowing primitive conversion</i> from {@code long} to
3814      * {@code int} as defined in section 5.1.3 of
3815      * <cite>The Java&trade; Language Specification</cite>:
3816      * if this BigInteger is too big to fit in an
3817      * {@code int}, only the low-order 32 bits are returned.
3818      * Note that this conversion can lose information about the
3819      * overall magnitude of the BigInteger value as well as return a
3820      * result with the opposite sign.
3821      *
3822      * @return this BigInteger converted to an {@code int}.
3823      * @see #intValueExact()
3824      */
3825     public int intValue() {
3826         int result = 0;
3827         result = getInt(0);
3828         return result;
3829     }
3830 
3831     /**
3832      * Converts this BigInteger to a {@code long}.  This
3833      * conversion is analogous to a
3834      * <i>narrowing primitive conversion</i> from {@code long} to
3835      * {@code int} as defined in section 5.1.3 of
3836      * <cite>The Java&trade; Language Specification</cite>:
3837      * if this BigInteger is too big to fit in a
3838      * {@code long}, only the low-order 64 bits are returned.
3839      * Note that this conversion can lose information about the
3840      * overall magnitude of the BigInteger value as well as return a
3841      * result with the opposite sign.
3842      *
3843      * @return this BigInteger converted to a {@code long}.
3844      * @see #longValueExact()
3845      */
3846     public long longValue() {
3847         long result = 0;
3848 
3849         for (int i=1; i >= 0; i--)
3850             result = (result << 32) + (getInt(i) & LONG_MASK);
3851         return result;
3852     }
3853 
3854     /**
3855      * Converts this BigInteger to a {@code float}.  This
3856      * conversion is similar to the
3857      * <i>narrowing primitive conversion</i> from {@code double} to
3858      * {@code float} as defined in section 5.1.3 of
3859      * <cite>The Java&trade; Language Specification</cite>:
3860      * if this BigInteger has too great a magnitude
3861      * to represent as a {@code float}, it will be converted to
3862      * {@link Float#NEGATIVE_INFINITY} or {@link
3863      * Float#POSITIVE_INFINITY} as appropriate.  Note that even when
3864      * the return value is finite, this conversion can lose
3865      * information about the precision of the BigInteger value.
3866      *
3867      * @return this BigInteger converted to a {@code float}.
3868      */
3869     public float floatValue() {
3870         if (signum == 0) {
3871             return 0.0f;
3872         }
3873 
3874         int exponent = ((mag.length - 1) << 5) + bitLengthForInt(mag[0]) - 1;
3875 
3876         // exponent == floor(log2(abs(this)))
3877         if (exponent < Long.SIZE - 1) {
3878             return longValue();
3879         } else if (exponent > Float.MAX_EXPONENT) {
3880             return signum > 0 ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3881         }
3882 
3883         /*
3884          * We need the top SIGNIFICAND_WIDTH bits, including the "implicit"
3885          * one bit. To make rounding easier, we pick out the top
3886          * SIGNIFICAND_WIDTH + 1 bits, so we have one to help us round up or
3887          * down. twiceSignifFloor will contain the top SIGNIFICAND_WIDTH + 1
3888          * bits, and signifFloor the top SIGNIFICAND_WIDTH.
3889          *
3890          * It helps to consider the real number signif = abs(this) *
3891          * 2^(SIGNIFICAND_WIDTH - 1 - exponent).
3892          */
3893         int shift = exponent - FloatConsts.SIGNIFICAND_WIDTH;
3894 
3895         int twiceSignifFloor;
3896         // twiceSignifFloor will be == abs().shiftRight(shift).intValue()
3897         // We do the shift into an int directly to improve performance.
3898 
3899         int nBits = shift & 0x1f;
3900         int nBits2 = 32 - nBits;
3901 
3902         if (nBits == 0) {
3903             twiceSignifFloor = mag[0];
3904         } else {
3905             twiceSignifFloor = mag[0] >>> nBits;
3906             if (twiceSignifFloor == 0) {
3907                 twiceSignifFloor = (mag[0] << nBits2) | (mag[1] >>> nBits);
3908             }
3909         }
3910 
3911         int signifFloor = twiceSignifFloor >> 1;
3912         signifFloor &= FloatConsts.SIGNIF_BIT_MASK; // remove the implied bit
3913 
3914         /*
3915          * We round up if either the fractional part of signif is strictly
3916          * greater than 0.5 (which is true if the 0.5 bit is set and any lower
3917          * bit is set), or if the fractional part of signif is >= 0.5 and
3918          * signifFloor is odd (which is true if both the 0.5 bit and the 1 bit
3919          * are set). This is equivalent to the desired HALF_EVEN rounding.
3920          */
3921         boolean increment = (twiceSignifFloor & 1) != 0
3922                 && ((signifFloor & 1) != 0 || abs().getLowestSetBit() < shift);
3923         int signifRounded = increment ? signifFloor + 1 : signifFloor;
3924         int bits = ((exponent + FloatConsts.EXP_BIAS))
3925                 << (FloatConsts.SIGNIFICAND_WIDTH - 1);
3926         bits += signifRounded;
3927         /*
3928          * If signifRounded == 2^24, we'd need to set all of the significand
3929          * bits to zero and add 1 to the exponent. This is exactly the behavior
3930          * we get from just adding signifRounded to bits directly. If the
3931          * exponent is Float.MAX_EXPONENT, we round up (correctly) to
3932          * Float.POSITIVE_INFINITY.
3933          */
3934         bits |= signum & FloatConsts.SIGN_BIT_MASK;
3935         return Float.intBitsToFloat(bits);
3936     }
3937 
3938     /**
3939      * Converts this BigInteger to a {@code double}.  This
3940      * conversion is similar to the
3941      * <i>narrowing primitive conversion</i> from {@code double} to
3942      * {@code float} as defined in section 5.1.3 of
3943      * <cite>The Java&trade; Language Specification</cite>:
3944      * if this BigInteger has too great a magnitude
3945      * to represent as a {@code double}, it will be converted to
3946      * {@link Double#NEGATIVE_INFINITY} or {@link
3947      * Double#POSITIVE_INFINITY} as appropriate.  Note that even when
3948      * the return value is finite, this conversion can lose
3949      * information about the precision of the BigInteger value.
3950      *
3951      * @return this BigInteger converted to a {@code double}.
3952      */
3953     public double doubleValue() {
3954         if (signum == 0) {
3955             return 0.0;
3956         }
3957 
3958         int exponent = ((mag.length - 1) << 5) + bitLengthForInt(mag[0]) - 1;
3959 
3960         // exponent == floor(log2(abs(this))Double)
3961         if (exponent < Long.SIZE - 1) {
3962             return longValue();
3963         } else if (exponent > Double.MAX_EXPONENT) {
3964             return signum > 0 ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3965         }
3966 
3967         /*
3968          * We need the top SIGNIFICAND_WIDTH bits, including the "implicit"
3969          * one bit. To make rounding easier, we pick out the top
3970          * SIGNIFICAND_WIDTH + 1 bits, so we have one to help us round up or
3971          * down. twiceSignifFloor will contain the top SIGNIFICAND_WIDTH + 1
3972          * bits, and signifFloor the top SIGNIFICAND_WIDTH.
3973          *
3974          * It helps to consider the real number signif = abs(this) *
3975          * 2^(SIGNIFICAND_WIDTH - 1 - exponent).
3976          */
3977         int shift = exponent - DoubleConsts.SIGNIFICAND_WIDTH;
3978 
3979         long twiceSignifFloor;
3980         // twiceSignifFloor will be == abs().shiftRight(shift).longValue()
3981         // We do the shift into a long directly to improve performance.
3982 
3983         int nBits = shift & 0x1f;
3984         int nBits2 = 32 - nBits;
3985 
3986         int highBits;
3987         int lowBits;
3988         if (nBits == 0) {
3989             highBits = mag[0];
3990             lowBits = mag[1];
3991         } else {
3992             highBits = mag[0] >>> nBits;
3993             lowBits = (mag[0] << nBits2) | (mag[1] >>> nBits);
3994             if (highBits == 0) {
3995                 highBits = lowBits;
3996                 lowBits = (mag[1] << nBits2) | (mag[2] >>> nBits);
3997             }
3998         }
3999 
4000         twiceSignifFloor = ((highBits & LONG_MASK) << 32)
4001                 | (lowBits & LONG_MASK);
4002 
4003         long signifFloor = twiceSignifFloor >> 1;
4004         signifFloor &= DoubleConsts.SIGNIF_BIT_MASK; // remove the implied bit
4005 
4006         /*
4007          * We round up if either the fractional part of signif is strictly
4008          * greater than 0.5 (which is true if the 0.5 bit is set and any lower
4009          * bit is set), or if the fractional part of signif is >= 0.5 and
4010          * signifFloor is odd (which is true if both the 0.5 bit and the 1 bit
4011          * are set). This is equivalent to the desired HALF_EVEN rounding.
4012          */
4013         boolean increment = (twiceSignifFloor & 1) != 0
4014                 && ((signifFloor & 1) != 0 || abs().getLowestSetBit() < shift);
4015         long signifRounded = increment ? signifFloor + 1 : signifFloor;
4016         long bits = (long) ((exponent + DoubleConsts.EXP_BIAS))
4017                 << (DoubleConsts.SIGNIFICAND_WIDTH - 1);
4018         bits += signifRounded;
4019         /*
4020          * If signifRounded == 2^53, we'd need to set all of the significand
4021          * bits to zero and add 1 to the exponent. This is exactly the behavior
4022          * we get from just adding signifRounded to bits directly. If the
4023          * exponent is Double.MAX_EXPONENT, we round up (correctly) to
4024          * Double.POSITIVE_INFINITY.
4025          */
4026         bits |= signum & DoubleConsts.SIGN_BIT_MASK;
4027         return Double.longBitsToDouble(bits);
4028     }
4029 
4030     /**
4031      * Returns a copy of the input array stripped of any leading zero bytes.
4032      */
4033     private static int[] stripLeadingZeroInts(int val[]) {
4034         int vlen = val.length;
4035         int keep;
4036 
4037         // Find first nonzero byte
4038         for (keep = 0; keep < vlen && val[keep] == 0; keep++)
4039             ;
4040         return java.util.Arrays.copyOfRange(val, keep, vlen);
4041     }
4042 
4043     /**
4044      * Returns the input array stripped of any leading zero bytes.
4045      * Since the source is trusted the copying may be skipped.
4046      */
4047     private static int[] trustedStripLeadingZeroInts(int val[]) {
4048         int vlen = val.length;
4049         int keep;
4050 
4051         // Find first nonzero byte
4052         for (keep = 0; keep < vlen && val[keep] == 0; keep++)
4053             ;
4054         return keep == 0 ? val : java.util.Arrays.copyOfRange(val, keep, vlen);
4055     }
4056 
4057     /**
4058      * Returns a copy of the input array stripped of any leading zero bytes.
4059      */
4060     private static int[] stripLeadingZeroBytes(byte a[], int off, int len) {
4061         int indexBound = off + len;
4062         int keep;
4063 
4064         // Find first nonzero byte
4065         for (keep = off; keep < indexBound && a[keep] == 0; keep++)
4066             ;
4067 
4068         // Allocate new array and copy relevant part of input array
4069         int intLength = ((indexBound - keep) + 3) >>> 2;
4070         int[] result = new int[intLength];
4071         int b = indexBound - 1;
4072         for (int i = intLength-1; i >= 0; i--) {
4073             result[i] = a[b--] & 0xff;
4074             int bytesRemaining = b - keep + 1;
4075             int bytesToTransfer = Math.min(3, bytesRemaining);
4076             for (int j=8; j <= (bytesToTransfer << 3); j += 8)
4077                 result[i] |= ((a[b--] & 0xff) << j);
4078         }
4079         return result;
4080     }
4081 
4082     /**
4083      * Takes an array a representing a negative 2's-complement number and
4084      * returns the minimal (no leading zero bytes) unsigned whose value is -a.
4085      */
4086     private static int[] makePositive(byte a[], int off, int len) {
4087         int keep, k;
4088         int indexBound = off + len;
4089 
4090         // Find first non-sign (0xff) byte of input
4091         for (keep=off; keep < indexBound && a[keep] == -1; keep++)
4092             ;
4093 
4094 
4095         /* Allocate output array.  If all non-sign bytes are 0x00, we must
4096          * allocate space for one extra output byte. */
4097         for (k=keep; k < indexBound && a[k] == 0; k++)
4098             ;
4099 
4100         int extraByte = (k == indexBound) ? 1 : 0;
4101         int intLength = ((indexBound - keep + extraByte) + 3) >>> 2;
4102         int result[] = new int[intLength];
4103 
4104         /* Copy one's complement of input into output, leaving extra
4105          * byte (if it exists) == 0x00 */
4106         int b = indexBound - 1;
4107         for (int i = intLength-1; i >= 0; i--) {
4108             result[i] = a[b--] & 0xff;
4109             int numBytesToTransfer = Math.min(3, b-keep+1);
4110             if (numBytesToTransfer < 0)
4111                 numBytesToTransfer = 0;
4112             for (int j=8; j <= 8*numBytesToTransfer; j += 8)
4113                 result[i] |= ((a[b--] & 0xff) << j);
4114 
4115             // Mask indicates which bits must be complemented
4116             int mask = -1 >>> (8*(3-numBytesToTransfer));
4117             result[i] = ~result[i] & mask;
4118         }
4119 
4120         // Add one to one's complement to generate two's complement
4121         for (int i=result.length-1; i >= 0; i--) {
4122             result[i] = (int)((result[i] & LONG_MASK) + 1);
4123             if (result[i] != 0)
4124                 break;
4125         }
4126 
4127         return result;
4128     }
4129 
4130     /**
4131      * Takes an array a representing a negative 2's-complement number and
4132      * returns the minimal (no leading zero ints) unsigned whose value is -a.
4133      */
4134     private static int[] makePositive(int a[]) {
4135         int keep, j;
4136 
4137         // Find first non-sign (0xffffffff) int of input
4138         for (keep=0; keep < a.length && a[keep] == -1; keep++)
4139             ;
4140 
4141         /* Allocate output array.  If all non-sign ints are 0x00, we must
4142          * allocate space for one extra output int. */
4143         for (j=keep; j < a.length && a[j] == 0; j++)
4144             ;
4145         int extraInt = (j == a.length ? 1 : 0);
4146         int result[] = new int[a.length - keep + extraInt];
4147 
4148         /* Copy one's complement of input into output, leaving extra
4149          * int (if it exists) == 0x00 */
4150         for (int i = keep; i < a.length; i++)
4151             result[i - keep + extraInt] = ~a[i];
4152 
4153         // Add one to one's complement to generate two's complement
4154         for (int i=result.length-1; ++result[i] == 0; i--)
4155             ;
4156 
4157         return result;
4158     }
4159 
4160     /*
4161      * The following two arrays are used for fast String conversions.  Both
4162      * are indexed by radix.  The first is the number of digits of the given
4163      * radix that can fit in a Java long without "going negative", i.e., the
4164      * highest integer n such that radix**n < 2**63.  The second is the
4165      * "long radix" that tears each number into "long digits", each of which
4166      * consists of the number of digits in the corresponding element in
4167      * digitsPerLong (longRadix[i] = i**digitPerLong[i]).  Both arrays have
4168      * nonsense values in their 0 and 1 elements, as radixes 0 and 1 are not
4169      * used.
4170      */
4171     private static int digitsPerLong[] = {0, 0,
4172         62, 39, 31, 27, 24, 22, 20, 19, 18, 18, 17, 17, 16, 16, 15, 15, 15, 14,
4173         14, 14, 14, 13, 13, 13, 13, 13, 13, 12, 12, 12, 12, 12, 12, 12, 12};
4174 
4175     private static BigInteger longRadix[] = {null, null,
4176         valueOf(0x4000000000000000L), valueOf(0x383d9170b85ff80bL),
4177         valueOf(0x4000000000000000L), valueOf(0x6765c793fa10079dL),
4178         valueOf(0x41c21cb8e1000000L), valueOf(0x3642798750226111L),
4179         valueOf(0x1000000000000000L), valueOf(0x12bf307ae81ffd59L),
4180         valueOf( 0xde0b6b3a7640000L), valueOf(0x4d28cb56c33fa539L),
4181         valueOf(0x1eca170c00000000L), valueOf(0x780c7372621bd74dL),
4182         valueOf(0x1e39a5057d810000L), valueOf(0x5b27ac993df97701L),
4183         valueOf(0x1000000000000000L), valueOf(0x27b95e997e21d9f1L),
4184         valueOf(0x5da0e1e53c5c8000L), valueOf( 0xb16a458ef403f19L),
4185         valueOf(0x16bcc41e90000000L), valueOf(0x2d04b7fdd9c0ef49L),
4186         valueOf(0x5658597bcaa24000L), valueOf( 0x6feb266931a75b7L),
4187         valueOf( 0xc29e98000000000L), valueOf(0x14adf4b7320334b9L),
4188         valueOf(0x226ed36478bfa000L), valueOf(0x383d9170b85ff80bL),
4189         valueOf(0x5a3c23e39c000000L), valueOf( 0x4e900abb53e6b71L),
4190         valueOf( 0x7600ec618141000L), valueOf( 0xaee5720ee830681L),
4191         valueOf(0x1000000000000000L), valueOf(0x172588ad4f5f0981L),
4192         valueOf(0x211e44f7d02c1000L), valueOf(0x2ee56725f06e5c71L),
4193         valueOf(0x41c21cb8e1000000L)};
4194 
4195     /*
4196      * These two arrays are the integer analogue of above.
4197      */
4198     private static int digitsPerInt[] = {0, 0, 30, 19, 15, 13, 11,
4199         11, 10, 9, 9, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6,
4200         6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5};
4201 
4202     private static int intRadix[] = {0, 0,
4203         0x40000000, 0x4546b3db, 0x40000000, 0x48c27395, 0x159fd800,
4204         0x75db9c97, 0x40000000, 0x17179149, 0x3b9aca00, 0xcc6db61,
4205         0x19a10000, 0x309f1021, 0x57f6c100, 0xa2f1b6f,  0x10000000,
4206         0x18754571, 0x247dbc80, 0x3547667b, 0x4c4b4000, 0x6b5a6e1d,
4207         0x6c20a40,  0x8d2d931,  0xb640000,  0xe8d4a51,  0x1269ae40,
4208         0x17179149, 0x1cb91000, 0x23744899, 0x2b73a840, 0x34e63b41,
4209         0x40000000, 0x4cfa3cc1, 0x5c13d840, 0x6d91b519, 0x39aa400
4210     };
4211 
4212     /**
4213      * These routines provide access to the two's complement representation
4214      * of BigIntegers.
4215      */
4216 
4217     /**
4218      * Returns the length of the two's complement representation in ints,
4219      * including space for at least one sign bit.
4220      */
4221     private int intLength() {
4222         return (bitLength() >>> 5) + 1;
4223     }
4224 
4225     /* Returns sign bit */
4226     private int signBit() {
4227         return signum < 0 ? 1 : 0;
4228     }
4229 
4230     /* Returns an int of sign bits */
4231     private int signInt() {
4232         return signum < 0 ? -1 : 0;
4233     }
4234 
4235     /**
4236      * Returns the specified int of the little-endian two's complement
4237      * representation (int 0 is the least significant).  The int number can
4238      * be arbitrarily high (values are logically preceded by infinitely many
4239      * sign ints).
4240      */
4241     private int getInt(int n) {
4242         if (n < 0)
4243             return 0;
4244         if (n >= mag.length)
4245             return signInt();
4246 
4247         int magInt = mag[mag.length-n-1];
4248 
4249         return (signum >= 0 ? magInt :
4250                 (n <= firstNonzeroIntNum() ? -magInt : ~magInt));
4251     }
4252 
4253     /**
4254     * Returns the index of the int that contains the first nonzero int in the
4255     * little-endian binary representation of the magnitude (int 0 is the
4256     * least significant). If the magnitude is zero, return value is undefined.
4257     *
4258     * <p>Note: never used for a BigInteger with a magnitude of zero.
4259     * @see #getInt.
4260     */
4261     private int firstNonzeroIntNum() {
4262         int fn = firstNonzeroIntNumPlusTwo - 2;
4263         if (fn == -2) { // firstNonzeroIntNum not initialized yet
4264             // Search for the first nonzero int
4265             int i;
4266             int mlen = mag.length;
4267             for (i = mlen - 1; i >= 0 && mag[i] == 0; i--)
4268                 ;
4269             fn = mlen - i - 1;
4270             firstNonzeroIntNumPlusTwo = fn + 2; // offset by two to initialize
4271         }
4272         return fn;
4273     }
4274 
4275     /** use serialVersionUID from JDK 1.1. for interoperability */
4276     private static final long serialVersionUID = -8287574255936472291L;
4277 
4278     /**
4279      * Serializable fields for BigInteger.
4280      *
4281      * @serialField signum  int
4282      *              signum of this BigInteger
4283      * @serialField magnitude byte[]
4284      *              magnitude array of this BigInteger
4285      * @serialField bitCount  int
4286      *              appears in the serialized form for backward compatibility
4287      * @serialField bitLength int
4288      *              appears in the serialized form for backward compatibility
4289      * @serialField firstNonzeroByteNum int
4290      *              appears in the serialized form for backward compatibility
4291      * @serialField lowestSetBit int
4292      *              appears in the serialized form for backward compatibility
4293      */
4294     private static final ObjectStreamField[] serialPersistentFields = {
4295         new ObjectStreamField("signum", Integer.TYPE),
4296         new ObjectStreamField("magnitude", byte[].class),
4297         new ObjectStreamField("bitCount", Integer.TYPE),
4298         new ObjectStreamField("bitLength", Integer.TYPE),
4299         new ObjectStreamField("firstNonzeroByteNum", Integer.TYPE),
4300         new ObjectStreamField("lowestSetBit", Integer.TYPE)
4301         };
4302 
4303     /**
4304      * Reconstitute the {@code BigInteger} instance from a stream (that is,
4305      * deserialize it). The magnitude is read in as an array of bytes
4306      * for historical reasons, but it is converted to an array of ints
4307      * and the byte array is discarded.
4308      * Note:
4309      * The current convention is to initialize the cache fields, bitCountPlusOne,
4310      * bitLengthPlusOne and lowestSetBitPlusTwo, to 0 rather than some other
4311      * marker value. Therefore, no explicit action to set these fields needs to
4312      * be taken in readObject because those fields already have a 0 value by
4313      * default since defaultReadObject is not being used.
4314      */
4315     private void readObject(java.io.ObjectInputStream s)
4316         throws java.io.IOException, ClassNotFoundException {
4317         // prepare to read the alternate persistent fields
4318         ObjectInputStream.GetField fields = s.readFields();
4319 
4320         // Read the alternate persistent fields that we care about
4321         int sign = fields.get("signum", -2);
4322         byte[] magnitude = (byte[])fields.get("magnitude", null);
4323 
4324         // Validate signum
4325         if (sign < -1 || sign > 1) {
4326             String message = "BigInteger: Invalid signum value";
4327             if (fields.defaulted("signum"))
4328                 message = "BigInteger: Signum not present in stream";
4329             throw new java.io.StreamCorruptedException(message);
4330         }
4331         int[] mag = stripLeadingZeroBytes(magnitude, 0, magnitude.length);
4332         if ((mag.length == 0) != (sign == 0)) {
4333             String message = "BigInteger: signum-magnitude mismatch";
4334             if (fields.defaulted("magnitude"))
4335                 message = "BigInteger: Magnitude not present in stream";
4336             throw new java.io.StreamCorruptedException(message);
4337         }
4338 
4339         // Commit final fields via Unsafe
4340         UnsafeHolder.putSign(this, sign);
4341 
4342         // Calculate mag field from magnitude and discard magnitude
4343         UnsafeHolder.putMag(this, mag);
4344         if (mag.length >= MAX_MAG_LENGTH) {
4345             try {
4346                 checkRange();
4347             } catch (ArithmeticException e) {
4348                 throw new java.io.StreamCorruptedException("BigInteger: Out of the supported range");
4349             }
4350         }
4351     }
4352 
4353     // Support for resetting final fields while deserializing
4354     private static class UnsafeHolder {
4355         private static final sun.misc.Unsafe unsafe;
4356         private static final long signumOffset;
4357         private static final long magOffset;
4358         static {
4359             try {
4360                 unsafe = sun.misc.Unsafe.getUnsafe();
4361                 signumOffset = unsafe.objectFieldOffset
4362                     (BigInteger.class.getDeclaredField("signum"));
4363                 magOffset = unsafe.objectFieldOffset
4364                     (BigInteger.class.getDeclaredField("mag"));
4365             } catch (Exception ex) {
4366                 throw new ExceptionInInitializerError(ex);
4367             }
4368         }
4369 
4370         static void putSign(BigInteger bi, int sign) {
4371             unsafe.putInt(bi, signumOffset, sign);
4372         }
4373 
4374         static void putMag(BigInteger bi, int[] magnitude) {
4375             unsafe.putObject(bi, magOffset, magnitude);
4376         }
4377     }
4378 
4379     /**
4380      * Save the {@code BigInteger} instance to a stream.  The magnitude of a
4381      * {@code BigInteger} is serialized as a byte array for historical reasons.
4382      * To maintain compatibility with older implementations, the integers
4383      * -1, -1, -2, and -2 are written as the values of the obsolete fields
4384      * {@code bitCount}, {@code bitLength}, {@code lowestSetBit}, and
4385      * {@code firstNonzeroByteNum}, respectively.  These values are compatible
4386      * with older implementations, but will be ignored by current
4387      * implementations.
4388      */
4389     private void writeObject(ObjectOutputStream s) throws IOException {
4390         // set the values of the Serializable fields
4391         ObjectOutputStream.PutField fields = s.putFields();
4392         fields.put("signum", signum);
4393         fields.put("magnitude", magSerializedForm());
4394         // The values written for cached fields are compatible with older
4395         // versions, but are ignored in readObject so don't otherwise matter.
4396         fields.put("bitCount", -1);
4397         fields.put("bitLength", -1);
4398         fields.put("lowestSetBit", -2);
4399         fields.put("firstNonzeroByteNum", -2);
4400 
4401         // save them
4402         s.writeFields();
4403     }
4404 
4405     /**
4406      * Returns the mag array as an array of bytes.
4407      */
4408     private byte[] magSerializedForm() {
4409         int len = mag.length;
4410 
4411         int bitLen = (len == 0 ? 0 : ((len - 1) << 5) + bitLengthForInt(mag[0]));
4412         int byteLen = (bitLen + 7) >>> 3;
4413         byte[] result = new byte[byteLen];
4414 
4415         for (int i = byteLen - 1, bytesCopied = 4, intIndex = len - 1, nextInt = 0;
4416              i >= 0; i--) {
4417             if (bytesCopied == 4) {
4418                 nextInt = mag[intIndex--];
4419                 bytesCopied = 1;
4420             } else {
4421                 nextInt >>>= 8;
4422                 bytesCopied++;
4423             }
4424             result[i] = (byte)nextInt;
4425         }
4426         return result;
4427     }
4428 
4429     /**
4430      * Converts this {@code BigInteger} to a {@code long}, checking
4431      * for lost information.  If the value of this {@code BigInteger}
4432      * is out of the range of the {@code long} type, then an
4433      * {@code ArithmeticException} is thrown.
4434      *
4435      * @return this {@code BigInteger} converted to a {@code long}.
4436      * @throws ArithmeticException if the value of {@code this} will
4437      * not exactly fit in a {@code long}.
4438      * @see BigInteger#longValue
4439      * @since  1.8
4440      */
4441     public long longValueExact() {
4442         if (mag.length <= 2 && bitLength() <= 63)
4443             return longValue();
4444         else
4445             throw new ArithmeticException("BigInteger out of long range");
4446     }
4447 
4448     /**
4449      * Converts this {@code BigInteger} to an {@code int}, checking
4450      * for lost information.  If the value of this {@code BigInteger}
4451      * is out of the range of the {@code int} type, then an
4452      * {@code ArithmeticException} is thrown.
4453      *
4454      * @return this {@code BigInteger} converted to an {@code int}.
4455      * @throws ArithmeticException if the value of {@code this} will
4456      * not exactly fit in a {@code int}.
4457      * @see BigInteger#intValue
4458      * @since  1.8
4459      */
4460     public int intValueExact() {
4461         if (mag.length <= 1 && bitLength() <= 31)
4462             return intValue();
4463         else
4464             throw new ArithmeticException("BigInteger out of int range");
4465     }
4466 
4467     /**
4468      * Converts this {@code BigInteger} to a {@code short}, checking
4469      * for lost information.  If the value of this {@code BigInteger}
4470      * is out of the range of the {@code short} type, then an
4471      * {@code ArithmeticException} is thrown.
4472      *
4473      * @return this {@code BigInteger} converted to a {@code short}.
4474      * @throws ArithmeticException if the value of {@code this} will
4475      * not exactly fit in a {@code short}.
4476      * @see BigInteger#shortValue
4477      * @since  1.8
4478      */
4479     public short shortValueExact() {
4480         if (mag.length <= 1 && bitLength() <= 31) {
4481             int value = intValue();
4482             if (value >= Short.MIN_VALUE && value <= Short.MAX_VALUE)
4483                 return shortValue();
4484         }
4485         throw new ArithmeticException("BigInteger out of short range");
4486     }
4487 
4488     /**
4489      * Converts this {@code BigInteger} to a {@code byte}, checking
4490      * for lost information.  If the value of this {@code BigInteger}
4491      * is out of the range of the {@code byte} type, then an
4492      * {@code ArithmeticException} is thrown.
4493      *
4494      * @return this {@code BigInteger} converted to a {@code byte}.
4495      * @throws ArithmeticException if the value of {@code this} will
4496      * not exactly fit in a {@code byte}.
4497      * @see BigInteger#byteValue
4498      * @since  1.8
4499      */
4500     public byte byteValueExact() {
4501         if (mag.length <= 1 && bitLength() <= 31) {
4502             int value = intValue();
4503             if (value >= Byte.MIN_VALUE && value <= Byte.MAX_VALUE)
4504                 return byteValue();
4505         }
4506         throw new ArithmeticException("BigInteger out of byte range");
4507     }
4508 }