8247782: typos in java.math
Reviewed-by: rriggs, lancea, darcy

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