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