rev 7571 : 8017540: Improve multi-threaded contention behavior of radix conversion cache
Summary: Replace array of ArrayList of BigIntegers with a volatile two-dimensional BigInteger array eliminate the synchronization of getRadixConversionCache()
Reviewed-by: plevart, shade, bpb, alanb
Contributed-by: Peter Levart <peter.levart@gmail.com>, Dmitry Nadezhin <dmitry.nadezhin@oracle.com>, Aleksey Shipilev <aleksey.shipilev@oracle.com>

   1 /*
   2  * Copyright (c) 1996, 2013, Oracle and/or its affiliates. All rights reserved.
   3  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
   4  *
   5  * This code is free software; you can redistribute it and/or modify it
   6  * under the terms of the GNU General Public License version 2 only, as
   7  * published by the Free Software Foundation.  Oracle designates this
   8  * particular file as subject to the "Classpath" exception as provided
   9  * by Oracle in the LICENSE file that accompanied this code.
  10  *
  11  * This code is distributed in the hope that it will be useful, but WITHOUT
  12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14  * version 2 for more details (a copy is included in the LICENSE file that
  15  * accompanied this code).
  16  *
  17  * You should have received a copy of the GNU General Public License version
  18  * 2 along with this work; if not, write to the Free Software Foundation,
  19  * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
  20  *
  21  * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
  22  * or visit www.oracle.com if you need additional information or have any
  23  * questions.
  24  */
  25 
  26 /*
  27  * Portions Copyright (c) 1995  Colin Plumb.  All rights reserved.
  28  */
  29 
  30 package java.math;
  31 
  32 import java.io.IOException;
  33 import java.io.ObjectInputStream;
  34 import java.io.ObjectOutputStream;
  35 import java.io.ObjectStreamField;
  36 import java.util.ArrayList;
  37 import java.util.Arrays;
  38 import java.util.Random;
  39 import sun.misc.DoubleConsts;
  40 import sun.misc.FloatConsts;
  41 
  42 /**
  43  * Immutable arbitrary-precision integers.  All operations behave as if
  44  * BigIntegers were represented in two's-complement notation (like Java's
  45  * primitive integer types).  BigInteger provides analogues to all of Java's
  46  * primitive integer operators, and all relevant methods from java.lang.Math.
  47  * Additionally, BigInteger provides operations for modular arithmetic, GCD
  48  * calculation, primality testing, prime generation, bit manipulation,
  49  * and a few other miscellaneous operations.
  50  *
  51  * <p>Semantics of arithmetic operations exactly mimic those of Java's integer
  52  * arithmetic operators, as defined in <i>The Java Language Specification</i>.
  53  * For example, division by zero throws an {@code ArithmeticException}, and
  54  * division of a negative by a positive yields a negative (or zero) remainder.
  55  * All of the details in the Spec concerning overflow are ignored, as
  56  * BigIntegers are made as large as necessary to accommodate the results of an
  57  * operation.
  58  *
  59  * <p>Semantics of shift operations extend those of Java's shift operators
  60  * to allow for negative shift distances.  A right-shift with a negative
  61  * shift distance results in a left shift, and vice-versa.  The unsigned
  62  * right shift operator ({@code >>>}) is omitted, as this operation makes
  63  * little sense in combination with the "infinite word size" abstraction
  64  * provided by this class.
  65  *
  66  * <p>Semantics of bitwise logical operations exactly mimic those of Java's
  67  * bitwise integer operators.  The binary operators ({@code and},
  68  * {@code or}, {@code xor}) implicitly perform sign extension on the shorter
  69  * of the two operands prior to performing the operation.
  70  *
  71  * <p>Comparison operations perform signed integer comparisons, analogous to
  72  * those performed by Java's relational and equality operators.
  73  *
  74  * <p>Modular arithmetic operations are provided to compute residues, perform
  75  * exponentiation, and compute multiplicative inverses.  These methods always
  76  * return a non-negative result, between {@code 0} and {@code (modulus - 1)},
  77  * inclusive.
  78  *
  79  * <p>Bit operations operate on a single bit of the two's-complement
  80  * representation of their operand.  If necessary, the operand is sign-
  81  * extended so that it contains the designated bit.  None of the single-bit
  82  * operations can produce a BigInteger with a different sign from the
  83  * BigInteger being operated on, as they affect only a single bit, and the
  84  * "infinite word size" abstraction provided by this class ensures that there
  85  * are infinitely many "virtual sign bits" preceding each BigInteger.
  86  *
  87  * <p>For the sake of brevity and clarity, pseudo-code is used throughout the
  88  * descriptions of BigInteger methods.  The pseudo-code expression
  89  * {@code (i + j)} is shorthand for "a BigInteger whose value is
  90  * that of the BigInteger {@code i} plus that of the BigInteger {@code j}."
  91  * The pseudo-code expression {@code (i == j)} is shorthand for
  92  * "{@code true} if and only if the BigInteger {@code i} represents the same
  93  * value as the BigInteger {@code j}."  Other pseudo-code expressions are
  94  * interpreted similarly.
  95  *
  96  * <p>All methods and constructors in this class throw
  97  * {@code NullPointerException} when passed
  98  * a null object reference for any input parameter.
  99  *
 100  * @see     BigDecimal
 101  * @author  Josh Bloch
 102  * @author  Michael McCloskey
 103  * @author  Alan Eliasen
 104  * @since JDK1.1
 105  */
 106 
 107 public class BigInteger extends Number implements Comparable<BigInteger> {
 108     /**
 109      * The signum of this BigInteger: -1 for negative, 0 for zero, or
 110      * 1 for positive.  Note that the BigInteger zero <i>must</i> have
 111      * a signum of 0.  This is necessary to ensures that there is exactly one
 112      * representation for each BigInteger value.
 113      *
 114      * @serial
 115      */
 116     final int signum;
 117 
 118     /**
 119      * The magnitude of this BigInteger, in <i>big-endian</i> order: the
 120      * zeroth element of this array is the most-significant int of the
 121      * magnitude.  The magnitude must be "minimal" in that the most-significant
 122      * int ({@code mag[0]}) must be non-zero.  This is necessary to
 123      * ensure that there is exactly one representation for each BigInteger
 124      * value.  Note that this implies that the BigInteger zero has a
 125      * zero-length mag array.
 126      */
 127     final int[] mag;
 128 
 129     // These "redundant fields" are initialized with recognizable nonsense
 130     // values, and cached the first time they are needed (or never, if they
 131     // aren't needed).
 132 
 133      /**
 134      * One plus the bitCount of this BigInteger. Zeros means unitialized.
 135      *
 136      * @serial
 137      * @see #bitCount
 138      * @deprecated Deprecated since logical value is offset from stored
 139      * value and correction factor is applied in accessor method.
 140      */
 141     @Deprecated
 142     private int bitCount;
 143 
 144     /**
 145      * One plus the bitLength of this BigInteger. Zeros means unitialized.
 146      * (either value is acceptable).
 147      *
 148      * @serial
 149      * @see #bitLength()
 150      * @deprecated Deprecated since logical value is offset from stored
 151      * value and correction factor is applied in accessor method.
 152      */
 153     @Deprecated
 154     private int bitLength;
 155 
 156     /**
 157      * Two plus the lowest set bit of this BigInteger, as returned by
 158      * getLowestSetBit().
 159      *
 160      * @serial
 161      * @see #getLowestSetBit
 162      * @deprecated Deprecated since logical value is offset from stored
 163      * value and correction factor is applied in accessor method.
 164      */
 165     @Deprecated
 166     private int lowestSetBit;
 167 
 168     /**
 169      * Two plus the index of the lowest-order int in the magnitude of this
 170      * BigInteger that contains a nonzero int, or -2 (either value is acceptable).
 171      * The least significant int has int-number 0, the next int in order of
 172      * increasing significance has int-number 1, and so forth.
 173      * @deprecated Deprecated since logical value is offset from stored
 174      * value and correction factor is applied in accessor method.
 175      */
 176     @Deprecated
 177     private int firstNonzeroIntNum;
 178 
 179     /**
 180      * This mask is used to obtain the value of an int as if it were unsigned.
 181      */
 182     final static long LONG_MASK = 0xffffffffL;
 183 
 184     /**
 185      * The threshold value for using Karatsuba multiplication.  If the number
 186      * of ints in both mag arrays are greater than this number, then
 187      * Karatsuba multiplication will be used.   This value is found
 188      * experimentally to work well.
 189      */
 190     private static final int KARATSUBA_THRESHOLD = 50;
 191 
 192     /**
 193      * The threshold value for using 3-way Toom-Cook multiplication.
 194      * If the number of ints in each mag array is greater than the
 195      * Karatsuba threshold, and the number of ints in at least one of
 196      * the mag arrays is greater than this threshold, then Toom-Cook
 197      * multiplication will be used.
 198      */
 199     private static final int TOOM_COOK_THRESHOLD = 75;
 200 
 201     /**
 202      * The threshold value for using Karatsuba squaring.  If the number
 203      * of ints in the number are larger than this value,
 204      * Karatsuba squaring will be used.   This value is found
 205      * experimentally to work well.
 206      */
 207     private static final int KARATSUBA_SQUARE_THRESHOLD = 90;
 208 
 209     /**
 210      * The threshold value for using Toom-Cook squaring.  If the number
 211      * of ints in the number are larger than this value,
 212      * Toom-Cook squaring will be used.   This value is found
 213      * experimentally to work well.
 214      */
 215     private static final int TOOM_COOK_SQUARE_THRESHOLD = 140;
 216 
 217     /**
 218      * The threshold value for using Schoenhage recursive base conversion. If
 219      * the number of ints in the number are larger than this value,
 220      * the Schoenhage algorithm will be used.  In practice, it appears that the
 221      * Schoenhage routine is faster for any threshold down to 2, and is
 222      * relatively flat for thresholds between 2-25, so this choice may be
 223      * varied within this range for very small effect.
 224      */
 225     private static final int SCHOENHAGE_BASE_CONVERSION_THRESHOLD = 8;
 226 
 227     //Constructors
 228 
 229     /**
 230      * Translates a byte array containing the two's-complement binary
 231      * representation of a BigInteger into a BigInteger.  The input array is
 232      * assumed to be in <i>big-endian</i> byte-order: the most significant
 233      * byte is in the zeroth element.
 234      *
 235      * @param  val big-endian two's-complement binary representation of
 236      *         BigInteger.
 237      * @throws NumberFormatException {@code val} is zero bytes long.
 238      */
 239     public BigInteger(byte[] val) {
 240         if (val.length == 0)
 241             throw new NumberFormatException("Zero length BigInteger");
 242 
 243         if (val[0] < 0) {
 244             mag = makePositive(val);
 245             signum = -1;
 246         } else {
 247             mag = stripLeadingZeroBytes(val);
 248             signum = (mag.length == 0 ? 0 : 1);
 249         }
 250     }
 251 
 252     /**
 253      * This private constructor translates an int array containing the
 254      * two's-complement binary representation of a BigInteger into a
 255      * BigInteger. The input array is assumed to be in <i>big-endian</i>
 256      * int-order: the most significant int is in the zeroth element.
 257      */
 258     private BigInteger(int[] val) {
 259         if (val.length == 0)
 260             throw new NumberFormatException("Zero length BigInteger");
 261 
 262         if (val[0] < 0) {
 263             mag = makePositive(val);
 264             signum = -1;
 265         } else {
 266             mag = trustedStripLeadingZeroInts(val);
 267             signum = (mag.length == 0 ? 0 : 1);
 268         }
 269     }
 270 
 271     /**
 272      * Translates the sign-magnitude representation of a BigInteger into a
 273      * BigInteger.  The sign is represented as an integer signum value: -1 for
 274      * negative, 0 for zero, or 1 for positive.  The magnitude is a byte array
 275      * in <i>big-endian</i> byte-order: the most significant byte is in the
 276      * zeroth element.  A zero-length magnitude array is permissible, and will
 277      * result in a BigInteger value of 0, whether signum is -1, 0 or 1.
 278      *
 279      * @param  signum signum of the number (-1 for negative, 0 for zero, 1
 280      *         for positive).
 281      * @param  magnitude big-endian binary representation of the magnitude of
 282      *         the number.
 283      * @throws NumberFormatException {@code signum} is not one of the three
 284      *         legal values (-1, 0, and 1), or {@code signum} is 0 and
 285      *         {@code magnitude} contains one or more non-zero bytes.
 286      */
 287     public BigInteger(int signum, byte[] magnitude) {
 288         this.mag = stripLeadingZeroBytes(magnitude);
 289 
 290         if (signum < -1 || signum > 1)
 291             throw(new NumberFormatException("Invalid signum value"));
 292 
 293         if (this.mag.length==0) {
 294             this.signum = 0;
 295         } else {
 296             if (signum == 0)
 297                 throw(new NumberFormatException("signum-magnitude mismatch"));
 298             this.signum = signum;
 299         }
 300     }
 301 
 302     /**
 303      * A constructor for internal use that translates the sign-magnitude
 304      * representation of a BigInteger into a BigInteger. It checks the
 305      * arguments and copies the magnitude so this constructor would be
 306      * safe for external use.
 307      */
 308     private BigInteger(int signum, int[] magnitude) {
 309         this.mag = stripLeadingZeroInts(magnitude);
 310 
 311         if (signum < -1 || signum > 1)
 312             throw(new NumberFormatException("Invalid signum value"));
 313 
 314         if (this.mag.length==0) {
 315             this.signum = 0;
 316         } else {
 317             if (signum == 0)
 318                 throw(new NumberFormatException("signum-magnitude mismatch"));
 319             this.signum = signum;
 320         }
 321     }
 322 
 323     /**
 324      * Translates the String representation of a BigInteger in the
 325      * specified radix into a BigInteger.  The String representation
 326      * consists of an optional minus or plus sign followed by a
 327      * sequence of one or more digits in the specified radix.  The
 328      * character-to-digit mapping is provided by {@code
 329      * Character.digit}.  The String may not contain any extraneous
 330      * characters (whitespace, for example).
 331      *
 332      * @param val String representation of BigInteger.
 333      * @param radix radix to be used in interpreting {@code val}.
 334      * @throws NumberFormatException {@code val} is not a valid representation
 335      *         of a BigInteger in the specified radix, or {@code radix} is
 336      *         outside the range from {@link Character#MIN_RADIX} to
 337      *         {@link Character#MAX_RADIX}, inclusive.
 338      * @see    Character#digit
 339      */
 340     public BigInteger(String val, int radix) {
 341         int cursor = 0, numDigits;
 342         final int len = val.length();
 343 
 344         if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
 345             throw new NumberFormatException("Radix out of range");
 346         if (len == 0)
 347             throw new NumberFormatException("Zero length BigInteger");
 348 
 349         // Check for at most one leading sign
 350         int sign = 1;
 351         int index1 = val.lastIndexOf('-');
 352         int index2 = val.lastIndexOf('+');
 353         if ((index1 + index2) <= -1) {
 354             // No leading sign character or at most one leading sign character
 355             if (index1 == 0 || index2 == 0) {
 356                 cursor = 1;
 357                 if (len == 1)
 358                     throw new NumberFormatException("Zero length BigInteger");
 359             }
 360             if (index1 == 0)
 361                 sign = -1;
 362         } else
 363             throw new NumberFormatException("Illegal embedded sign character");
 364 
 365         // Skip leading zeros and compute number of digits in magnitude
 366         while (cursor < len &&
 367                Character.digit(val.charAt(cursor), radix) == 0)
 368             cursor++;
 369         if (cursor == len) {
 370             signum = 0;
 371             mag = ZERO.mag;
 372             return;
 373         }
 374 
 375         numDigits = len - cursor;
 376         signum = sign;
 377 
 378         // Pre-allocate array of expected size. May be too large but can
 379         // never be too small. Typically exact.
 380         int numBits = (int)(((numDigits * bitsPerDigit[radix]) >>> 10) + 1);
 381         int numWords = (numBits + 31) >>> 5;
 382         int[] magnitude = new int[numWords];
 383 
 384         // Process first (potentially short) digit group
 385         int firstGroupLen = numDigits % digitsPerInt[radix];
 386         if (firstGroupLen == 0)
 387             firstGroupLen = digitsPerInt[radix];
 388         String group = val.substring(cursor, cursor += firstGroupLen);
 389         magnitude[numWords - 1] = Integer.parseInt(group, radix);
 390         if (magnitude[numWords - 1] < 0)
 391             throw new NumberFormatException("Illegal digit");
 392 
 393         // Process remaining digit groups
 394         int superRadix = intRadix[radix];
 395         int groupVal = 0;
 396         while (cursor < len) {
 397             group = val.substring(cursor, cursor += digitsPerInt[radix]);
 398             groupVal = Integer.parseInt(group, radix);
 399             if (groupVal < 0)
 400                 throw new NumberFormatException("Illegal digit");
 401             destructiveMulAdd(magnitude, superRadix, groupVal);
 402         }
 403         // Required for cases where the array was overallocated.
 404         mag = trustedStripLeadingZeroInts(magnitude);
 405     }
 406 
 407     /*
 408      * Constructs a new BigInteger using a char array with radix=10.
 409      * Sign is precalculated outside and not allowed in the val.
 410      */
 411     BigInteger(char[] val, int sign, int len) {
 412         int cursor = 0, numDigits;
 413 
 414         // Skip leading zeros and compute number of digits in magnitude
 415         while (cursor < len && Character.digit(val[cursor], 10) == 0) {
 416             cursor++;
 417         }
 418         if (cursor == len) {
 419             signum = 0;
 420             mag = ZERO.mag;
 421             return;
 422         }
 423 
 424         numDigits = len - cursor;
 425         signum = sign;
 426         // Pre-allocate array of expected size
 427         int numWords;
 428         if (len < 10) {
 429             numWords = 1;
 430         } else {
 431             int numBits = (int)(((numDigits * bitsPerDigit[10]) >>> 10) + 1);
 432             numWords = (numBits + 31) >>> 5;
 433         }
 434         int[] magnitude = new int[numWords];
 435 
 436         // Process first (potentially short) digit group
 437         int firstGroupLen = numDigits % digitsPerInt[10];
 438         if (firstGroupLen == 0)
 439             firstGroupLen = digitsPerInt[10];
 440         magnitude[numWords - 1] = parseInt(val, cursor,  cursor += firstGroupLen);
 441 
 442         // Process remaining digit groups
 443         while (cursor < len) {
 444             int groupVal = parseInt(val, cursor, cursor += digitsPerInt[10]);
 445             destructiveMulAdd(magnitude, intRadix[10], groupVal);
 446         }
 447         mag = trustedStripLeadingZeroInts(magnitude);
 448     }
 449 
 450     // Create an integer with the digits between the two indexes
 451     // Assumes start < end. The result may be negative, but it
 452     // is to be treated as an unsigned value.
 453     private int parseInt(char[] source, int start, int end) {
 454         int result = Character.digit(source[start++], 10);
 455         if (result == -1)
 456             throw new NumberFormatException(new String(source));
 457 
 458         for (int index = start; index<end; index++) {
 459             int nextVal = Character.digit(source[index], 10);
 460             if (nextVal == -1)
 461                 throw new NumberFormatException(new String(source));
 462             result = 10*result + nextVal;
 463         }
 464 
 465         return result;
 466     }
 467 
 468     // bitsPerDigit in the given radix times 1024
 469     // Rounded up to avoid underallocation.
 470     private static long bitsPerDigit[] = { 0, 0,
 471         1024, 1624, 2048, 2378, 2648, 2875, 3072, 3247, 3402, 3543, 3672,
 472         3790, 3899, 4001, 4096, 4186, 4271, 4350, 4426, 4498, 4567, 4633,
 473         4696, 4756, 4814, 4870, 4923, 4975, 5025, 5074, 5120, 5166, 5210,
 474                                            5253, 5295};
 475 
 476     // Multiply x array times word y in place, and add word z
 477     private static void destructiveMulAdd(int[] x, int y, int z) {
 478         // Perform the multiplication word by word
 479         long ylong = y & LONG_MASK;
 480         long zlong = z & LONG_MASK;
 481         int len = x.length;
 482 
 483         long product = 0;
 484         long carry = 0;
 485         for (int i = len-1; i >= 0; i--) {
 486             product = ylong * (x[i] & LONG_MASK) + carry;
 487             x[i] = (int)product;
 488             carry = product >>> 32;
 489         }
 490 
 491         // Perform the addition
 492         long sum = (x[len-1] & LONG_MASK) + zlong;
 493         x[len-1] = (int)sum;
 494         carry = sum >>> 32;
 495         for (int i = len-2; i >= 0; i--) {
 496             sum = (x[i] & LONG_MASK) + carry;
 497             x[i] = (int)sum;
 498             carry = sum >>> 32;
 499         }
 500     }
 501 
 502     /**
 503      * Translates the decimal String representation of a BigInteger into a
 504      * BigInteger.  The String representation consists of an optional minus
 505      * sign followed by a sequence of one or more decimal digits.  The
 506      * character-to-digit mapping is provided by {@code Character.digit}.
 507      * The String may not contain any extraneous characters (whitespace, for
 508      * example).
 509      *
 510      * @param val decimal String representation of BigInteger.
 511      * @throws NumberFormatException {@code val} is not a valid representation
 512      *         of a BigInteger.
 513      * @see    Character#digit
 514      */
 515     public BigInteger(String val) {
 516         this(val, 10);
 517     }
 518 
 519     /**
 520      * Constructs a randomly generated BigInteger, uniformly distributed over
 521      * the range 0 to (2<sup>{@code numBits}</sup> - 1), inclusive.
 522      * The uniformity of the distribution assumes that a fair source of random
 523      * bits is provided in {@code rnd}.  Note that this constructor always
 524      * constructs a non-negative BigInteger.
 525      *
 526      * @param  numBits maximum bitLength of the new BigInteger.
 527      * @param  rnd source of randomness to be used in computing the new
 528      *         BigInteger.
 529      * @throws IllegalArgumentException {@code numBits} is negative.
 530      * @see #bitLength()
 531      */
 532     public BigInteger(int numBits, Random rnd) {
 533         this(1, randomBits(numBits, rnd));
 534     }
 535 
 536     private static byte[] randomBits(int numBits, Random rnd) {
 537         if (numBits < 0)
 538             throw new IllegalArgumentException("numBits must be non-negative");
 539         int numBytes = (int)(((long)numBits+7)/8); // avoid overflow
 540         byte[] randomBits = new byte[numBytes];
 541 
 542         // Generate random bytes and mask out any excess bits
 543         if (numBytes > 0) {
 544             rnd.nextBytes(randomBits);
 545             int excessBits = 8*numBytes - numBits;
 546             randomBits[0] &= (1 << (8-excessBits)) - 1;
 547         }
 548         return randomBits;
 549     }
 550 
 551     /**
 552      * Constructs a randomly generated positive BigInteger that is probably
 553      * prime, with the specified bitLength.
 554      *
 555      * <p>It is recommended that the {@link #probablePrime probablePrime}
 556      * method be used in preference to this constructor unless there
 557      * is a compelling need to specify a certainty.
 558      *
 559      * @param  bitLength bitLength of the returned BigInteger.
 560      * @param  certainty a measure of the uncertainty that the caller is
 561      *         willing to tolerate.  The probability that the new BigInteger
 562      *         represents a prime number will exceed
 563      *         (1 - 1/2<sup>{@code certainty}</sup>).  The execution time of
 564      *         this constructor is proportional to the value of this parameter.
 565      * @param  rnd source of random bits used to select candidates to be
 566      *         tested for primality.
 567      * @throws ArithmeticException {@code bitLength < 2}.
 568      * @see    #bitLength()
 569      */
 570     public BigInteger(int bitLength, int certainty, Random rnd) {
 571         BigInteger prime;
 572 
 573         if (bitLength < 2)
 574             throw new ArithmeticException("bitLength < 2");
 575         prime = (bitLength < SMALL_PRIME_THRESHOLD
 576                                 ? smallPrime(bitLength, certainty, rnd)
 577                                 : largePrime(bitLength, certainty, rnd));
 578         signum = 1;
 579         mag = prime.mag;
 580     }
 581 
 582     // Minimum size in bits that the requested prime number has
 583     // before we use the large prime number generating algorithms.
 584     // The cutoff of 95 was chosen empirically for best performance.
 585     private static final int SMALL_PRIME_THRESHOLD = 95;
 586 
 587     // Certainty required to meet the spec of probablePrime
 588     private static final int DEFAULT_PRIME_CERTAINTY = 100;
 589 
 590     /**
 591      * Returns a positive BigInteger that is probably prime, with the
 592      * specified bitLength. The probability that a BigInteger returned
 593      * by this method is composite does not exceed 2<sup>-100</sup>.
 594      *
 595      * @param  bitLength bitLength of the returned BigInteger.
 596      * @param  rnd source of random bits used to select candidates to be
 597      *         tested for primality.
 598      * @return a BigInteger of {@code bitLength} bits that is probably prime
 599      * @throws ArithmeticException {@code bitLength < 2}.
 600      * @see    #bitLength()
 601      * @since 1.4
 602      */
 603     public static BigInteger probablePrime(int bitLength, Random rnd) {
 604         if (bitLength < 2)
 605             throw new ArithmeticException("bitLength < 2");
 606 
 607         return (bitLength < SMALL_PRIME_THRESHOLD ?
 608                 smallPrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd) :
 609                 largePrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd));
 610     }
 611 
 612     /**
 613      * Find a random number of the specified bitLength that is probably prime.
 614      * This method is used for smaller primes, its performance degrades on
 615      * larger bitlengths.
 616      *
 617      * This method assumes bitLength > 1.
 618      */
 619     private static BigInteger smallPrime(int bitLength, int certainty, Random rnd) {
 620         int magLen = (bitLength + 31) >>> 5;
 621         int temp[] = new int[magLen];
 622         int highBit = 1 << ((bitLength+31) & 0x1f);  // High bit of high int
 623         int highMask = (highBit << 1) - 1;  // Bits to keep in high int
 624 
 625         while(true) {
 626             // Construct a candidate
 627             for (int i=0; i<magLen; i++)
 628                 temp[i] = rnd.nextInt();
 629             temp[0] = (temp[0] & highMask) | highBit;  // Ensure exact length
 630             if (bitLength > 2)
 631                 temp[magLen-1] |= 1;  // Make odd if bitlen > 2
 632 
 633             BigInteger p = new BigInteger(temp, 1);
 634 
 635             // Do cheap "pre-test" if applicable
 636             if (bitLength > 6) {
 637                 long r = p.remainder(SMALL_PRIME_PRODUCT).longValue();
 638                 if ((r%3==0)  || (r%5==0)  || (r%7==0)  || (r%11==0) ||
 639                     (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) ||
 640                     (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0))
 641                     continue; // Candidate is composite; try another
 642             }
 643 
 644             // All candidates of bitLength 2 and 3 are prime by this point
 645             if (bitLength < 4)
 646                 return p;
 647 
 648             // Do expensive test if we survive pre-test (or it's inapplicable)
 649             if (p.primeToCertainty(certainty, rnd))
 650                 return p;
 651         }
 652     }
 653 
 654     private static final BigInteger SMALL_PRIME_PRODUCT
 655                        = valueOf(3L*5*7*11*13*17*19*23*29*31*37*41);
 656 
 657     /**
 658      * Find a random number of the specified bitLength that is probably prime.
 659      * This method is more appropriate for larger bitlengths since it uses
 660      * a sieve to eliminate most composites before using a more expensive
 661      * test.
 662      */
 663     private static BigInteger largePrime(int bitLength, int certainty, Random rnd) {
 664         BigInteger p;
 665         p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
 666         p.mag[p.mag.length-1] &= 0xfffffffe;
 667 
 668         // Use a sieve length likely to contain the next prime number
 669         int searchLen = (bitLength / 20) * 64;
 670         BitSieve searchSieve = new BitSieve(p, searchLen);
 671         BigInteger candidate = searchSieve.retrieve(p, certainty, rnd);
 672 
 673         while ((candidate == null) || (candidate.bitLength() != bitLength)) {
 674             p = p.add(BigInteger.valueOf(2*searchLen));
 675             if (p.bitLength() != bitLength)
 676                 p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
 677             p.mag[p.mag.length-1] &= 0xfffffffe;
 678             searchSieve = new BitSieve(p, searchLen);
 679             candidate = searchSieve.retrieve(p, certainty, rnd);
 680         }
 681         return candidate;
 682     }
 683 
 684    /**
 685     * Returns the first integer greater than this {@code BigInteger} that
 686     * is probably prime.  The probability that the number returned by this
 687     * method is composite does not exceed 2<sup>-100</sup>. This method will
 688     * never skip over a prime when searching: if it returns {@code p}, there
 689     * is no prime {@code q} such that {@code this < q < p}.
 690     *
 691     * @return the first integer greater than this {@code BigInteger} that
 692     *         is probably prime.
 693     * @throws ArithmeticException {@code this < 0}.
 694     * @since 1.5
 695     */
 696     public BigInteger nextProbablePrime() {
 697         if (this.signum < 0)
 698             throw new ArithmeticException("start < 0: " + this);
 699 
 700         // Handle trivial cases
 701         if ((this.signum == 0) || this.equals(ONE))
 702             return TWO;
 703 
 704         BigInteger result = this.add(ONE);
 705 
 706         // Fastpath for small numbers
 707         if (result.bitLength() < SMALL_PRIME_THRESHOLD) {
 708 
 709             // Ensure an odd number
 710             if (!result.testBit(0))
 711                 result = result.add(ONE);
 712 
 713             while(true) {
 714                 // Do cheap "pre-test" if applicable
 715                 if (result.bitLength() > 6) {
 716                     long r = result.remainder(SMALL_PRIME_PRODUCT).longValue();
 717                     if ((r%3==0)  || (r%5==0)  || (r%7==0)  || (r%11==0) ||
 718                         (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) ||
 719                         (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0)) {
 720                         result = result.add(TWO);
 721                         continue; // Candidate is composite; try another
 722                     }
 723                 }
 724 
 725                 // All candidates of bitLength 2 and 3 are prime by this point
 726                 if (result.bitLength() < 4)
 727                     return result;
 728 
 729                 // The expensive test
 730                 if (result.primeToCertainty(DEFAULT_PRIME_CERTAINTY, null))
 731                     return result;
 732 
 733                 result = result.add(TWO);
 734             }
 735         }
 736 
 737         // Start at previous even number
 738         if (result.testBit(0))
 739             result = result.subtract(ONE);
 740 
 741         // Looking for the next large prime
 742         int searchLen = (result.bitLength() / 20) * 64;
 743 
 744         while(true) {
 745            BitSieve searchSieve = new BitSieve(result, searchLen);
 746            BigInteger candidate = searchSieve.retrieve(result,
 747                                                  DEFAULT_PRIME_CERTAINTY, null);
 748            if (candidate != null)
 749                return candidate;
 750            result = result.add(BigInteger.valueOf(2 * searchLen));
 751         }
 752     }
 753 
 754     /**
 755      * Returns {@code true} if this BigInteger is probably prime,
 756      * {@code false} if it's definitely composite.
 757      *
 758      * This method assumes bitLength > 2.
 759      *
 760      * @param  certainty a measure of the uncertainty that the caller is
 761      *         willing to tolerate: if the call returns {@code true}
 762      *         the probability that this BigInteger is prime exceeds
 763      *         {@code (1 - 1/2<sup>certainty</sup>)}.  The execution time of
 764      *         this method is proportional to the value of this parameter.
 765      * @return {@code true} if this BigInteger is probably prime,
 766      *         {@code false} if it's definitely composite.
 767      */
 768     boolean primeToCertainty(int certainty, Random random) {
 769         int rounds = 0;
 770         int n = (Math.min(certainty, Integer.MAX_VALUE-1)+1)/2;
 771 
 772         // The relationship between the certainty and the number of rounds
 773         // we perform is given in the draft standard ANSI X9.80, "PRIME
 774         // NUMBER GENERATION, PRIMALITY TESTING, AND PRIMALITY CERTIFICATES".
 775         int sizeInBits = this.bitLength();
 776         if (sizeInBits < 100) {
 777             rounds = 50;
 778             rounds = n < rounds ? n : rounds;
 779             return passesMillerRabin(rounds, random);
 780         }
 781 
 782         if (sizeInBits < 256) {
 783             rounds = 27;
 784         } else if (sizeInBits < 512) {
 785             rounds = 15;
 786         } else if (sizeInBits < 768) {
 787             rounds = 8;
 788         } else if (sizeInBits < 1024) {
 789             rounds = 4;
 790         } else {
 791             rounds = 2;
 792         }
 793         rounds = n < rounds ? n : rounds;
 794 
 795         return passesMillerRabin(rounds, random) && passesLucasLehmer();
 796     }
 797 
 798     /**
 799      * Returns true iff this BigInteger is a Lucas-Lehmer probable prime.
 800      *
 801      * The following assumptions are made:
 802      * This BigInteger is a positive, odd number.
 803      */
 804     private boolean passesLucasLehmer() {
 805         BigInteger thisPlusOne = this.add(ONE);
 806 
 807         // Step 1
 808         int d = 5;
 809         while (jacobiSymbol(d, this) != -1) {
 810             // 5, -7, 9, -11, ...
 811             d = (d<0) ? Math.abs(d)+2 : -(d+2);
 812         }
 813 
 814         // Step 2
 815         BigInteger u = lucasLehmerSequence(d, thisPlusOne, this);
 816 
 817         // Step 3
 818         return u.mod(this).equals(ZERO);
 819     }
 820 
 821     /**
 822      * Computes Jacobi(p,n).
 823      * Assumes n positive, odd, n>=3.
 824      */
 825     private static int jacobiSymbol(int p, BigInteger n) {
 826         if (p == 0)
 827             return 0;
 828 
 829         // Algorithm and comments adapted from Colin Plumb's C library.
 830         int j = 1;
 831         int u = n.mag[n.mag.length-1];
 832 
 833         // Make p positive
 834         if (p < 0) {
 835             p = -p;
 836             int n8 = u & 7;
 837             if ((n8 == 3) || (n8 == 7))
 838                 j = -j; // 3 (011) or 7 (111) mod 8
 839         }
 840 
 841         // Get rid of factors of 2 in p
 842         while ((p & 3) == 0)
 843             p >>= 2;
 844         if ((p & 1) == 0) {
 845             p >>= 1;
 846             if (((u ^ (u>>1)) & 2) != 0)
 847                 j = -j; // 3 (011) or 5 (101) mod 8
 848         }
 849         if (p == 1)
 850             return j;
 851         // Then, apply quadratic reciprocity
 852         if ((p & u & 2) != 0)   // p = u = 3 (mod 4)?
 853             j = -j;
 854         // And reduce u mod p
 855         u = n.mod(BigInteger.valueOf(p)).intValue();
 856 
 857         // Now compute Jacobi(u,p), u < p
 858         while (u != 0) {
 859             while ((u & 3) == 0)
 860                 u >>= 2;
 861             if ((u & 1) == 0) {
 862                 u >>= 1;
 863                 if (((p ^ (p>>1)) & 2) != 0)
 864                     j = -j;     // 3 (011) or 5 (101) mod 8
 865             }
 866             if (u == 1)
 867                 return j;
 868             // Now both u and p are odd, so use quadratic reciprocity
 869             assert (u < p);
 870             int t = u; u = p; p = t;
 871             if ((u & p & 2) != 0) // u = p = 3 (mod 4)?
 872                 j = -j;
 873             // Now u >= p, so it can be reduced
 874             u %= p;
 875         }
 876         return 0;
 877     }
 878 
 879     private static BigInteger lucasLehmerSequence(int z, BigInteger k, BigInteger n) {
 880         BigInteger d = BigInteger.valueOf(z);
 881         BigInteger u = ONE; BigInteger u2;
 882         BigInteger v = ONE; BigInteger v2;
 883 
 884         for (int i=k.bitLength()-2; i>=0; i--) {
 885             u2 = u.multiply(v).mod(n);
 886 
 887             v2 = v.square().add(d.multiply(u.square())).mod(n);
 888             if (v2.testBit(0))
 889                 v2 = v2.subtract(n);
 890 
 891             v2 = v2.shiftRight(1);
 892 
 893             u = u2; v = v2;
 894             if (k.testBit(i)) {
 895                 u2 = u.add(v).mod(n);
 896                 if (u2.testBit(0))
 897                     u2 = u2.subtract(n);
 898 
 899                 u2 = u2.shiftRight(1);
 900                 v2 = v.add(d.multiply(u)).mod(n);
 901                 if (v2.testBit(0))
 902                     v2 = v2.subtract(n);
 903                 v2 = v2.shiftRight(1);
 904 
 905                 u = u2; v = v2;
 906             }
 907         }
 908         return u;
 909     }
 910 
 911     private static volatile Random staticRandom;
 912 
 913     private static Random getSecureRandom() {
 914         if (staticRandom == null) {
 915             staticRandom = new java.security.SecureRandom();
 916         }
 917         return staticRandom;
 918     }
 919 
 920     /**
 921      * Returns true iff this BigInteger passes the specified number of
 922      * Miller-Rabin tests. This test is taken from the DSA spec (NIST FIPS
 923      * 186-2).
 924      *
 925      * The following assumptions are made:
 926      * This BigInteger is a positive, odd number greater than 2.
 927      * iterations<=50.
 928      */
 929     private boolean passesMillerRabin(int iterations, Random rnd) {
 930         // Find a and m such that m is odd and this == 1 + 2**a * m
 931         BigInteger thisMinusOne = this.subtract(ONE);
 932         BigInteger m = thisMinusOne;
 933         int a = m.getLowestSetBit();
 934         m = m.shiftRight(a);
 935 
 936         // Do the tests
 937         if (rnd == null) {
 938             rnd = getSecureRandom();
 939         }
 940         for (int i=0; i<iterations; i++) {
 941             // Generate a uniform random on (1, this)
 942             BigInteger b;
 943             do {
 944                 b = new BigInteger(this.bitLength(), rnd);
 945             } while (b.compareTo(ONE) <= 0 || b.compareTo(this) >= 0);
 946 
 947             int j = 0;
 948             BigInteger z = b.modPow(m, this);
 949             while(!((j==0 && z.equals(ONE)) || z.equals(thisMinusOne))) {
 950                 if (j>0 && z.equals(ONE) || ++j==a)
 951                     return false;
 952                 z = z.modPow(TWO, this);
 953             }
 954         }
 955         return true;
 956     }
 957 
 958     /**
 959      * This internal constructor differs from its public cousin
 960      * with the arguments reversed in two ways: it assumes that its
 961      * arguments are correct, and it doesn't copy the magnitude array.
 962      */
 963     BigInteger(int[] magnitude, int signum) {
 964         this.signum = (magnitude.length==0 ? 0 : signum);
 965         this.mag = magnitude;
 966     }
 967 
 968     /**
 969      * This private constructor is for internal use and assumes that its
 970      * arguments are correct.
 971      */
 972     private BigInteger(byte[] magnitude, int signum) {
 973         this.signum = (magnitude.length==0 ? 0 : signum);
 974         this.mag = stripLeadingZeroBytes(magnitude);
 975     }
 976 
 977     //Static Factory Methods
 978 
 979     /**
 980      * Returns a BigInteger whose value is equal to that of the
 981      * specified {@code long}.  This "static factory method" is
 982      * provided in preference to a ({@code long}) constructor
 983      * because it allows for reuse of frequently used BigIntegers.
 984      *
 985      * @param  val value of the BigInteger to return.
 986      * @return a BigInteger with the specified value.
 987      */
 988     public static BigInteger valueOf(long val) {
 989         // If -MAX_CONSTANT < val < MAX_CONSTANT, return stashed constant
 990         if (val == 0)
 991             return ZERO;
 992         if (val > 0 && val <= MAX_CONSTANT)
 993             return posConst[(int) val];
 994         else if (val < 0 && val >= -MAX_CONSTANT)
 995             return negConst[(int) -val];
 996 
 997         return new BigInteger(val);
 998     }
 999 
1000     /**
1001      * Constructs a BigInteger with the specified value, which may not be zero.
1002      */
1003     private BigInteger(long val) {
1004         if (val < 0) {
1005             val = -val;
1006             signum = -1;
1007         } else {
1008             signum = 1;
1009         }
1010 
1011         int highWord = (int)(val >>> 32);
1012         if (highWord==0) {
1013             mag = new int[1];
1014             mag[0] = (int)val;
1015         } else {
1016             mag = new int[2];
1017             mag[0] = highWord;
1018             mag[1] = (int)val;
1019         }
1020     }
1021 
1022     /**
1023      * Returns a BigInteger with the given two's complement representation.
1024      * Assumes that the input array will not be modified (the returned
1025      * BigInteger will reference the input array if feasible).
1026      */
1027     private static BigInteger valueOf(int val[]) {
1028         return (val[0]>0 ? new BigInteger(val, 1) : new BigInteger(val));
1029     }
1030 
1031     // Constants
1032 
1033     /**
1034      * Initialize static constant array when class is loaded.
1035      */
1036     private final static int MAX_CONSTANT = 16;
1037     private static BigInteger posConst[] = new BigInteger[MAX_CONSTANT+1];
1038     private static BigInteger negConst[] = new BigInteger[MAX_CONSTANT+1];
1039 
1040     /**
1041      * The cache of powers of each radix.  This allows us to not have to
1042      * recalculate powers of radix^(2^n) more than once.  This speeds
1043      * Schoenhage recursive base conversion significantly.
1044      */
1045     private static volatile BigInteger[][] powerCache;
1046 
1047     /** The cache of logarithms of radices for base conversion. */
1048     private static final double[] logCache;
1049 
1050     /** The natural log of 2.  This is used in computing cache indices. */
1051     private static final double LOG_TWO = Math.log(2.0);
1052 
1053     static {
1054         for (int i = 1; i <= MAX_CONSTANT; i++) {
1055             int[] magnitude = new int[1];
1056             magnitude[0] = i;
1057             posConst[i] = new BigInteger(magnitude,  1);
1058             negConst[i] = new BigInteger(magnitude, -1);
1059         }
1060 
1061         /*
1062          * Initialize the cache of radix^(2^x) values used for base conversion
1063          * with just the very first value.  Additional values will be created
1064          * on demand.
1065          */
1066         powerCache = new BigInteger[Character.MAX_RADIX+1][];

1067         logCache = new double[Character.MAX_RADIX+1];
1068 
1069         for (int i=Character.MIN_RADIX; i<=Character.MAX_RADIX; i++)
1070         {
1071             powerCache[i] = new BigInteger[] { BigInteger.valueOf(i) };

1072             logCache[i] = Math.log(i);
1073         }
1074     }
1075 
1076     /**
1077      * The BigInteger constant zero.
1078      *
1079      * @since   1.2
1080      */
1081     public static final BigInteger ZERO = new BigInteger(new int[0], 0);
1082 
1083     /**
1084      * The BigInteger constant one.
1085      *
1086      * @since   1.2
1087      */
1088     public static final BigInteger ONE = valueOf(1);
1089 
1090     /**
1091      * The BigInteger constant two.  (Not exported.)
1092      */
1093     private static final BigInteger TWO = valueOf(2);
1094 
1095     /**
1096      * The BigInteger constant -1.  (Not exported.)
1097      */
1098     private static final BigInteger NEGATIVE_ONE = valueOf(-1);
1099 
1100     /**
1101      * The BigInteger constant ten.
1102      *
1103      * @since   1.5
1104      */
1105     public static final BigInteger TEN = valueOf(10);
1106 
1107     // Arithmetic Operations
1108 
1109     /**
1110      * Returns a BigInteger whose value is {@code (this + val)}.
1111      *
1112      * @param  val value to be added to this BigInteger.
1113      * @return {@code this + val}
1114      */
1115     public BigInteger add(BigInteger val) {
1116         if (val.signum == 0)
1117             return this;
1118         if (signum == 0)
1119             return val;
1120         if (val.signum == signum)
1121             return new BigInteger(add(mag, val.mag), signum);
1122 
1123         int cmp = compareMagnitude(val);
1124         if (cmp == 0)
1125             return ZERO;
1126         int[] resultMag = (cmp > 0 ? subtract(mag, val.mag)
1127                            : subtract(val.mag, mag));
1128         resultMag = trustedStripLeadingZeroInts(resultMag);
1129 
1130         return new BigInteger(resultMag, cmp == signum ? 1 : -1);
1131     }
1132 
1133     /**
1134      * Package private methods used by BigDecimal code to add a BigInteger
1135      * with a long. Assumes val is not equal to INFLATED.
1136      */
1137     BigInteger add(long val) {
1138         if (val == 0)
1139             return this;
1140         if (signum == 0)
1141             return valueOf(val);
1142         if (Long.signum(val) == signum)
1143             return new BigInteger(add(mag, Math.abs(val)), signum);
1144         int cmp = compareMagnitude(val);
1145         if (cmp == 0)
1146             return ZERO;
1147         int[] resultMag = (cmp > 0 ? subtract(mag, Math.abs(val)) : subtract(Math.abs(val), mag));
1148         resultMag = trustedStripLeadingZeroInts(resultMag);
1149         return new BigInteger(resultMag, cmp == signum ? 1 : -1);
1150     }
1151 
1152     /**
1153      * Adds the contents of the int array x and long value val. This
1154      * method allocates a new int array to hold the answer and returns
1155      * a reference to that array.  Assumes x.length &gt; 0 and val is
1156      * non-negative
1157      */
1158     private static int[] add(int[] x, long val) {
1159         int[] y;
1160         long sum = 0;
1161         int xIndex = x.length;
1162         int[] result;
1163         int highWord = (int)(val >>> 32);
1164         if (highWord==0) {
1165             result = new int[xIndex];
1166             sum = (x[--xIndex] & LONG_MASK) + val;
1167             result[xIndex] = (int)sum;
1168         } else {
1169             if (xIndex == 1) {
1170                 result = new int[2];
1171                 sum = val  + (x[0] & LONG_MASK);
1172                 result[1] = (int)sum;
1173                 result[0] = (int)(sum >>> 32);
1174                 return result;
1175             } else {
1176                 result = new int[xIndex];
1177                 sum = (x[--xIndex] & LONG_MASK) + (val & LONG_MASK);
1178                 result[xIndex] = (int)sum;
1179                 sum = (x[--xIndex] & LONG_MASK) + (highWord & LONG_MASK) + (sum >>> 32);
1180                 result[xIndex] = (int)sum;
1181             }
1182         }
1183         // Copy remainder of longer number while carry propagation is required
1184         boolean carry = (sum >>> 32 != 0);
1185         while (xIndex > 0 && carry)
1186             carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
1187         // Copy remainder of longer number
1188         while (xIndex > 0)
1189             result[--xIndex] = x[xIndex];
1190         // Grow result if necessary
1191         if (carry) {
1192             int bigger[] = new int[result.length + 1];
1193             System.arraycopy(result, 0, bigger, 1, result.length);
1194             bigger[0] = 0x01;
1195             return bigger;
1196         }
1197         return result;
1198     }
1199 
1200     /**
1201      * Adds the contents of the int arrays x and y. This method allocates
1202      * a new int array to hold the answer and returns a reference to that
1203      * array.
1204      */
1205     private static int[] add(int[] x, int[] y) {
1206         // If x is shorter, swap the two arrays
1207         if (x.length < y.length) {
1208             int[] tmp = x;
1209             x = y;
1210             y = tmp;
1211         }
1212 
1213         int xIndex = x.length;
1214         int yIndex = y.length;
1215         int result[] = new int[xIndex];
1216         long sum = 0;
1217         if(yIndex==1) {
1218             sum = (x[--xIndex] & LONG_MASK) + (y[0] & LONG_MASK) ;
1219             result[xIndex] = (int)sum;
1220         } else {
1221             // Add common parts of both numbers
1222             while(yIndex > 0) {
1223                 sum = (x[--xIndex] & LONG_MASK) +
1224                       (y[--yIndex] & LONG_MASK) + (sum >>> 32);
1225                 result[xIndex] = (int)sum;
1226             }
1227         }
1228         // Copy remainder of longer number while carry propagation is required
1229         boolean carry = (sum >>> 32 != 0);
1230         while (xIndex > 0 && carry)
1231             carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
1232 
1233         // Copy remainder of longer number
1234         while (xIndex > 0)
1235             result[--xIndex] = x[xIndex];
1236 
1237         // Grow result if necessary
1238         if (carry) {
1239             int bigger[] = new int[result.length + 1];
1240             System.arraycopy(result, 0, bigger, 1, result.length);
1241             bigger[0] = 0x01;
1242             return bigger;
1243         }
1244         return result;
1245     }
1246 
1247     private static int[] subtract(long val, int[] little) {
1248         int highWord = (int)(val >>> 32);
1249         if (highWord==0) {
1250             int result[] = new int[1];
1251             result[0] = (int)(val - (little[0] & LONG_MASK));
1252             return result;
1253         } else {
1254             int result[] = new int[2];
1255             if(little.length==1) {
1256                 long difference = ((int)val & LONG_MASK) - (little[0] & LONG_MASK);
1257                 result[1] = (int)difference;
1258                 // Subtract remainder of longer number while borrow propagates
1259                 boolean borrow = (difference >> 32 != 0);
1260                 if(borrow) {
1261                     result[0] = highWord - 1;
1262                 } else {        // Copy remainder of longer number
1263                     result[0] = highWord;
1264                 }
1265                 return result;
1266             } else { // little.length==2
1267                 long difference = ((int)val & LONG_MASK) - (little[1] & LONG_MASK);
1268                 result[1] = (int)difference;
1269                 difference = (highWord & LONG_MASK) - (little[0] & LONG_MASK) + (difference >> 32);
1270                 result[0] = (int)difference;
1271                 return result;
1272             }
1273         }
1274     }
1275 
1276     /**
1277      * Subtracts the contents of the second argument (val) from the
1278      * first (big).  The first int array (big) must represent a larger number
1279      * than the second.  This method allocates the space necessary to hold the
1280      * answer.
1281      * assumes val &gt;= 0
1282      */
1283     private static int[] subtract(int[] big, long val) {
1284         int highWord = (int)(val >>> 32);
1285         int bigIndex = big.length;
1286         int result[] = new int[bigIndex];
1287         long difference = 0;
1288 
1289         if (highWord==0) {
1290             difference = (big[--bigIndex] & LONG_MASK) - val;
1291             result[bigIndex] = (int)difference;
1292         } else {
1293             difference = (big[--bigIndex] & LONG_MASK) - (val & LONG_MASK);
1294             result[bigIndex] = (int)difference;
1295             difference = (big[--bigIndex] & LONG_MASK) - (highWord & LONG_MASK) + (difference >> 32);
1296             result[bigIndex] = (int)difference;
1297         }
1298 
1299 
1300         // Subtract remainder of longer number while borrow propagates
1301         boolean borrow = (difference >> 32 != 0);
1302         while (bigIndex > 0 && borrow)
1303             borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
1304 
1305         // Copy remainder of longer number
1306         while (bigIndex > 0)
1307             result[--bigIndex] = big[bigIndex];
1308 
1309         return result;
1310     }
1311 
1312     /**
1313      * Returns a BigInteger whose value is {@code (this - val)}.
1314      *
1315      * @param  val value to be subtracted from this BigInteger.
1316      * @return {@code this - val}
1317      */
1318     public BigInteger subtract(BigInteger val) {
1319         if (val.signum == 0)
1320             return this;
1321         if (signum == 0)
1322             return val.negate();
1323         if (val.signum != signum)
1324             return new BigInteger(add(mag, val.mag), signum);
1325 
1326         int cmp = compareMagnitude(val);
1327         if (cmp == 0)
1328             return ZERO;
1329         int[] resultMag = (cmp > 0 ? subtract(mag, val.mag)
1330                            : subtract(val.mag, mag));
1331         resultMag = trustedStripLeadingZeroInts(resultMag);
1332         return new BigInteger(resultMag, cmp == signum ? 1 : -1);
1333     }
1334 
1335     /**
1336      * Subtracts the contents of the second int arrays (little) from the
1337      * first (big).  The first int array (big) must represent a larger number
1338      * than the second.  This method allocates the space necessary to hold the
1339      * answer.
1340      */
1341     private static int[] subtract(int[] big, int[] little) {
1342         int bigIndex = big.length;
1343         int result[] = new int[bigIndex];
1344         int littleIndex = little.length;
1345         long difference = 0;
1346 
1347         // Subtract common parts of both numbers
1348         while(littleIndex > 0) {
1349             difference = (big[--bigIndex] & LONG_MASK) -
1350                          (little[--littleIndex] & LONG_MASK) +
1351                          (difference >> 32);
1352             result[bigIndex] = (int)difference;
1353         }
1354 
1355         // Subtract remainder of longer number while borrow propagates
1356         boolean borrow = (difference >> 32 != 0);
1357         while (bigIndex > 0 && borrow)
1358             borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
1359 
1360         // Copy remainder of longer number
1361         while (bigIndex > 0)
1362             result[--bigIndex] = big[bigIndex];
1363 
1364         return result;
1365     }
1366 
1367     /**
1368      * Returns a BigInteger whose value is {@code (this * val)}.
1369      *
1370      * @param  val value to be multiplied by this BigInteger.
1371      * @return {@code this * val}
1372      */
1373     public BigInteger multiply(BigInteger val) {
1374         if (val.signum == 0 || signum == 0)
1375             return ZERO;
1376 
1377         int xlen = mag.length;
1378         int ylen = val.mag.length;
1379 
1380         if ((xlen < KARATSUBA_THRESHOLD) || (ylen < KARATSUBA_THRESHOLD))
1381         {
1382             int resultSign = signum == val.signum ? 1 : -1;
1383             if (val.mag.length == 1) {
1384                 return multiplyByInt(mag,val.mag[0], resultSign);
1385             }
1386             if(mag.length == 1) {
1387                 return multiplyByInt(val.mag,mag[0], resultSign);
1388             }
1389             int[] result = multiplyToLen(mag, xlen,
1390                                          val.mag, ylen, null);
1391             result = trustedStripLeadingZeroInts(result);
1392             return new BigInteger(result, resultSign);
1393         }
1394         else
1395             if ((xlen < TOOM_COOK_THRESHOLD) && (ylen < TOOM_COOK_THRESHOLD))
1396                 return multiplyKaratsuba(this, val);
1397             else
1398                 return multiplyToomCook3(this, val);
1399     }
1400 
1401     private static BigInteger multiplyByInt(int[] x, int y, int sign) {
1402         if(Integer.bitCount(y)==1) {
1403             return new BigInteger(shiftLeft(x,Integer.numberOfTrailingZeros(y)), sign);
1404         }
1405         int xlen = x.length;
1406         int[] rmag =  new int[xlen + 1];
1407         long carry = 0;
1408         long yl = y & LONG_MASK;
1409         int rstart = rmag.length - 1;
1410         for (int i = xlen - 1; i >= 0; i--) {
1411             long product = (x[i] & LONG_MASK) * yl + carry;
1412             rmag[rstart--] = (int)product;
1413             carry = product >>> 32;
1414         }
1415         if (carry == 0L) {
1416             rmag = java.util.Arrays.copyOfRange(rmag, 1, rmag.length);
1417         } else {
1418             rmag[rstart] = (int)carry;
1419         }
1420         return new BigInteger(rmag, sign);
1421     }
1422 
1423     /**
1424      * Package private methods used by BigDecimal code to multiply a BigInteger
1425      * with a long. Assumes v is not equal to INFLATED.
1426      */
1427     BigInteger multiply(long v) {
1428         if (v == 0 || signum == 0)
1429           return ZERO;
1430         if (v == BigDecimal.INFLATED)
1431             return multiply(BigInteger.valueOf(v));
1432         int rsign = (v > 0 ? signum : -signum);
1433         if (v < 0)
1434             v = -v;
1435         long dh = v >>> 32;      // higher order bits
1436         long dl = v & LONG_MASK; // lower order bits
1437 
1438         int xlen = mag.length;
1439         int[] value = mag;
1440         int[] rmag = (dh == 0L) ? (new int[xlen + 1]) : (new int[xlen + 2]);
1441         long carry = 0;
1442         int rstart = rmag.length - 1;
1443         for (int i = xlen - 1; i >= 0; i--) {
1444             long product = (value[i] & LONG_MASK) * dl + carry;
1445             rmag[rstart--] = (int)product;
1446             carry = product >>> 32;
1447         }
1448         rmag[rstart] = (int)carry;
1449         if (dh != 0L) {
1450             carry = 0;
1451             rstart = rmag.length - 2;
1452             for (int i = xlen - 1; i >= 0; i--) {
1453                 long product = (value[i] & LONG_MASK) * dh +
1454                     (rmag[rstart] & LONG_MASK) + carry;
1455                 rmag[rstart--] = (int)product;
1456                 carry = product >>> 32;
1457             }
1458             rmag[0] = (int)carry;
1459         }
1460         if (carry == 0L)
1461             rmag = java.util.Arrays.copyOfRange(rmag, 1, rmag.length);
1462         return new BigInteger(rmag, rsign);
1463     }
1464 
1465     /**
1466      * Multiplies int arrays x and y to the specified lengths and places
1467      * the result into z. There will be no leading zeros in the resultant array.
1468      */
1469     private int[] multiplyToLen(int[] x, int xlen, int[] y, int ylen, int[] z) {
1470         int xstart = xlen - 1;
1471         int ystart = ylen - 1;
1472 
1473         if (z == null || z.length < (xlen+ ylen))
1474             z = new int[xlen+ylen];
1475 
1476         long carry = 0;
1477         for (int j=ystart, k=ystart+1+xstart; j>=0; j--, k--) {
1478             long product = (y[j] & LONG_MASK) *
1479                            (x[xstart] & LONG_MASK) + carry;
1480             z[k] = (int)product;
1481             carry = product >>> 32;
1482         }
1483         z[xstart] = (int)carry;
1484 
1485         for (int i = xstart-1; i >= 0; i--) {
1486             carry = 0;
1487             for (int j=ystart, k=ystart+1+i; j>=0; j--, k--) {
1488                 long product = (y[j] & LONG_MASK) *
1489                                (x[i] & LONG_MASK) +
1490                                (z[k] & LONG_MASK) + carry;
1491                 z[k] = (int)product;
1492                 carry = product >>> 32;
1493             }
1494             z[i] = (int)carry;
1495         }
1496         return z;
1497     }
1498 
1499     /**
1500      * Multiplies two BigIntegers using the Karatsuba multiplication
1501      * algorithm.  This is a recursive divide-and-conquer algorithm which is
1502      * more efficient for large numbers than what is commonly called the
1503      * "grade-school" algorithm used in multiplyToLen.  If the numbers to be
1504      * multiplied have length n, the "grade-school" algorithm has an
1505      * asymptotic complexity of O(n^2).  In contrast, the Karatsuba algorithm
1506      * has complexity of O(n^(log2(3))), or O(n^1.585).  It achieves this
1507      * increased performance by doing 3 multiplies instead of 4 when
1508      * evaluating the product.  As it has some overhead, should be used when
1509      * both numbers are larger than a certain threshold (found
1510      * experimentally).
1511      *
1512      * See:  http://en.wikipedia.org/wiki/Karatsuba_algorithm
1513      */
1514     private static BigInteger multiplyKaratsuba(BigInteger x, BigInteger y)
1515     {
1516         int xlen = x.mag.length;
1517         int ylen = y.mag.length;
1518 
1519         // The number of ints in each half of the number.
1520         int half = (Math.max(xlen, ylen)+1) / 2;
1521 
1522         // xl and yl are the lower halves of x and y respectively,
1523         // xh and yh are the upper halves.
1524         BigInteger xl = x.getLower(half);
1525         BigInteger xh = x.getUpper(half);
1526         BigInteger yl = y.getLower(half);
1527         BigInteger yh = y.getUpper(half);
1528 
1529         BigInteger p1 = xh.multiply(yh);  // p1 = xh*yh
1530         BigInteger p2 = xl.multiply(yl);  // p2 = xl*yl
1531 
1532         // p3=(xh+xl)*(yh+yl)
1533         BigInteger p3 = xh.add(xl).multiply(yh.add(yl));
1534 
1535         // result = p1 * 2^(32*2*half) + (p3 - p1 - p2) * 2^(32*half) + p2
1536         BigInteger result = p1.shiftLeft(32*half).add(p3.subtract(p1).subtract(p2)).shiftLeft(32*half).add(p2);
1537 
1538         if (x.signum != y.signum)
1539             return result.negate();
1540         else
1541             return result;
1542     }
1543 
1544     /**
1545      * Multiplies two BigIntegers using a 3-way Toom-Cook multiplication
1546      * algorithm.  This is a recursive divide-and-conquer algorithm which is
1547      * more efficient for large numbers than what is commonly called the
1548      * "grade-school" algorithm used in multiplyToLen.  If the numbers to be
1549      * multiplied have length n, the "grade-school" algorithm has an
1550      * asymptotic complexity of O(n^2).  In contrast, 3-way Toom-Cook has a
1551      * complexity of about O(n^1.465).  It achieves this increased asymptotic
1552      * performance by breaking each number into three parts and by doing 5
1553      * multiplies instead of 9 when evaluating the product.  Due to overhead
1554      * (additions, shifts, and one division) in the Toom-Cook algorithm, it
1555      * should only be used when both numbers are larger than a certain
1556      * threshold (found experimentally).  This threshold is generally larger
1557      * than that for Karatsuba multiplication, so this algorithm is generally
1558      * only used when numbers become significantly larger.
1559      *
1560      * The algorithm used is the "optimal" 3-way Toom-Cook algorithm outlined
1561      * by Marco Bodrato.
1562      *
1563      *  See: http://bodrato.it/toom-cook/
1564      *       http://bodrato.it/papers/#WAIFI2007
1565      *
1566      * "Towards Optimal Toom-Cook Multiplication for Univariate and
1567      * Multivariate Polynomials in Characteristic 2 and 0." by Marco BODRATO;
1568      * In C.Carlet and B.Sunar, Eds., "WAIFI'07 proceedings", p. 116-133,
1569      * LNCS #4547. Springer, Madrid, Spain, June 21-22, 2007.
1570      *
1571      */
1572     private static BigInteger multiplyToomCook3(BigInteger a, BigInteger b)
1573     {
1574         int alen = a.mag.length;
1575         int blen = b.mag.length;
1576 
1577         int largest = Math.max(alen, blen);
1578 
1579         // k is the size (in ints) of the lower-order slices.
1580         int k = (largest+2)/3;   // Equal to ceil(largest/3)
1581 
1582         // r is the size (in ints) of the highest-order slice.
1583         int r = largest - 2*k;
1584 
1585         // Obtain slices of the numbers. a2 and b2 are the most significant
1586         // bits of the numbers a and b, and a0 and b0 the least significant.
1587         BigInteger a0, a1, a2, b0, b1, b2;
1588         a2 = a.getToomSlice(k, r, 0, largest);
1589         a1 = a.getToomSlice(k, r, 1, largest);
1590         a0 = a.getToomSlice(k, r, 2, largest);
1591         b2 = b.getToomSlice(k, r, 0, largest);
1592         b1 = b.getToomSlice(k, r, 1, largest);
1593         b0 = b.getToomSlice(k, r, 2, largest);
1594 
1595         BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1, db1;
1596 
1597         v0 = a0.multiply(b0);
1598         da1 = a2.add(a0);
1599         db1 = b2.add(b0);
1600         vm1 = da1.subtract(a1).multiply(db1.subtract(b1));
1601         da1 = da1.add(a1);
1602         db1 = db1.add(b1);
1603         v1 = da1.multiply(db1);
1604         v2 = da1.add(a2).shiftLeft(1).subtract(a0).multiply(
1605              db1.add(b2).shiftLeft(1).subtract(b0));
1606         vinf = a2.multiply(b2);
1607 
1608         /* The algorithm requires two divisions by 2 and one by 3.
1609            All divisions are known to be exact, that is, they do not produce
1610            remainders, and all results are positive.  The divisions by 2 are
1611            implemented as right shifts which are relatively efficient, leaving
1612            only an exact division by 3, which is done by a specialized
1613            linear-time algorithm. */
1614         t2 = v2.subtract(vm1).exactDivideBy3();
1615         tm1 = v1.subtract(vm1).shiftRight(1);
1616         t1 = v1.subtract(v0);
1617         t2 = t2.subtract(t1).shiftRight(1);
1618         t1 = t1.subtract(tm1).subtract(vinf);
1619         t2 = t2.subtract(vinf.shiftLeft(1));
1620         tm1 = tm1.subtract(t2);
1621 
1622         // Number of bits to shift left.
1623         int ss = k*32;
1624 
1625         BigInteger result = vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
1626 
1627         if (a.signum != b.signum)
1628             return result.negate();
1629         else
1630             return result;
1631     }
1632 
1633 
1634     /**
1635      * Returns a slice of a BigInteger for use in Toom-Cook multiplication.
1636      *
1637      * @param lowerSize The size of the lower-order bit slices.
1638      * @param upperSize The size of the higher-order bit slices.
1639      * @param slice The index of which slice is requested, which must be a
1640      * number from 0 to size-1. Slice 0 is the highest-order bits, and slice
1641      * size-1 are the lowest-order bits. Slice 0 may be of different size than
1642      * the other slices.
1643      * @param fullsize The size of the larger integer array, used to align
1644      * slices to the appropriate position when multiplying different-sized
1645      * numbers.
1646      */
1647     private BigInteger getToomSlice(int lowerSize, int upperSize, int slice,
1648                                     int fullsize)
1649     {
1650         int start, end, sliceSize, len, offset;
1651 
1652         len = mag.length;
1653         offset = fullsize - len;
1654 
1655         if (slice == 0)
1656         {
1657             start = 0 - offset;
1658             end = upperSize - 1 - offset;
1659         }
1660         else
1661         {
1662             start = upperSize + (slice-1)*lowerSize - offset;
1663             end = start + lowerSize - 1;
1664         }
1665 
1666         if (start < 0)
1667             start = 0;
1668         if (end < 0)
1669            return ZERO;
1670 
1671         sliceSize = (end-start) + 1;
1672 
1673         if (sliceSize <= 0)
1674             return ZERO;
1675 
1676         // While performing Toom-Cook, all slices are positive and
1677         // the sign is adjusted when the final number is composed.
1678         if (start==0 && sliceSize >= len)
1679             return this.abs();
1680 
1681         int intSlice[] = new int[sliceSize];
1682         System.arraycopy(mag, start, intSlice, 0, sliceSize);
1683 
1684         return new BigInteger(trustedStripLeadingZeroInts(intSlice), 1);
1685     }
1686 
1687     /**
1688      * Does an exact division (that is, the remainder is known to be zero)
1689      * of the specified number by 3.  This is used in Toom-Cook
1690      * multiplication.  This is an efficient algorithm that runs in linear
1691      * time.  If the argument is not exactly divisible by 3, results are
1692      * undefined.  Note that this is expected to be called with positive
1693      * arguments only.
1694      */
1695     private BigInteger exactDivideBy3()
1696     {
1697         int len = mag.length;
1698         int[] result = new int[len];
1699         long x, w, q, borrow;
1700         borrow = 0L;
1701         for (int i=len-1; i>=0; i--)
1702         {
1703             x = (mag[i] & LONG_MASK);
1704             w = x - borrow;
1705             if (borrow > x)       // Did we make the number go negative?
1706                 borrow = 1L;
1707             else
1708                 borrow = 0L;
1709 
1710             // 0xAAAAAAAB is the modular inverse of 3 (mod 2^32).  Thus,
1711             // the effect of this is to divide by 3 (mod 2^32).
1712             // This is much faster than division on most architectures.
1713             q = (w * 0xAAAAAAABL) & LONG_MASK;
1714             result[i] = (int) q;
1715 
1716             // Now check the borrow. The second check can of course be
1717             // eliminated if the first fails.
1718             if (q >= 0x55555556L)
1719             {
1720                 borrow++;
1721                 if (q >= 0xAAAAAAABL)
1722                     borrow++;
1723             }
1724         }
1725         result = trustedStripLeadingZeroInts(result);
1726         return new BigInteger(result, signum);
1727     }
1728 
1729     /**
1730      * Returns a new BigInteger representing n lower ints of the number.
1731      * This is used by Karatsuba multiplication and Karatsuba squaring.
1732      */
1733     private BigInteger getLower(int n) {
1734         int len = mag.length;
1735 
1736         if (len <= n)
1737             return this;
1738 
1739         int lowerInts[] = new int[n];
1740         System.arraycopy(mag, len-n, lowerInts, 0, n);
1741 
1742         return new BigInteger(trustedStripLeadingZeroInts(lowerInts), 1);
1743     }
1744 
1745     /**
1746      * Returns a new BigInteger representing mag.length-n upper
1747      * ints of the number.  This is used by Karatsuba multiplication and
1748      * Karatsuba squaring.
1749      */
1750     private BigInteger getUpper(int n) {
1751         int len = mag.length;
1752 
1753         if (len <= n)
1754             return ZERO;
1755 
1756         int upperLen = len - n;
1757         int upperInts[] = new int[upperLen];
1758         System.arraycopy(mag, 0, upperInts, 0, upperLen);
1759 
1760         return new BigInteger(trustedStripLeadingZeroInts(upperInts), 1);
1761     }
1762 
1763     // Squaring
1764 
1765     /**
1766      * Returns a BigInteger whose value is {@code (this<sup>2</sup>)}.
1767      *
1768      * @return {@code this<sup>2</sup>}
1769      */
1770     private BigInteger square() {
1771         if (signum == 0)
1772             return ZERO;
1773         int len = mag.length;
1774 
1775         if (len < KARATSUBA_SQUARE_THRESHOLD)
1776         {
1777             int[] z = squareToLen(mag, len, null);
1778             return new BigInteger(trustedStripLeadingZeroInts(z), 1);
1779         }
1780         else
1781             if (len < TOOM_COOK_SQUARE_THRESHOLD)
1782                 return squareKaratsuba();
1783             else
1784                return squareToomCook3();
1785     }
1786 
1787     /**
1788      * Squares the contents of the int array x. The result is placed into the
1789      * int array z.  The contents of x are not changed.
1790      */
1791     private static final int[] squareToLen(int[] x, int len, int[] z) {
1792         /*
1793          * The algorithm used here is adapted from Colin Plumb's C library.
1794          * Technique: Consider the partial products in the multiplication
1795          * of "abcde" by itself:
1796          *
1797          *               a  b  c  d  e
1798          *            *  a  b  c  d  e
1799          *          ==================
1800          *              ae be ce de ee
1801          *           ad bd cd dd de
1802          *        ac bc cc cd ce
1803          *     ab bb bc bd be
1804          *  aa ab ac ad ae
1805          *
1806          * Note that everything above the main diagonal:
1807          *              ae be ce de = (abcd) * e
1808          *           ad bd cd       = (abc) * d
1809          *        ac bc             = (ab) * c
1810          *     ab                   = (a) * b
1811          *
1812          * is a copy of everything below the main diagonal:
1813          *                       de
1814          *                 cd ce
1815          *           bc bd be
1816          *     ab ac ad ae
1817          *
1818          * Thus, the sum is 2 * (off the diagonal) + diagonal.
1819          *
1820          * This is accumulated beginning with the diagonal (which
1821          * consist of the squares of the digits of the input), which is then
1822          * divided by two, the off-diagonal added, and multiplied by two
1823          * again.  The low bit is simply a copy of the low bit of the
1824          * input, so it doesn't need special care.
1825          */
1826         int zlen = len << 1;
1827         if (z == null || z.length < zlen)
1828             z = new int[zlen];
1829 
1830         // Store the squares, right shifted one bit (i.e., divided by 2)
1831         int lastProductLowWord = 0;
1832         for (int j=0, i=0; j<len; j++) {
1833             long piece = (x[j] & LONG_MASK);
1834             long product = piece * piece;
1835             z[i++] = (lastProductLowWord << 31) | (int)(product >>> 33);
1836             z[i++] = (int)(product >>> 1);
1837             lastProductLowWord = (int)product;
1838         }
1839 
1840         // Add in off-diagonal sums
1841         for (int i=len, offset=1; i>0; i--, offset+=2) {
1842             int t = x[i-1];
1843             t = mulAdd(z, x, offset, i-1, t);
1844             addOne(z, offset-1, i, t);
1845         }
1846 
1847         // Shift back up and set low bit
1848         primitiveLeftShift(z, zlen, 1);
1849         z[zlen-1] |= x[len-1] & 1;
1850 
1851         return z;
1852     }
1853 
1854     /**
1855      * Squares a BigInteger using the Karatsuba squaring algorithm.  It should
1856      * be used when both numbers are larger than a certain threshold (found
1857      * experimentally).  It is a recursive divide-and-conquer algorithm that
1858      * has better asymptotic performance than the algorithm used in
1859      * squareToLen.
1860      */
1861     private BigInteger squareKaratsuba()
1862     {
1863         int half = (mag.length+1) / 2;
1864 
1865         BigInteger xl = getLower(half);
1866         BigInteger xh = getUpper(half);
1867 
1868         BigInteger xhs = xh.square();  // xhs = xh^2
1869         BigInteger xls = xl.square();  // xls = xl^2
1870 
1871         // xh^2 << 64  +  (((xl+xh)^2 - (xh^2 + xl^2)) << 32) + xl^2
1872         return xhs.shiftLeft(half*32).add(xl.add(xh).square().subtract(xhs.add(xls))).shiftLeft(half*32).add(xls);
1873     }
1874 
1875     /**
1876      * Squares a BigInteger using the 3-way Toom-Cook squaring algorithm.  It
1877      * should be used when both numbers are larger than a certain threshold
1878      * (found experimentally).  It is a recursive divide-and-conquer algorithm
1879      * that has better asymptotic performance than the algorithm used in
1880      * squareToLen or squareKaratsuba.
1881      */
1882     private BigInteger squareToomCook3()
1883     {
1884         int len = mag.length;
1885 
1886         // k is the size (in ints) of the lower-order slices.
1887         int k = (len+2)/3;   // Equal to ceil(largest/3)
1888 
1889         // r is the size (in ints) of the highest-order slice.
1890         int r = len - 2*k;
1891 
1892         // Obtain slices of the numbers. a2 is the most significant
1893         // bits of the number, and a0 the least significant.
1894         BigInteger a0, a1, a2;
1895         a2 = getToomSlice(k, r, 0, len);
1896         a1 = getToomSlice(k, r, 1, len);
1897         a0 = getToomSlice(k, r, 2, len);
1898         BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1;
1899 
1900         v0 = a0.square();
1901         da1 = a2.add(a0);
1902         vm1 = da1.subtract(a1).square();
1903         da1 = da1.add(a1);
1904         v1 = da1.square();
1905         vinf = a2.square();
1906         v2 = da1.add(a2).shiftLeft(1).subtract(a0).square();
1907 
1908         /* The algorithm requires two divisions by 2 and one by 3.
1909            All divisions are known to be exact, that is, they do not produce
1910            remainders, and all results are positive.  The divisions by 2 are
1911            implemented as right shifts which are relatively efficient, leaving
1912            only a division by 3.
1913            The division by 3 is done by an optimized algorithm for this case.
1914         */
1915         t2 = v2.subtract(vm1).exactDivideBy3();
1916         tm1 = v1.subtract(vm1).shiftRight(1);
1917         t1 = v1.subtract(v0);
1918         t2 = t2.subtract(t1).shiftRight(1);
1919         t1 = t1.subtract(tm1).subtract(vinf);
1920         t2 = t2.subtract(vinf.shiftLeft(1));
1921         tm1 = tm1.subtract(t2);
1922 
1923         // Number of bits to shift left.
1924         int ss = k*32;
1925 
1926         return vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
1927     }
1928 
1929     // Division
1930 
1931     /**
1932      * Returns a BigInteger whose value is {@code (this / val)}.
1933      *
1934      * @param  val value by which this BigInteger is to be divided.
1935      * @return {@code this / val}
1936      * @throws ArithmeticException if {@code val} is zero.
1937      */
1938     public BigInteger divide(BigInteger val) {
1939         MutableBigInteger q = new MutableBigInteger(),
1940                           a = new MutableBigInteger(this.mag),
1941                           b = new MutableBigInteger(val.mag);
1942 
1943         a.divide(b, q, false);
1944         return q.toBigInteger(this.signum * val.signum);
1945     }
1946 
1947     /**
1948      * Returns an array of two BigIntegers containing {@code (this / val)}
1949      * followed by {@code (this % val)}.
1950      *
1951      * @param  val value by which this BigInteger is to be divided, and the
1952      *         remainder computed.
1953      * @return an array of two BigIntegers: the quotient {@code (this / val)}
1954      *         is the initial element, and the remainder {@code (this % val)}
1955      *         is the final element.
1956      * @throws ArithmeticException if {@code val} is zero.
1957      */
1958     public BigInteger[] divideAndRemainder(BigInteger val) {
1959         BigInteger[] result = new BigInteger[2];
1960         MutableBigInteger q = new MutableBigInteger(),
1961                           a = new MutableBigInteger(this.mag),
1962                           b = new MutableBigInteger(val.mag);
1963         MutableBigInteger r = a.divide(b, q);
1964         result[0] = q.toBigInteger(this.signum == val.signum ? 1 : -1);
1965         result[1] = r.toBigInteger(this.signum);
1966         return result;
1967     }
1968 
1969     /**
1970      * Returns a BigInteger whose value is {@code (this % val)}.
1971      *
1972      * @param  val value by which this BigInteger is to be divided, and the
1973      *         remainder computed.
1974      * @return {@code this % val}
1975      * @throws ArithmeticException if {@code val} is zero.
1976      */
1977     public BigInteger remainder(BigInteger val) {
1978         MutableBigInteger q = new MutableBigInteger(),
1979                           a = new MutableBigInteger(this.mag),
1980                           b = new MutableBigInteger(val.mag);
1981 
1982         return a.divide(b, q).toBigInteger(this.signum);
1983     }
1984 
1985     /**
1986      * Returns a BigInteger whose value is <tt>(this<sup>exponent</sup>)</tt>.
1987      * Note that {@code exponent} is an integer rather than a BigInteger.
1988      *
1989      * @param  exponent exponent to which this BigInteger is to be raised.
1990      * @return <tt>this<sup>exponent</sup></tt>
1991      * @throws ArithmeticException {@code exponent} is negative.  (This would
1992      *         cause the operation to yield a non-integer value.)
1993      */
1994     public BigInteger pow(int exponent) {
1995         if (exponent < 0)
1996             throw new ArithmeticException("Negative exponent");
1997         if (signum==0)
1998             return (exponent==0 ? ONE : this);
1999 
2000         BigInteger partToSquare = this.abs();
2001 
2002         // Factor out powers of two from the base, as the exponentiation of
2003         // these can be done by left shifts only.
2004         // The remaining part can then be exponentiated faster.  The
2005         // powers of two will be multiplied back at the end.
2006         int powersOfTwo = partToSquare.getLowestSetBit();
2007 
2008         int remainingBits;
2009 
2010         // Factor the powers of two out quickly by shifting right, if needed.
2011         if (powersOfTwo > 0)
2012         {
2013             partToSquare = partToSquare.shiftRight(powersOfTwo);
2014             remainingBits = partToSquare.bitLength();
2015             if (remainingBits == 1)  // Nothing left but +/- 1?
2016                 if (signum<0 && (exponent&1)==1)
2017                     return NEGATIVE_ONE.shiftLeft(powersOfTwo*exponent);
2018                 else
2019                     return ONE.shiftLeft(powersOfTwo*exponent);
2020         }
2021         else
2022         {
2023             remainingBits = partToSquare.bitLength();
2024             if (remainingBits == 1)  // Nothing left but +/- 1?
2025                 if (signum<0 && (exponent&1)==1)
2026                     return NEGATIVE_ONE;
2027                 else
2028                     return ONE;
2029         }
2030 
2031         // This is a quick way to approximate the size of the result,
2032         // similar to doing log2[n] * exponent.  This will give an upper bound
2033         // of how big the result can be, and which algorithm to use.
2034         int scaleFactor = remainingBits * exponent;
2035 
2036         // Use slightly different algorithms for small and large operands.
2037         // See if the result will safely fit into a long. (Largest 2^63-1)
2038         if (partToSquare.mag.length==1 && scaleFactor <= 62)
2039         {
2040             // Small number algorithm.  Everything fits into a long.
2041             int newSign = (signum<0 && (exponent&1)==1 ? -1 : 1);
2042             long result = 1;
2043             long baseToPow2 = partToSquare.mag[0] & LONG_MASK;
2044 
2045             int workingExponent = exponent;
2046 
2047             // Perform exponentiation using repeated squaring trick
2048             while (workingExponent != 0) {
2049                 if ((workingExponent & 1)==1)
2050                     result = result * baseToPow2;
2051 
2052                 if ((workingExponent >>>= 1) != 0)
2053                     baseToPow2 = baseToPow2 * baseToPow2;
2054             }
2055 
2056             // Multiply back the powers of two (quickly, by shifting left)
2057             if (powersOfTwo > 0)
2058             {
2059                 int bitsToShift = powersOfTwo*exponent;
2060                 if (bitsToShift + scaleFactor <= 62) // Fits in long?
2061                     return valueOf((result << bitsToShift) * newSign);
2062                 else
2063                     return valueOf(result*newSign).shiftLeft(bitsToShift);
2064             }
2065             else
2066                 return valueOf(result*newSign);
2067         }
2068         else
2069         {
2070             // Large number algorithm.  This is basically identical to
2071             // the algorithm above, but calls multiply() and square()
2072             // which may use more efficient algorithms for large numbers.
2073             BigInteger answer = ONE;
2074 
2075             int workingExponent = exponent;
2076             // Perform exponentiation using repeated squaring trick
2077             while (workingExponent != 0) {
2078                 if ((workingExponent & 1)==1)
2079                     answer = answer.multiply(partToSquare);
2080 
2081                 if ((workingExponent >>>= 1) != 0)
2082                     partToSquare = partToSquare.square();
2083             }
2084             // Multiply back the (exponentiated) powers of two (quickly,
2085             // by shifting left)
2086             if (powersOfTwo > 0)
2087                 answer = answer.shiftLeft(powersOfTwo*exponent);
2088 
2089             if (signum<0 && (exponent&1)==1)
2090                 return answer.negate();
2091             else
2092                 return answer;
2093         }
2094     }
2095 
2096     /**
2097      * Returns a BigInteger whose value is the greatest common divisor of
2098      * {@code abs(this)} and {@code abs(val)}.  Returns 0 if
2099      * {@code this==0 && val==0}.
2100      *
2101      * @param  val value with which the GCD is to be computed.
2102      * @return {@code GCD(abs(this), abs(val))}
2103      */
2104     public BigInteger gcd(BigInteger val) {
2105         if (val.signum == 0)
2106             return this.abs();
2107         else if (this.signum == 0)
2108             return val.abs();
2109 
2110         MutableBigInteger a = new MutableBigInteger(this);
2111         MutableBigInteger b = new MutableBigInteger(val);
2112 
2113         MutableBigInteger result = a.hybridGCD(b);
2114 
2115         return result.toBigInteger(1);
2116     }
2117 
2118     /**
2119      * Package private method to return bit length for an integer.
2120      */
2121     static int bitLengthForInt(int n) {
2122         return 32 - Integer.numberOfLeadingZeros(n);
2123     }
2124 
2125     /**
2126      * Left shift int array a up to len by n bits. Returns the array that
2127      * results from the shift since space may have to be reallocated.
2128      */
2129     private static int[] leftShift(int[] a, int len, int n) {
2130         int nInts = n >>> 5;
2131         int nBits = n&0x1F;
2132         int bitsInHighWord = bitLengthForInt(a[0]);
2133 
2134         // If shift can be done without recopy, do so
2135         if (n <= (32-bitsInHighWord)) {
2136             primitiveLeftShift(a, len, nBits);
2137             return a;
2138         } else { // Array must be resized
2139             if (nBits <= (32-bitsInHighWord)) {
2140                 int result[] = new int[nInts+len];
2141                 System.arraycopy(a, 0, result, 0, len);
2142                 primitiveLeftShift(result, result.length, nBits);
2143                 return result;
2144             } else {
2145                 int result[] = new int[nInts+len+1];
2146                 System.arraycopy(a, 0, result, 0, len);
2147                 primitiveRightShift(result, result.length, 32 - nBits);
2148                 return result;
2149             }
2150         }
2151     }
2152 
2153     // shifts a up to len right n bits assumes no leading zeros, 0<n<32
2154     static void primitiveRightShift(int[] a, int len, int n) {
2155         int n2 = 32 - n;
2156         for (int i=len-1, c=a[i]; i>0; i--) {
2157             int b = c;
2158             c = a[i-1];
2159             a[i] = (c << n2) | (b >>> n);
2160         }
2161         a[0] >>>= n;
2162     }
2163 
2164     // shifts a up to len left n bits assumes no leading zeros, 0<=n<32
2165     static void primitiveLeftShift(int[] a, int len, int n) {
2166         if (len == 0 || n == 0)
2167             return;
2168 
2169         int n2 = 32 - n;
2170         for (int i=0, c=a[i], m=i+len-1; i<m; i++) {
2171             int b = c;
2172             c = a[i+1];
2173             a[i] = (b << n) | (c >>> n2);
2174         }
2175         a[len-1] <<= n;
2176     }
2177 
2178     /**
2179      * Calculate bitlength of contents of the first len elements an int array,
2180      * assuming there are no leading zero ints.
2181      */
2182     private static int bitLength(int[] val, int len) {
2183         if (len == 0)
2184             return 0;
2185         return ((len - 1) << 5) + bitLengthForInt(val[0]);
2186     }
2187 
2188     /**
2189      * Returns a BigInteger whose value is the absolute value of this
2190      * BigInteger.
2191      *
2192      * @return {@code abs(this)}
2193      */
2194     public BigInteger abs() {
2195         return (signum >= 0 ? this : this.negate());
2196     }
2197 
2198     /**
2199      * Returns a BigInteger whose value is {@code (-this)}.
2200      *
2201      * @return {@code -this}
2202      */
2203     public BigInteger negate() {
2204         return new BigInteger(this.mag, -this.signum);
2205     }
2206 
2207     /**
2208      * Returns the signum function of this BigInteger.
2209      *
2210      * @return -1, 0 or 1 as the value of this BigInteger is negative, zero or
2211      *         positive.
2212      */
2213     public int signum() {
2214         return this.signum;
2215     }
2216 
2217     // Modular Arithmetic Operations
2218 
2219     /**
2220      * Returns a BigInteger whose value is {@code (this mod m}).  This method
2221      * differs from {@code remainder} in that it always returns a
2222      * <i>non-negative</i> BigInteger.
2223      *
2224      * @param  m the modulus.
2225      * @return {@code this mod m}
2226      * @throws ArithmeticException {@code m} &le; 0
2227      * @see    #remainder
2228      */
2229     public BigInteger mod(BigInteger m) {
2230         if (m.signum <= 0)
2231             throw new ArithmeticException("BigInteger: modulus not positive");
2232 
2233         BigInteger result = this.remainder(m);
2234         return (result.signum >= 0 ? result : result.add(m));
2235     }
2236 
2237     /**
2238      * Returns a BigInteger whose value is
2239      * <tt>(this<sup>exponent</sup> mod m)</tt>.  (Unlike {@code pow}, this
2240      * method permits negative exponents.)
2241      *
2242      * @param  exponent the exponent.
2243      * @param  m the modulus.
2244      * @return <tt>this<sup>exponent</sup> mod m</tt>
2245      * @throws ArithmeticException {@code m} &le; 0 or the exponent is
2246      *         negative and this BigInteger is not <i>relatively
2247      *         prime</i> to {@code m}.
2248      * @see    #modInverse
2249      */
2250     public BigInteger modPow(BigInteger exponent, BigInteger m) {
2251         if (m.signum <= 0)
2252             throw new ArithmeticException("BigInteger: modulus not positive");
2253 
2254         // Trivial cases
2255         if (exponent.signum == 0)
2256             return (m.equals(ONE) ? ZERO : ONE);
2257 
2258         if (this.equals(ONE))
2259             return (m.equals(ONE) ? ZERO : ONE);
2260 
2261         if (this.equals(ZERO) && exponent.signum >= 0)
2262             return ZERO;
2263 
2264         if (this.equals(negConst[1]) && (!exponent.testBit(0)))
2265             return (m.equals(ONE) ? ZERO : ONE);
2266 
2267         boolean invertResult;
2268         if ((invertResult = (exponent.signum < 0)))
2269             exponent = exponent.negate();
2270 
2271         BigInteger base = (this.signum < 0 || this.compareTo(m) >= 0
2272                            ? this.mod(m) : this);
2273         BigInteger result;
2274         if (m.testBit(0)) { // odd modulus
2275             result = base.oddModPow(exponent, m);
2276         } else {
2277             /*
2278              * Even modulus.  Tear it into an "odd part" (m1) and power of two
2279              * (m2), exponentiate mod m1, manually exponentiate mod m2, and
2280              * use Chinese Remainder Theorem to combine results.
2281              */
2282 
2283             // Tear m apart into odd part (m1) and power of 2 (m2)
2284             int p = m.getLowestSetBit();   // Max pow of 2 that divides m
2285 
2286             BigInteger m1 = m.shiftRight(p);  // m/2**p
2287             BigInteger m2 = ONE.shiftLeft(p); // 2**p
2288 
2289             // Calculate new base from m1
2290             BigInteger base2 = (this.signum < 0 || this.compareTo(m1) >= 0
2291                                 ? this.mod(m1) : this);
2292 
2293             // Caculate (base ** exponent) mod m1.
2294             BigInteger a1 = (m1.equals(ONE) ? ZERO :
2295                              base2.oddModPow(exponent, m1));
2296 
2297             // Calculate (this ** exponent) mod m2
2298             BigInteger a2 = base.modPow2(exponent, p);
2299 
2300             // Combine results using Chinese Remainder Theorem
2301             BigInteger y1 = m2.modInverse(m1);
2302             BigInteger y2 = m1.modInverse(m2);
2303 
2304             result = a1.multiply(m2).multiply(y1).add
2305                      (a2.multiply(m1).multiply(y2)).mod(m);
2306         }
2307 
2308         return (invertResult ? result.modInverse(m) : result);
2309     }
2310 
2311     static int[] bnExpModThreshTable = {7, 25, 81, 241, 673, 1793,
2312                                                 Integer.MAX_VALUE}; // Sentinel
2313 
2314     /**
2315      * Returns a BigInteger whose value is x to the power of y mod z.
2316      * Assumes: z is odd && x < z.
2317      */
2318     private BigInteger oddModPow(BigInteger y, BigInteger z) {
2319     /*
2320      * The algorithm is adapted from Colin Plumb's C library.
2321      *
2322      * The window algorithm:
2323      * The idea is to keep a running product of b1 = n^(high-order bits of exp)
2324      * and then keep appending exponent bits to it.  The following patterns
2325      * apply to a 3-bit window (k = 3):
2326      * To append   0: square
2327      * To append   1: square, multiply by n^1
2328      * To append  10: square, multiply by n^1, square
2329      * To append  11: square, square, multiply by n^3
2330      * To append 100: square, multiply by n^1, square, square
2331      * To append 101: square, square, square, multiply by n^5
2332      * To append 110: square, square, multiply by n^3, square
2333      * To append 111: square, square, square, multiply by n^7
2334      *
2335      * Since each pattern involves only one multiply, the longer the pattern
2336      * the better, except that a 0 (no multiplies) can be appended directly.
2337      * We precompute a table of odd powers of n, up to 2^k, and can then
2338      * multiply k bits of exponent at a time.  Actually, assuming random
2339      * exponents, there is on average one zero bit between needs to
2340      * multiply (1/2 of the time there's none, 1/4 of the time there's 1,
2341      * 1/8 of the time, there's 2, 1/32 of the time, there's 3, etc.), so
2342      * you have to do one multiply per k+1 bits of exponent.
2343      *
2344      * The loop walks down the exponent, squaring the result buffer as
2345      * it goes.  There is a wbits+1 bit lookahead buffer, buf, that is
2346      * filled with the upcoming exponent bits.  (What is read after the
2347      * end of the exponent is unimportant, but it is filled with zero here.)
2348      * When the most-significant bit of this buffer becomes set, i.e.
2349      * (buf & tblmask) != 0, we have to decide what pattern to multiply
2350      * by, and when to do it.  We decide, remember to do it in future
2351      * after a suitable number of squarings have passed (e.g. a pattern
2352      * of "100" in the buffer requires that we multiply by n^1 immediately;
2353      * a pattern of "110" calls for multiplying by n^3 after one more
2354      * squaring), clear the buffer, and continue.
2355      *
2356      * When we start, there is one more optimization: the result buffer
2357      * is implcitly one, so squaring it or multiplying by it can be
2358      * optimized away.  Further, if we start with a pattern like "100"
2359      * in the lookahead window, rather than placing n into the buffer
2360      * and then starting to square it, we have already computed n^2
2361      * to compute the odd-powers table, so we can place that into
2362      * the buffer and save a squaring.
2363      *
2364      * This means that if you have a k-bit window, to compute n^z,
2365      * where z is the high k bits of the exponent, 1/2 of the time
2366      * it requires no squarings.  1/4 of the time, it requires 1
2367      * squaring, ... 1/2^(k-1) of the time, it reqires k-2 squarings.
2368      * And the remaining 1/2^(k-1) of the time, the top k bits are a
2369      * 1 followed by k-1 0 bits, so it again only requires k-2
2370      * squarings, not k-1.  The average of these is 1.  Add that
2371      * to the one squaring we have to do to compute the table,
2372      * and you'll see that a k-bit window saves k-2 squarings
2373      * as well as reducing the multiplies.  (It actually doesn't
2374      * hurt in the case k = 1, either.)
2375      */
2376         // Special case for exponent of one
2377         if (y.equals(ONE))
2378             return this;
2379 
2380         // Special case for base of zero
2381         if (signum==0)
2382             return ZERO;
2383 
2384         int[] base = mag.clone();
2385         int[] exp = y.mag;
2386         int[] mod = z.mag;
2387         int modLen = mod.length;
2388 
2389         // Select an appropriate window size
2390         int wbits = 0;
2391         int ebits = bitLength(exp, exp.length);
2392         // if exponent is 65537 (0x10001), use minimum window size
2393         if ((ebits != 17) || (exp[0] != 65537)) {
2394             while (ebits > bnExpModThreshTable[wbits]) {
2395                 wbits++;
2396             }
2397         }
2398 
2399         // Calculate appropriate table size
2400         int tblmask = 1 << wbits;
2401 
2402         // Allocate table for precomputed odd powers of base in Montgomery form
2403         int[][] table = new int[tblmask][];
2404         for (int i=0; i<tblmask; i++)
2405             table[i] = new int[modLen];
2406 
2407         // Compute the modular inverse
2408         int inv = -MutableBigInteger.inverseMod32(mod[modLen-1]);
2409 
2410         // Convert base to Montgomery form
2411         int[] a = leftShift(base, base.length, modLen << 5);
2412 
2413         MutableBigInteger q = new MutableBigInteger(),
2414                           a2 = new MutableBigInteger(a),
2415                           b2 = new MutableBigInteger(mod);
2416 
2417         MutableBigInteger r= a2.divide(b2, q);
2418         table[0] = r.toIntArray();
2419 
2420         // Pad table[0] with leading zeros so its length is at least modLen
2421         if (table[0].length < modLen) {
2422            int offset = modLen - table[0].length;
2423            int[] t2 = new int[modLen];
2424            for (int i=0; i<table[0].length; i++)
2425                t2[i+offset] = table[0][i];
2426            table[0] = t2;
2427         }
2428 
2429         // Set b to the square of the base
2430         int[] b = squareToLen(table[0], modLen, null);
2431         b = montReduce(b, mod, modLen, inv);
2432 
2433         // Set t to high half of b
2434         int[] t = Arrays.copyOf(b, modLen);
2435 
2436         // Fill in the table with odd powers of the base
2437         for (int i=1; i<tblmask; i++) {
2438             int[] prod = multiplyToLen(t, modLen, table[i-1], modLen, null);
2439             table[i] = montReduce(prod, mod, modLen, inv);
2440         }
2441 
2442         // Pre load the window that slides over the exponent
2443         int bitpos = 1 << ((ebits-1) & (32-1));
2444 
2445         int buf = 0;
2446         int elen = exp.length;
2447         int eIndex = 0;
2448         for (int i = 0; i <= wbits; i++) {
2449             buf = (buf << 1) | (((exp[eIndex] & bitpos) != 0)?1:0);
2450             bitpos >>>= 1;
2451             if (bitpos == 0) {
2452                 eIndex++;
2453                 bitpos = 1 << (32-1);
2454                 elen--;
2455             }
2456         }
2457 
2458         int multpos = ebits;
2459 
2460         // The first iteration, which is hoisted out of the main loop
2461         ebits--;
2462         boolean isone = true;
2463 
2464         multpos = ebits - wbits;
2465         while ((buf & 1) == 0) {
2466             buf >>>= 1;
2467             multpos++;
2468         }
2469 
2470         int[] mult = table[buf >>> 1];
2471 
2472         buf = 0;
2473         if (multpos == ebits)
2474             isone = false;
2475 
2476         // The main loop
2477         while(true) {
2478             ebits--;
2479             // Advance the window
2480             buf <<= 1;
2481 
2482             if (elen != 0) {
2483                 buf |= ((exp[eIndex] & bitpos) != 0) ? 1 : 0;
2484                 bitpos >>>= 1;
2485                 if (bitpos == 0) {
2486                     eIndex++;
2487                     bitpos = 1 << (32-1);
2488                     elen--;
2489                 }
2490             }
2491 
2492             // Examine the window for pending multiplies
2493             if ((buf & tblmask) != 0) {
2494                 multpos = ebits - wbits;
2495                 while ((buf & 1) == 0) {
2496                     buf >>>= 1;
2497                     multpos++;
2498                 }
2499                 mult = table[buf >>> 1];
2500                 buf = 0;
2501             }
2502 
2503             // Perform multiply
2504             if (ebits == multpos) {
2505                 if (isone) {
2506                     b = mult.clone();
2507                     isone = false;
2508                 } else {
2509                     t = b;
2510                     a = multiplyToLen(t, modLen, mult, modLen, a);
2511                     a = montReduce(a, mod, modLen, inv);
2512                     t = a; a = b; b = t;
2513                 }
2514             }
2515 
2516             // Check if done
2517             if (ebits == 0)
2518                 break;
2519 
2520             // Square the input
2521             if (!isone) {
2522                 t = b;
2523                 a = squareToLen(t, modLen, a);
2524                 a = montReduce(a, mod, modLen, inv);
2525                 t = a; a = b; b = t;
2526             }
2527         }
2528 
2529         // Convert result out of Montgomery form and return
2530         int[] t2 = new int[2*modLen];
2531         System.arraycopy(b, 0, t2, modLen, modLen);
2532 
2533         b = montReduce(t2, mod, modLen, inv);
2534 
2535         t2 = Arrays.copyOf(b, modLen);
2536 
2537         return new BigInteger(1, t2);
2538     }
2539 
2540     /**
2541      * Montgomery reduce n, modulo mod.  This reduces modulo mod and divides
2542      * by 2^(32*mlen). Adapted from Colin Plumb's C library.
2543      */
2544     private static int[] montReduce(int[] n, int[] mod, int mlen, int inv) {
2545         int c=0;
2546         int len = mlen;
2547         int offset=0;
2548 
2549         do {
2550             int nEnd = n[n.length-1-offset];
2551             int carry = mulAdd(n, mod, offset, mlen, inv * nEnd);
2552             c += addOne(n, offset, mlen, carry);
2553             offset++;
2554         } while(--len > 0);
2555 
2556         while(c>0)
2557             c += subN(n, mod, mlen);
2558 
2559         while (intArrayCmpToLen(n, mod, mlen) >= 0)
2560             subN(n, mod, mlen);
2561 
2562         return n;
2563     }
2564 
2565 
2566     /*
2567      * Returns -1, 0 or +1 as big-endian unsigned int array arg1 is less than,
2568      * equal to, or greater than arg2 up to length len.
2569      */
2570     private static int intArrayCmpToLen(int[] arg1, int[] arg2, int len) {
2571         for (int i=0; i<len; i++) {
2572             long b1 = arg1[i] & LONG_MASK;
2573             long b2 = arg2[i] & LONG_MASK;
2574             if (b1 < b2)
2575                 return -1;
2576             if (b1 > b2)
2577                 return 1;
2578         }
2579         return 0;
2580     }
2581 
2582     /**
2583      * Subtracts two numbers of same length, returning borrow.
2584      */
2585     private static int subN(int[] a, int[] b, int len) {
2586         long sum = 0;
2587 
2588         while(--len >= 0) {
2589             sum = (a[len] & LONG_MASK) -
2590                  (b[len] & LONG_MASK) + (sum >> 32);
2591             a[len] = (int)sum;
2592         }
2593 
2594         return (int)(sum >> 32);
2595     }
2596 
2597     /**
2598      * Multiply an array by one word k and add to result, return the carry
2599      */
2600     static int mulAdd(int[] out, int[] in, int offset, int len, int k) {
2601         long kLong = k & LONG_MASK;
2602         long carry = 0;
2603 
2604         offset = out.length-offset - 1;
2605         for (int j=len-1; j >= 0; j--) {
2606             long product = (in[j] & LONG_MASK) * kLong +
2607                            (out[offset] & LONG_MASK) + carry;
2608             out[offset--] = (int)product;
2609             carry = product >>> 32;
2610         }
2611         return (int)carry;
2612     }
2613 
2614     /**
2615      * Add one word to the number a mlen words into a. Return the resulting
2616      * carry.
2617      */
2618     static int addOne(int[] a, int offset, int mlen, int carry) {
2619         offset = a.length-1-mlen-offset;
2620         long t = (a[offset] & LONG_MASK) + (carry & LONG_MASK);
2621 
2622         a[offset] = (int)t;
2623         if ((t >>> 32) == 0)
2624             return 0;
2625         while (--mlen >= 0) {
2626             if (--offset < 0) { // Carry out of number
2627                 return 1;
2628             } else {
2629                 a[offset]++;
2630                 if (a[offset] != 0)
2631                     return 0;
2632             }
2633         }
2634         return 1;
2635     }
2636 
2637     /**
2638      * Returns a BigInteger whose value is (this ** exponent) mod (2**p)
2639      */
2640     private BigInteger modPow2(BigInteger exponent, int p) {
2641         /*
2642          * Perform exponentiation using repeated squaring trick, chopping off
2643          * high order bits as indicated by modulus.
2644          */
2645         BigInteger result = ONE;
2646         BigInteger baseToPow2 = this.mod2(p);
2647         int expOffset = 0;
2648 
2649         int limit = exponent.bitLength();
2650 
2651         if (this.testBit(0))
2652            limit = (p-1) < limit ? (p-1) : limit;
2653 
2654         while (expOffset < limit) {
2655             if (exponent.testBit(expOffset))
2656                 result = result.multiply(baseToPow2).mod2(p);
2657             expOffset++;
2658             if (expOffset < limit)
2659                 baseToPow2 = baseToPow2.square().mod2(p);
2660         }
2661 
2662         return result;
2663     }
2664 
2665     /**
2666      * Returns a BigInteger whose value is this mod(2**p).
2667      * Assumes that this {@code BigInteger >= 0} and {@code p > 0}.
2668      */
2669     private BigInteger mod2(int p) {
2670         if (bitLength() <= p)
2671             return this;
2672 
2673         // Copy remaining ints of mag
2674         int numInts = (p + 31) >>> 5;
2675         int[] mag = new int[numInts];
2676         System.arraycopy(this.mag, (this.mag.length - numInts), mag, 0, numInts);
2677 
2678         // Mask out any excess bits
2679         int excessBits = (numInts << 5) - p;
2680         mag[0] &= (1L << (32-excessBits)) - 1;
2681 
2682         return (mag[0]==0 ? new BigInteger(1, mag) : new BigInteger(mag, 1));
2683     }
2684 
2685     /**
2686      * Returns a BigInteger whose value is {@code (this}<sup>-1</sup> {@code mod m)}.
2687      *
2688      * @param  m the modulus.
2689      * @return {@code this}<sup>-1</sup> {@code mod m}.
2690      * @throws ArithmeticException {@code  m} &le; 0, or this BigInteger
2691      *         has no multiplicative inverse mod m (that is, this BigInteger
2692      *         is not <i>relatively prime</i> to m).
2693      */
2694     public BigInteger modInverse(BigInteger m) {
2695         if (m.signum != 1)
2696             throw new ArithmeticException("BigInteger: modulus not positive");
2697 
2698         if (m.equals(ONE))
2699             return ZERO;
2700 
2701         // Calculate (this mod m)
2702         BigInteger modVal = this;
2703         if (signum < 0 || (this.compareMagnitude(m) >= 0))
2704             modVal = this.mod(m);
2705 
2706         if (modVal.equals(ONE))
2707             return ONE;
2708 
2709         MutableBigInteger a = new MutableBigInteger(modVal);
2710         MutableBigInteger b = new MutableBigInteger(m);
2711 
2712         MutableBigInteger result = a.mutableModInverse(b);
2713         return result.toBigInteger(1);
2714     }
2715 
2716     // Shift Operations
2717 
2718     /**
2719      * Returns a BigInteger whose value is {@code (this << n)}.
2720      * The shift distance, {@code n}, may be negative, in which case
2721      * this method performs a right shift.
2722      * (Computes <tt>floor(this * 2<sup>n</sup>)</tt>.)
2723      *
2724      * @param  n shift distance, in bits.
2725      * @return {@code this << n}
2726      * @throws ArithmeticException if the shift distance is {@code
2727      *         Integer.MIN_VALUE}.
2728      * @see #shiftRight
2729      */
2730     public BigInteger shiftLeft(int n) {
2731         if (signum == 0)
2732             return ZERO;
2733         if (n==0)
2734             return this;
2735         if (n<0) {
2736             if (n == Integer.MIN_VALUE) {
2737                 throw new ArithmeticException("Shift distance of Integer.MIN_VALUE not supported.");
2738             } else {
2739                 return shiftRight(-n);
2740             }
2741         }
2742         int[] newMag = shiftLeft(mag, n);
2743 
2744         return new BigInteger(newMag, signum);
2745     }
2746 
2747     private static int[] shiftLeft(int[] mag, int n) {
2748         int nInts = n >>> 5;
2749         int nBits = n & 0x1f;
2750         int magLen = mag.length;
2751         int newMag[] = null;
2752 
2753         if (nBits == 0) {
2754             newMag = new int[magLen + nInts];
2755             System.arraycopy(mag, 0, newMag, 0, magLen);
2756         } else {
2757             int i = 0;
2758             int nBits2 = 32 - nBits;
2759             int highBits = mag[0] >>> nBits2;
2760             if (highBits != 0) {
2761                 newMag = new int[magLen + nInts + 1];
2762                 newMag[i++] = highBits;
2763             } else {
2764                 newMag = new int[magLen + nInts];
2765             }
2766             int j=0;
2767             while (j < magLen-1)
2768                 newMag[i++] = mag[j++] << nBits | mag[j] >>> nBits2;
2769             newMag[i] = mag[j] << nBits;
2770         }
2771         return newMag;
2772     }
2773 
2774     /**
2775      * Returns a BigInteger whose value is {@code (this >> n)}.  Sign
2776      * extension is performed.  The shift distance, {@code n}, may be
2777      * negative, in which case this method performs a left shift.
2778      * (Computes <tt>floor(this / 2<sup>n</sup>)</tt>.)
2779      *
2780      * @param  n shift distance, in bits.
2781      * @return {@code this >> n}
2782      * @throws ArithmeticException if the shift distance is {@code
2783      *         Integer.MIN_VALUE}.
2784      * @see #shiftLeft
2785      */
2786     public BigInteger shiftRight(int n) {
2787         if (n==0)
2788             return this;
2789         if (n<0) {
2790             if (n == Integer.MIN_VALUE) {
2791                 throw new ArithmeticException("Shift distance of Integer.MIN_VALUE not supported.");
2792             } else {
2793                 return shiftLeft(-n);
2794             }
2795         }
2796 
2797         int nInts = n >>> 5;
2798         int nBits = n & 0x1f;
2799         int magLen = mag.length;
2800         int newMag[] = null;
2801 
2802         // Special case: entire contents shifted off the end
2803         if (nInts >= magLen)
2804             return (signum >= 0 ? ZERO : negConst[1]);
2805 
2806         if (nBits == 0) {
2807             int newMagLen = magLen - nInts;
2808             newMag = Arrays.copyOf(mag, newMagLen);
2809         } else {
2810             int i = 0;
2811             int highBits = mag[0] >>> nBits;
2812             if (highBits != 0) {
2813                 newMag = new int[magLen - nInts];
2814                 newMag[i++] = highBits;
2815             } else {
2816                 newMag = new int[magLen - nInts -1];
2817             }
2818 
2819             int nBits2 = 32 - nBits;
2820             int j=0;
2821             while (j < magLen - nInts - 1)
2822                 newMag[i++] = (mag[j++] << nBits2) | (mag[j] >>> nBits);
2823         }
2824 
2825         if (signum < 0) {
2826             // Find out whether any one-bits were shifted off the end.
2827             boolean onesLost = false;
2828             for (int i=magLen-1, j=magLen-nInts; i>=j && !onesLost; i--)
2829                 onesLost = (mag[i] != 0);
2830             if (!onesLost && nBits != 0)
2831                 onesLost = (mag[magLen - nInts - 1] << (32 - nBits) != 0);
2832 
2833             if (onesLost)
2834                 newMag = javaIncrement(newMag);
2835         }
2836 
2837         return new BigInteger(newMag, signum);
2838     }
2839 
2840     int[] javaIncrement(int[] val) {
2841         int lastSum = 0;
2842         for (int i=val.length-1;  i >= 0 && lastSum == 0; i--)
2843             lastSum = (val[i] += 1);
2844         if (lastSum == 0) {
2845             val = new int[val.length+1];
2846             val[0] = 1;
2847         }
2848         return val;
2849     }
2850 
2851     // Bitwise Operations
2852 
2853     /**
2854      * Returns a BigInteger whose value is {@code (this & val)}.  (This
2855      * method returns a negative BigInteger if and only if this and val are
2856      * both negative.)
2857      *
2858      * @param val value to be AND'ed with this BigInteger.
2859      * @return {@code this & val}
2860      */
2861     public BigInteger and(BigInteger val) {
2862         int[] result = new int[Math.max(intLength(), val.intLength())];
2863         for (int i=0; i<result.length; i++)
2864             result[i] = (getInt(result.length-i-1)
2865                          & val.getInt(result.length-i-1));
2866 
2867         return valueOf(result);
2868     }
2869 
2870     /**
2871      * Returns a BigInteger whose value is {@code (this | val)}.  (This method
2872      * returns a negative BigInteger if and only if either this or val is
2873      * negative.)
2874      *
2875      * @param val value to be OR'ed with this BigInteger.
2876      * @return {@code this | val}
2877      */
2878     public BigInteger or(BigInteger val) {
2879         int[] result = new int[Math.max(intLength(), val.intLength())];
2880         for (int i=0; i<result.length; i++)
2881             result[i] = (getInt(result.length-i-1)
2882                          | val.getInt(result.length-i-1));
2883 
2884         return valueOf(result);
2885     }
2886 
2887     /**
2888      * Returns a BigInteger whose value is {@code (this ^ val)}.  (This method
2889      * returns a negative BigInteger if and only if exactly one of this and
2890      * val are negative.)
2891      *
2892      * @param val value to be XOR'ed with this BigInteger.
2893      * @return {@code this ^ val}
2894      */
2895     public BigInteger xor(BigInteger val) {
2896         int[] result = new int[Math.max(intLength(), val.intLength())];
2897         for (int i=0; i<result.length; i++)
2898             result[i] = (getInt(result.length-i-1)
2899                          ^ val.getInt(result.length-i-1));
2900 
2901         return valueOf(result);
2902     }
2903 
2904     /**
2905      * Returns a BigInteger whose value is {@code (~this)}.  (This method
2906      * returns a negative value if and only if this BigInteger is
2907      * non-negative.)
2908      *
2909      * @return {@code ~this}
2910      */
2911     public BigInteger not() {
2912         int[] result = new int[intLength()];
2913         for (int i=0; i<result.length; i++)
2914             result[i] = ~getInt(result.length-i-1);
2915 
2916         return valueOf(result);
2917     }
2918 
2919     /**
2920      * Returns a BigInteger whose value is {@code (this & ~val)}.  This
2921      * method, which is equivalent to {@code and(val.not())}, is provided as
2922      * a convenience for masking operations.  (This method returns a negative
2923      * BigInteger if and only if {@code this} is negative and {@code val} is
2924      * positive.)
2925      *
2926      * @param val value to be complemented and AND'ed with this BigInteger.
2927      * @return {@code this & ~val}
2928      */
2929     public BigInteger andNot(BigInteger val) {
2930         int[] result = new int[Math.max(intLength(), val.intLength())];
2931         for (int i=0; i<result.length; i++)
2932             result[i] = (getInt(result.length-i-1)
2933                          & ~val.getInt(result.length-i-1));
2934 
2935         return valueOf(result);
2936     }
2937 
2938 
2939     // Single Bit Operations
2940 
2941     /**
2942      * Returns {@code true} if and only if the designated bit is set.
2943      * (Computes {@code ((this & (1<<n)) != 0)}.)
2944      *
2945      * @param  n index of bit to test.
2946      * @return {@code true} if and only if the designated bit is set.
2947      * @throws ArithmeticException {@code n} is negative.
2948      */
2949     public boolean testBit(int n) {
2950         if (n<0)
2951             throw new ArithmeticException("Negative bit address");
2952 
2953         return (getInt(n >>> 5) & (1 << (n & 31))) != 0;
2954     }
2955 
2956     /**
2957      * Returns a BigInteger whose value is equivalent to this BigInteger
2958      * with the designated bit set.  (Computes {@code (this | (1<<n))}.)
2959      *
2960      * @param  n index of bit to set.
2961      * @return {@code this | (1<<n)}
2962      * @throws ArithmeticException {@code n} is negative.
2963      */
2964     public BigInteger setBit(int n) {
2965         if (n<0)
2966             throw new ArithmeticException("Negative bit address");
2967 
2968         int intNum = n >>> 5;
2969         int[] result = new int[Math.max(intLength(), intNum+2)];
2970 
2971         for (int i=0; i<result.length; i++)
2972             result[result.length-i-1] = getInt(i);
2973 
2974         result[result.length-intNum-1] |= (1 << (n & 31));
2975 
2976         return valueOf(result);
2977     }
2978 
2979     /**
2980      * Returns a BigInteger whose value is equivalent to this BigInteger
2981      * with the designated bit cleared.
2982      * (Computes {@code (this & ~(1<<n))}.)
2983      *
2984      * @param  n index of bit to clear.
2985      * @return {@code this & ~(1<<n)}
2986      * @throws ArithmeticException {@code n} is negative.
2987      */
2988     public BigInteger clearBit(int n) {
2989         if (n<0)
2990             throw new ArithmeticException("Negative bit address");
2991 
2992         int intNum = n >>> 5;
2993         int[] result = new int[Math.max(intLength(), ((n + 1) >>> 5) + 1)];
2994 
2995         for (int i=0; i<result.length; i++)
2996             result[result.length-i-1] = getInt(i);
2997 
2998         result[result.length-intNum-1] &= ~(1 << (n & 31));
2999 
3000         return valueOf(result);
3001     }
3002 
3003     /**
3004      * Returns a BigInteger whose value is equivalent to this BigInteger
3005      * with the designated bit flipped.
3006      * (Computes {@code (this ^ (1<<n))}.)
3007      *
3008      * @param  n index of bit to flip.
3009      * @return {@code this ^ (1<<n)}
3010      * @throws ArithmeticException {@code n} is negative.
3011      */
3012     public BigInteger flipBit(int n) {
3013         if (n<0)
3014             throw new ArithmeticException("Negative bit address");
3015 
3016         int intNum = n >>> 5;
3017         int[] result = new int[Math.max(intLength(), intNum+2)];
3018 
3019         for (int i=0; i<result.length; i++)
3020             result[result.length-i-1] = getInt(i);
3021 
3022         result[result.length-intNum-1] ^= (1 << (n & 31));
3023 
3024         return valueOf(result);
3025     }
3026 
3027     /**
3028      * Returns the index of the rightmost (lowest-order) one bit in this
3029      * BigInteger (the number of zero bits to the right of the rightmost
3030      * one bit).  Returns -1 if this BigInteger contains no one bits.
3031      * (Computes {@code (this==0? -1 : log2(this & -this))}.)
3032      *
3033      * @return index of the rightmost one bit in this BigInteger.
3034      */
3035     public int getLowestSetBit() {
3036         @SuppressWarnings("deprecation") int lsb = lowestSetBit - 2;
3037         if (lsb == -2) {  // lowestSetBit not initialized yet
3038             lsb = 0;
3039             if (signum == 0) {
3040                 lsb -= 1;
3041             } else {
3042                 // Search for lowest order nonzero int
3043                 int i,b;
3044                 for (i=0; (b = getInt(i))==0; i++)
3045                     ;
3046                 lsb += (i << 5) + Integer.numberOfTrailingZeros(b);
3047             }
3048             lowestSetBit = lsb + 2;
3049         }
3050         return lsb;
3051     }
3052 
3053 
3054     // Miscellaneous Bit Operations
3055 
3056     /**
3057      * Returns the number of bits in the minimal two's-complement
3058      * representation of this BigInteger, <i>excluding</i> a sign bit.
3059      * For positive BigIntegers, this is equivalent to the number of bits in
3060      * the ordinary binary representation.  (Computes
3061      * {@code (ceil(log2(this < 0 ? -this : this+1)))}.)
3062      *
3063      * @return number of bits in the minimal two's-complement
3064      *         representation of this BigInteger, <i>excluding</i> a sign bit.
3065      */
3066     public int bitLength() {
3067         @SuppressWarnings("deprecation") int n = bitLength - 1;
3068         if (n == -1) { // bitLength not initialized yet
3069             int[] m = mag;
3070             int len = m.length;
3071             if (len == 0) {
3072                 n = 0; // offset by one to initialize
3073             }  else {
3074                 // Calculate the bit length of the magnitude
3075                 int magBitLength = ((len - 1) << 5) + bitLengthForInt(mag[0]);
3076                  if (signum < 0) {
3077                      // Check if magnitude is a power of two
3078                      boolean pow2 = (Integer.bitCount(mag[0]) == 1);
3079                      for (int i=1; i< len && pow2; i++)
3080                          pow2 = (mag[i] == 0);
3081 
3082                      n = (pow2 ? magBitLength -1 : magBitLength);
3083                  } else {
3084                      n = magBitLength;
3085                  }
3086             }
3087             bitLength = n + 1;
3088         }
3089         return n;
3090     }
3091 
3092     /**
3093      * Returns the number of bits in the two's complement representation
3094      * of this BigInteger that differ from its sign bit.  This method is
3095      * useful when implementing bit-vector style sets atop BigIntegers.
3096      *
3097      * @return number of bits in the two's complement representation
3098      *         of this BigInteger that differ from its sign bit.
3099      */
3100     public int bitCount() {
3101         @SuppressWarnings("deprecation") int bc = bitCount - 1;
3102         if (bc == -1) {  // bitCount not initialized yet
3103             bc = 0;      // offset by one to initialize
3104             // Count the bits in the magnitude
3105             for (int i=0; i<mag.length; i++)
3106                 bc += Integer.bitCount(mag[i]);
3107             if (signum < 0) {
3108                 // Count the trailing zeros in the magnitude
3109                 int magTrailingZeroCount = 0, j;
3110                 for (j=mag.length-1; mag[j]==0; j--)
3111                     magTrailingZeroCount += 32;
3112                 magTrailingZeroCount += Integer.numberOfTrailingZeros(mag[j]);
3113                 bc += magTrailingZeroCount - 1;
3114             }
3115             bitCount = bc + 1;
3116         }
3117         return bc;
3118     }
3119 
3120     // Primality Testing
3121 
3122     /**
3123      * Returns {@code true} if this BigInteger is probably prime,
3124      * {@code false} if it's definitely composite.  If
3125      * {@code certainty} is &le; 0, {@code true} is
3126      * returned.
3127      *
3128      * @param  certainty a measure of the uncertainty that the caller is
3129      *         willing to tolerate: if the call returns {@code true}
3130      *         the probability that this BigInteger is prime exceeds
3131      *         (1 - 1/2<sup>{@code certainty}</sup>).  The execution time of
3132      *         this method is proportional to the value of this parameter.
3133      * @return {@code true} if this BigInteger is probably prime,
3134      *         {@code false} if it's definitely composite.
3135      */
3136     public boolean isProbablePrime(int certainty) {
3137         if (certainty <= 0)
3138             return true;
3139         BigInteger w = this.abs();
3140         if (w.equals(TWO))
3141             return true;
3142         if (!w.testBit(0) || w.equals(ONE))
3143             return false;
3144 
3145         return w.primeToCertainty(certainty, null);
3146     }
3147 
3148     // Comparison Operations
3149 
3150     /**
3151      * Compares this BigInteger with the specified BigInteger.  This
3152      * method is provided in preference to individual methods for each
3153      * of the six boolean comparison operators ({@literal <}, ==,
3154      * {@literal >}, {@literal >=}, !=, {@literal <=}).  The suggested
3155      * idiom for performing these comparisons is: {@code
3156      * (x.compareTo(y)} &lt;<i>op</i>&gt; {@code 0)}, where
3157      * &lt;<i>op</i>&gt; is one of the six comparison operators.
3158      *
3159      * @param  val BigInteger to which this BigInteger is to be compared.
3160      * @return -1, 0 or 1 as this BigInteger is numerically less than, equal
3161      *         to, or greater than {@code val}.
3162      */
3163     public int compareTo(BigInteger val) {
3164         if (signum == val.signum) {
3165             switch (signum) {
3166             case 1:
3167                 return compareMagnitude(val);
3168             case -1:
3169                 return val.compareMagnitude(this);
3170             default:
3171                 return 0;
3172             }
3173         }
3174         return signum > val.signum ? 1 : -1;
3175     }
3176 
3177     /**
3178      * Compares the magnitude array of this BigInteger with the specified
3179      * BigInteger's. This is the version of compareTo ignoring sign.
3180      *
3181      * @param val BigInteger whose magnitude array to be compared.
3182      * @return -1, 0 or 1 as this magnitude array is less than, equal to or
3183      *         greater than the magnitude aray for the specified BigInteger's.
3184      */
3185     final int compareMagnitude(BigInteger val) {
3186         int[] m1 = mag;
3187         int len1 = m1.length;
3188         int[] m2 = val.mag;
3189         int len2 = m2.length;
3190         if (len1 < len2)
3191             return -1;
3192         if (len1 > len2)
3193             return 1;
3194         for (int i = 0; i < len1; i++) {
3195             int a = m1[i];
3196             int b = m2[i];
3197             if (a != b)
3198                 return ((a & LONG_MASK) < (b & LONG_MASK)) ? -1 : 1;
3199         }
3200         return 0;
3201     }
3202 
3203     /**
3204      * Version of compareMagnitude that compares magnitude with long value.
3205      * val can't be Long.MIN_VALUE.
3206      */
3207     final int compareMagnitude(long val) {
3208         assert val != Long.MIN_VALUE;
3209         int[] m1 = mag;
3210         int len = m1.length;
3211         if(len > 2) {
3212             return 1;
3213         }
3214         if (val < 0) {
3215             val = -val;
3216         }
3217         int highWord = (int)(val >>> 32);
3218         if (highWord==0) {
3219             if (len < 1)
3220                 return -1;
3221             if (len > 1)
3222                 return 1;
3223             int a = m1[0];
3224             int b = (int)val;
3225             if (a != b) {
3226                 return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
3227             }
3228             return 0;
3229         } else {
3230             if (len < 2)
3231                 return -1;
3232             int a = m1[0];
3233             int b = highWord;
3234             if (a != b) {
3235                 return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
3236             }
3237             a = m1[1];
3238             b = (int)val;
3239             if (a != b) {
3240                 return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
3241             }
3242             return 0;
3243         }
3244     }
3245 
3246     /**
3247      * Compares this BigInteger with the specified Object for equality.
3248      *
3249      * @param  x Object to which this BigInteger is to be compared.
3250      * @return {@code true} if and only if the specified Object is a
3251      *         BigInteger whose value is numerically equal to this BigInteger.
3252      */
3253     public boolean equals(Object x) {
3254         // This test is just an optimization, which may or may not help
3255         if (x == this)
3256             return true;
3257 
3258         if (!(x instanceof BigInteger))
3259             return false;
3260 
3261         BigInteger xInt = (BigInteger) x;
3262         if (xInt.signum != signum)
3263             return false;
3264 
3265         int[] m = mag;
3266         int len = m.length;
3267         int[] xm = xInt.mag;
3268         if (len != xm.length)
3269             return false;
3270 
3271         for (int i = 0; i < len; i++)
3272             if (xm[i] != m[i])
3273                 return false;
3274 
3275         return true;
3276     }
3277 
3278     /**
3279      * Returns the minimum of this BigInteger and {@code val}.
3280      *
3281      * @param  val value with which the minimum is to be computed.
3282      * @return the BigInteger whose value is the lesser of this BigInteger and
3283      *         {@code val}.  If they are equal, either may be returned.
3284      */
3285     public BigInteger min(BigInteger val) {
3286         return (compareTo(val)<0 ? this : val);
3287     }
3288 
3289     /**
3290      * Returns the maximum of this BigInteger and {@code val}.
3291      *
3292      * @param  val value with which the maximum is to be computed.
3293      * @return the BigInteger whose value is the greater of this and
3294      *         {@code val}.  If they are equal, either may be returned.
3295      */
3296     public BigInteger max(BigInteger val) {
3297         return (compareTo(val)>0 ? this : val);
3298     }
3299 
3300 
3301     // Hash Function
3302 
3303     /**
3304      * Returns the hash code for this BigInteger.
3305      *
3306      * @return hash code for this BigInteger.
3307      */
3308     public int hashCode() {
3309         int hashCode = 0;
3310 
3311         for (int i=0; i<mag.length; i++)
3312             hashCode = (int)(31*hashCode + (mag[i] & LONG_MASK));
3313 
3314         return hashCode * signum;
3315     }
3316 
3317     /**
3318      * Returns the String representation of this BigInteger in the
3319      * given radix.  If the radix is outside the range from {@link
3320      * Character#MIN_RADIX} to {@link Character#MAX_RADIX} inclusive,
3321      * it will default to 10 (as is the case for
3322      * {@code Integer.toString}).  The digit-to-character mapping
3323      * provided by {@code Character.forDigit} is used, and a minus
3324      * sign is prepended if appropriate.  (This representation is
3325      * compatible with the {@link #BigInteger(String, int) (String,
3326      * int)} constructor.)
3327      *
3328      * @param  radix  radix of the String representation.
3329      * @return String representation of this BigInteger in the given radix.
3330      * @see    Integer#toString
3331      * @see    Character#forDigit
3332      * @see    #BigInteger(java.lang.String, int)
3333      */
3334     public String toString(int radix) {
3335         if (signum == 0)
3336             return "0";
3337         if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
3338             radix = 10;
3339 
3340         // If it's small enough, use smallToString.
3341         if (mag.length <= SCHOENHAGE_BASE_CONVERSION_THRESHOLD)
3342            return smallToString(radix);
3343 
3344         // Otherwise use recursive toString, which requires positive arguments.
3345         // The results will be concatenated into this StringBuilder
3346         StringBuilder sb = new StringBuilder();
3347         if (signum < 0) {
3348             toString(this.negate(), sb, radix, 0);
3349             sb.insert(0, '-');
3350         }
3351         else
3352             toString(this, sb, radix, 0);
3353 
3354         return sb.toString();
3355     }
3356 
3357     /** This method is used to perform toString when arguments are small. */
3358     private String smallToString(int radix) {
3359         if (signum == 0)
3360             return "0";
3361 
3362         // Compute upper bound on number of digit groups and allocate space
3363         int maxNumDigitGroups = (4*mag.length + 6)/7;
3364         String digitGroup[] = new String[maxNumDigitGroups];
3365 
3366         // Translate number to string, a digit group at a time
3367         BigInteger tmp = this.abs();
3368         int numGroups = 0;
3369         while (tmp.signum != 0) {
3370             BigInteger d = longRadix[radix];
3371 
3372             MutableBigInteger q = new MutableBigInteger(),
3373                               a = new MutableBigInteger(tmp.mag),
3374                               b = new MutableBigInteger(d.mag);
3375             MutableBigInteger r = a.divide(b, q);
3376             BigInteger q2 = q.toBigInteger(tmp.signum * d.signum);
3377             BigInteger r2 = r.toBigInteger(tmp.signum * d.signum);
3378 
3379             digitGroup[numGroups++] = Long.toString(r2.longValue(), radix);
3380             tmp = q2;
3381         }
3382 
3383         // Put sign (if any) and first digit group into result buffer
3384         StringBuilder buf = new StringBuilder(numGroups*digitsPerLong[radix]+1);
3385         if (signum<0)
3386             buf.append('-');
3387         buf.append(digitGroup[numGroups-1]);
3388 
3389         // Append remaining digit groups padded with leading zeros
3390         for (int i=numGroups-2; i>=0; i--) {
3391             // Prepend (any) leading zeros for this digit group
3392             int numLeadingZeros = digitsPerLong[radix]-digitGroup[i].length();
3393             if (numLeadingZeros != 0)
3394                 buf.append(zeros[numLeadingZeros]);
3395             buf.append(digitGroup[i]);
3396         }
3397         return buf.toString();
3398     }
3399 
3400     /**
3401      * Converts the specified BigInteger to a string and appends to
3402      * <code>sb</code>.  This implements the recursive Schoenhage algorithm
3403      * for base conversions.
3404      * <p/>
3405      * See Knuth, Donald,  _The Art of Computer Programming_, Vol. 2,
3406      * Answers to Exercises (4.4) Question 14.
3407      *
3408      * @param u      The number to convert to a string.
3409      * @param sb     The StringBuilder that will be appended to in place.
3410      * @param radix  The base to convert to.
3411      * @param digits The minimum number of digits to pad to.
3412      */
3413     private static void toString(BigInteger u, StringBuilder sb, int radix,
3414                                  int digits) {
3415         /* If we're smaller than a certain threshold, use the smallToString
3416            method, padding with leading zeroes when necessary. */
3417         if (u.mag.length <= SCHOENHAGE_BASE_CONVERSION_THRESHOLD) {
3418             String s = u.smallToString(radix);
3419 
3420             // Pad with internal zeros if necessary.
3421             // Don't pad if we're at the beginning of the string.
3422             if ((s.length() < digits) && (sb.length() > 0))
3423                 for (int i=s.length(); i<digits; i++) // May be a faster way to
3424                     sb.append('0');                    // do this?
3425 
3426             sb.append(s);
3427             return;
3428         }
3429 
3430         int b, n;
3431         b = u.bitLength();
3432 
3433         // Calculate a value for n in the equation radix^(2^n) = u
3434         // and subtract 1 from that value.  This is used to find the
3435         // cache index that contains the best value to divide u.
3436         n = (int) Math.round(Math.log(b * LOG_TWO / logCache[radix]) / LOG_TWO - 1.0);
3437         BigInteger v = getRadixConversionCache(radix, n);
3438         BigInteger[] results;
3439         results = u.divideAndRemainder(v);
3440 
3441         int expectedDigits = 1 << n;
3442 
3443         // Now recursively build the two halves of each number.
3444         toString(results[0], sb, radix, digits-expectedDigits);
3445         toString(results[1], sb, radix, expectedDigits);
3446     }
3447 
3448     /**
3449      * Returns the value radix^(2^exponent) from the cache.
3450      * If this value doesn't already exist in the cache, it is added.
3451      * <p/>
3452      * This could be changed to a more complicated caching method using
3453      * <code>Future</code>.
3454      */
3455     private static BigInteger getRadixConversionCache(int radix, int exponent) {
3456         BigInteger[] cacheLine = powerCache[radix]; // volatile read
3457         if (exponent < cacheLine.length) {
3458             return cacheLine[exponent];






3459         }
3460 
3461         int oldLength = cacheLine.length;
3462         cacheLine = Arrays.copyOf(cacheLine, exponent + 1);
3463         for (int i = oldLength; i <= exponent; i++) {
3464             cacheLine[i] = cacheLine[i - 1].pow(2);
3465         }


3466 
3467         BigInteger[][] pc = powerCache; // volatile read again
3468         if (exponent >= pc[radix].length) {
3469             pc = pc.clone();
3470             pc[radix] = cacheLine;
3471             powerCache = pc; // volatile write, publish
3472         }
3473         return cacheLine[exponent];
3474     }
3475 
3476     /* zero[i] is a string of i consecutive zeros. */
3477     private static String zeros[] = new String[64];
3478     static {
3479         zeros[63] =
3480             "000000000000000000000000000000000000000000000000000000000000000";
3481         for (int i=0; i<63; i++)
3482             zeros[i] = zeros[63].substring(0, i);
3483     }
3484 
3485     /**
3486      * Returns the decimal String representation of this BigInteger.
3487      * The digit-to-character mapping provided by
3488      * {@code Character.forDigit} is used, and a minus sign is
3489      * prepended if appropriate.  (This representation is compatible
3490      * with the {@link #BigInteger(String) (String)} constructor, and
3491      * allows for String concatenation with Java's + operator.)
3492      *
3493      * @return decimal String representation of this BigInteger.
3494      * @see    Character#forDigit
3495      * @see    #BigInteger(java.lang.String)
3496      */
3497     public String toString() {
3498         return toString(10);
3499     }
3500 
3501     /**
3502      * Returns a byte array containing the two's-complement
3503      * representation of this BigInteger.  The byte array will be in
3504      * <i>big-endian</i> byte-order: the most significant byte is in
3505      * the zeroth element.  The array will contain the minimum number
3506      * of bytes required to represent this BigInteger, including at
3507      * least one sign bit, which is {@code (ceil((this.bitLength() +
3508      * 1)/8))}.  (This representation is compatible with the
3509      * {@link #BigInteger(byte[]) (byte[])} constructor.)
3510      *
3511      * @return a byte array containing the two's-complement representation of
3512      *         this BigInteger.
3513      * @see    #BigInteger(byte[])
3514      */
3515     public byte[] toByteArray() {
3516         int byteLen = bitLength()/8 + 1;
3517         byte[] byteArray = new byte[byteLen];
3518 
3519         for (int i=byteLen-1, bytesCopied=4, nextInt=0, intIndex=0; i>=0; i--) {
3520             if (bytesCopied == 4) {
3521                 nextInt = getInt(intIndex++);
3522                 bytesCopied = 1;
3523             } else {
3524                 nextInt >>>= 8;
3525                 bytesCopied++;
3526             }
3527             byteArray[i] = (byte)nextInt;
3528         }
3529         return byteArray;
3530     }
3531 
3532     /**
3533      * Converts this BigInteger to an {@code int}.  This
3534      * conversion is analogous to a
3535      * <i>narrowing primitive conversion</i> from {@code long} to
3536      * {@code int} as defined in section 5.1.3 of
3537      * <cite>The Java&trade; Language Specification</cite>:
3538      * if this BigInteger is too big to fit in an
3539      * {@code int}, only the low-order 32 bits are returned.
3540      * Note that this conversion can lose information about the
3541      * overall magnitude of the BigInteger value as well as return a
3542      * result with the opposite sign.
3543      *
3544      * @return this BigInteger converted to an {@code int}.
3545      * @see #intValueExact()
3546      */
3547     public int intValue() {
3548         int result = 0;
3549         result = getInt(0);
3550         return result;
3551     }
3552 
3553     /**
3554      * Converts this BigInteger to a {@code long}.  This
3555      * conversion is analogous to a
3556      * <i>narrowing primitive conversion</i> from {@code long} to
3557      * {@code int} as defined in section 5.1.3 of
3558      * <cite>The Java&trade; Language Specification</cite>:
3559      * if this BigInteger is too big to fit in a
3560      * {@code long}, only the low-order 64 bits are returned.
3561      * Note that this conversion can lose information about the
3562      * overall magnitude of the BigInteger value as well as return a
3563      * result with the opposite sign.
3564      *
3565      * @return this BigInteger converted to a {@code long}.
3566      * @see #longValueExact()
3567      */
3568     public long longValue() {
3569         long result = 0;
3570 
3571         for (int i=1; i>=0; i--)
3572             result = (result << 32) + (getInt(i) & LONG_MASK);
3573         return result;
3574     }
3575 
3576     /**
3577      * Converts this BigInteger to a {@code float}.  This
3578      * conversion is similar to the
3579      * <i>narrowing primitive conversion</i> from {@code double} to
3580      * {@code float} as defined in section 5.1.3 of
3581      * <cite>The Java&trade; Language Specification</cite>:
3582      * if this BigInteger has too great a magnitude
3583      * to represent as a {@code float}, it will be converted to
3584      * {@link Float#NEGATIVE_INFINITY} or {@link
3585      * Float#POSITIVE_INFINITY} as appropriate.  Note that even when
3586      * the return value is finite, this conversion can lose
3587      * information about the precision of the BigInteger value.
3588      *
3589      * @return this BigInteger converted to a {@code float}.
3590      */
3591     public float floatValue() {
3592         if (signum == 0) {
3593             return 0.0f;
3594         }
3595 
3596         int exponent = ((mag.length - 1) << 5) + bitLengthForInt(mag[0]) - 1;
3597 
3598         // exponent == floor(log2(abs(this)))
3599         if (exponent < Long.SIZE - 1) {
3600             return longValue();
3601         } else if (exponent > Float.MAX_EXPONENT) {
3602             return signum > 0 ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3603         }
3604 
3605         /*
3606          * We need the top SIGNIFICAND_WIDTH bits, including the "implicit"
3607          * one bit. To make rounding easier, we pick out the top
3608          * SIGNIFICAND_WIDTH + 1 bits, so we have one to help us round up or
3609          * down. twiceSignifFloor will contain the top SIGNIFICAND_WIDTH + 1
3610          * bits, and signifFloor the top SIGNIFICAND_WIDTH.
3611          *
3612          * It helps to consider the real number signif = abs(this) *
3613          * 2^(SIGNIFICAND_WIDTH - 1 - exponent).
3614          */
3615         int shift = exponent - FloatConsts.SIGNIFICAND_WIDTH;
3616 
3617         int twiceSignifFloor;
3618         // twiceSignifFloor will be == abs().shiftRight(shift).intValue()
3619         // We do the shift into an int directly to improve performance.
3620 
3621         int nBits = shift & 0x1f;
3622         int nBits2 = 32 - nBits;
3623 
3624         if (nBits == 0) {
3625             twiceSignifFloor = mag[0];
3626         } else {
3627             twiceSignifFloor = mag[0] >>> nBits;
3628             if (twiceSignifFloor == 0) {
3629                 twiceSignifFloor = (mag[0] << nBits2) | (mag[1] >>> nBits);
3630             }
3631         }
3632 
3633         int signifFloor = twiceSignifFloor >> 1;
3634         signifFloor &= FloatConsts.SIGNIF_BIT_MASK; // remove the implied bit
3635 
3636         /*
3637          * We round up if either the fractional part of signif is strictly
3638          * greater than 0.5 (which is true if the 0.5 bit is set and any lower
3639          * bit is set), or if the fractional part of signif is >= 0.5 and
3640          * signifFloor is odd (which is true if both the 0.5 bit and the 1 bit
3641          * are set). This is equivalent to the desired HALF_EVEN rounding.
3642          */
3643         boolean increment = (twiceSignifFloor & 1) != 0
3644                 && ((signifFloor & 1) != 0 || abs().getLowestSetBit() < shift);
3645         int signifRounded = increment ? signifFloor + 1 : signifFloor;
3646         int bits = ((exponent + FloatConsts.EXP_BIAS))
3647                 << (FloatConsts.SIGNIFICAND_WIDTH - 1);
3648         bits += signifRounded;
3649         /*
3650          * If signifRounded == 2^24, we'd need to set all of the significand
3651          * bits to zero and add 1 to the exponent. This is exactly the behavior
3652          * we get from just adding signifRounded to bits directly. If the
3653          * exponent is Float.MAX_EXPONENT, we round up (correctly) to
3654          * Float.POSITIVE_INFINITY.
3655          */
3656         bits |= signum & FloatConsts.SIGN_BIT_MASK;
3657         return Float.intBitsToFloat(bits);
3658     }
3659 
3660     /**
3661      * Converts this BigInteger to a {@code double}.  This
3662      * conversion is similar to the
3663      * <i>narrowing primitive conversion</i> from {@code double} to
3664      * {@code float} as defined in section 5.1.3 of
3665      * <cite>The Java&trade; Language Specification</cite>:
3666      * if this BigInteger has too great a magnitude
3667      * to represent as a {@code double}, it will be converted to
3668      * {@link Double#NEGATIVE_INFINITY} or {@link
3669      * Double#POSITIVE_INFINITY} as appropriate.  Note that even when
3670      * the return value is finite, this conversion can lose
3671      * information about the precision of the BigInteger value.
3672      *
3673      * @return this BigInteger converted to a {@code double}.
3674      */
3675     public double doubleValue() {
3676         if (signum == 0) {
3677             return 0.0;
3678         }
3679 
3680         int exponent = ((mag.length - 1) << 5) + bitLengthForInt(mag[0]) - 1;
3681 
3682         // exponent == floor(log2(abs(this))Double)
3683         if (exponent < Long.SIZE - 1) {
3684             return longValue();
3685         } else if (exponent > Double.MAX_EXPONENT) {
3686             return signum > 0 ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3687         }
3688 
3689         /*
3690          * We need the top SIGNIFICAND_WIDTH bits, including the "implicit"
3691          * one bit. To make rounding easier, we pick out the top
3692          * SIGNIFICAND_WIDTH + 1 bits, so we have one to help us round up or
3693          * down. twiceSignifFloor will contain the top SIGNIFICAND_WIDTH + 1
3694          * bits, and signifFloor the top SIGNIFICAND_WIDTH.
3695          *
3696          * It helps to consider the real number signif = abs(this) *
3697          * 2^(SIGNIFICAND_WIDTH - 1 - exponent).
3698          */
3699         int shift = exponent - DoubleConsts.SIGNIFICAND_WIDTH;
3700 
3701         long twiceSignifFloor;
3702         // twiceSignifFloor will be == abs().shiftRight(shift).longValue()
3703         // We do the shift into a long directly to improve performance.
3704 
3705         int nBits = shift & 0x1f;
3706         int nBits2 = 32 - nBits;
3707 
3708         int highBits;
3709         int lowBits;
3710         if (nBits == 0) {
3711             highBits = mag[0];
3712             lowBits = mag[1];
3713         } else {
3714             highBits = mag[0] >>> nBits;
3715             lowBits = (mag[0] << nBits2) | (mag[1] >>> nBits);
3716             if (highBits == 0) {
3717                 highBits = lowBits;
3718                 lowBits = (mag[1] << nBits2) | (mag[2] >>> nBits);
3719             }
3720         }
3721 
3722         twiceSignifFloor = ((highBits & LONG_MASK) << 32)
3723                 | (lowBits & LONG_MASK);
3724 
3725         long signifFloor = twiceSignifFloor >> 1;
3726         signifFloor &= DoubleConsts.SIGNIF_BIT_MASK; // remove the implied bit
3727 
3728         /*
3729          * We round up if either the fractional part of signif is strictly
3730          * greater than 0.5 (which is true if the 0.5 bit is set and any lower
3731          * bit is set), or if the fractional part of signif is >= 0.5 and
3732          * signifFloor is odd (which is true if both the 0.5 bit and the 1 bit
3733          * are set). This is equivalent to the desired HALF_EVEN rounding.
3734          */
3735         boolean increment = (twiceSignifFloor & 1) != 0
3736                 && ((signifFloor & 1) != 0 || abs().getLowestSetBit() < shift);
3737         long signifRounded = increment ? signifFloor + 1 : signifFloor;
3738         long bits = (long) ((exponent + DoubleConsts.EXP_BIAS))
3739                 << (DoubleConsts.SIGNIFICAND_WIDTH - 1);
3740         bits += signifRounded;
3741         /*
3742          * If signifRounded == 2^53, we'd need to set all of the significand
3743          * bits to zero and add 1 to the exponent. This is exactly the behavior
3744          * we get from just adding signifRounded to bits directly. If the
3745          * exponent is Double.MAX_EXPONENT, we round up (correctly) to
3746          * Double.POSITIVE_INFINITY.
3747          */
3748         bits |= signum & DoubleConsts.SIGN_BIT_MASK;
3749         return Double.longBitsToDouble(bits);
3750     }
3751 
3752     /**
3753      * Returns a copy of the input array stripped of any leading zero bytes.
3754      */
3755     private static int[] stripLeadingZeroInts(int val[]) {
3756         int vlen = val.length;
3757         int keep;
3758 
3759         // Find first nonzero byte
3760         for (keep = 0; keep < vlen && val[keep] == 0; keep++)
3761             ;
3762         return java.util.Arrays.copyOfRange(val, keep, vlen);
3763     }
3764 
3765     /**
3766      * Returns the input array stripped of any leading zero bytes.
3767      * Since the source is trusted the copying may be skipped.
3768      */
3769     private static int[] trustedStripLeadingZeroInts(int val[]) {
3770         int vlen = val.length;
3771         int keep;
3772 
3773         // Find first nonzero byte
3774         for (keep = 0; keep < vlen && val[keep] == 0; keep++)
3775             ;
3776         return keep == 0 ? val : java.util.Arrays.copyOfRange(val, keep, vlen);
3777     }
3778 
3779     /**
3780      * Returns a copy of the input array stripped of any leading zero bytes.
3781      */
3782     private static int[] stripLeadingZeroBytes(byte a[]) {
3783         int byteLength = a.length;
3784         int keep;
3785 
3786         // Find first nonzero byte
3787         for (keep = 0; keep < byteLength && a[keep]==0; keep++)
3788             ;
3789 
3790         // Allocate new array and copy relevant part of input array
3791         int intLength = ((byteLength - keep) + 3) >>> 2;
3792         int[] result = new int[intLength];
3793         int b = byteLength - 1;
3794         for (int i = intLength-1; i >= 0; i--) {
3795             result[i] = a[b--] & 0xff;
3796             int bytesRemaining = b - keep + 1;
3797             int bytesToTransfer = Math.min(3, bytesRemaining);
3798             for (int j=8; j <= (bytesToTransfer << 3); j += 8)
3799                 result[i] |= ((a[b--] & 0xff) << j);
3800         }
3801         return result;
3802     }
3803 
3804     /**
3805      * Takes an array a representing a negative 2's-complement number and
3806      * returns the minimal (no leading zero bytes) unsigned whose value is -a.
3807      */
3808     private static int[] makePositive(byte a[]) {
3809         int keep, k;
3810         int byteLength = a.length;
3811 
3812         // Find first non-sign (0xff) byte of input
3813         for (keep=0; keep<byteLength && a[keep]==-1; keep++)
3814             ;
3815 
3816 
3817         /* Allocate output array.  If all non-sign bytes are 0x00, we must
3818          * allocate space for one extra output byte. */
3819         for (k=keep; k<byteLength && a[k]==0; k++)
3820             ;
3821 
3822         int extraByte = (k==byteLength) ? 1 : 0;
3823         int intLength = ((byteLength - keep + extraByte) + 3)/4;
3824         int result[] = new int[intLength];
3825 
3826         /* Copy one's complement of input into output, leaving extra
3827          * byte (if it exists) == 0x00 */
3828         int b = byteLength - 1;
3829         for (int i = intLength-1; i >= 0; i--) {
3830             result[i] = a[b--] & 0xff;
3831             int numBytesToTransfer = Math.min(3, b-keep+1);
3832             if (numBytesToTransfer < 0)
3833                 numBytesToTransfer = 0;
3834             for (int j=8; j <= 8*numBytesToTransfer; j += 8)
3835                 result[i] |= ((a[b--] & 0xff) << j);
3836 
3837             // Mask indicates which bits must be complemented
3838             int mask = -1 >>> (8*(3-numBytesToTransfer));
3839             result[i] = ~result[i] & mask;
3840         }
3841 
3842         // Add one to one's complement to generate two's complement
3843         for (int i=result.length-1; i>=0; i--) {
3844             result[i] = (int)((result[i] & LONG_MASK) + 1);
3845             if (result[i] != 0)
3846                 break;
3847         }
3848 
3849         return result;
3850     }
3851 
3852     /**
3853      * Takes an array a representing a negative 2's-complement number and
3854      * returns the minimal (no leading zero ints) unsigned whose value is -a.
3855      */
3856     private static int[] makePositive(int a[]) {
3857         int keep, j;
3858 
3859         // Find first non-sign (0xffffffff) int of input
3860         for (keep=0; keep<a.length && a[keep]==-1; keep++)
3861             ;
3862 
3863         /* Allocate output array.  If all non-sign ints are 0x00, we must
3864          * allocate space for one extra output int. */
3865         for (j=keep; j<a.length && a[j]==0; j++)
3866             ;
3867         int extraInt = (j==a.length ? 1 : 0);
3868         int result[] = new int[a.length - keep + extraInt];
3869 
3870         /* Copy one's complement of input into output, leaving extra
3871          * int (if it exists) == 0x00 */
3872         for (int i = keep; i<a.length; i++)
3873             result[i - keep + extraInt] = ~a[i];
3874 
3875         // Add one to one's complement to generate two's complement
3876         for (int i=result.length-1; ++result[i]==0; i--)
3877             ;
3878 
3879         return result;
3880     }
3881 
3882     /*
3883      * The following two arrays are used for fast String conversions.  Both
3884      * are indexed by radix.  The first is the number of digits of the given
3885      * radix that can fit in a Java long without "going negative", i.e., the
3886      * highest integer n such that radix**n < 2**63.  The second is the
3887      * "long radix" that tears each number into "long digits", each of which
3888      * consists of the number of digits in the corresponding element in
3889      * digitsPerLong (longRadix[i] = i**digitPerLong[i]).  Both arrays have
3890      * nonsense values in their 0 and 1 elements, as radixes 0 and 1 are not
3891      * used.
3892      */
3893     private static int digitsPerLong[] = {0, 0,
3894         62, 39, 31, 27, 24, 22, 20, 19, 18, 18, 17, 17, 16, 16, 15, 15, 15, 14,
3895         14, 14, 14, 13, 13, 13, 13, 13, 13, 12, 12, 12, 12, 12, 12, 12, 12};
3896 
3897     private static BigInteger longRadix[] = {null, null,
3898         valueOf(0x4000000000000000L), valueOf(0x383d9170b85ff80bL),
3899         valueOf(0x4000000000000000L), valueOf(0x6765c793fa10079dL),
3900         valueOf(0x41c21cb8e1000000L), valueOf(0x3642798750226111L),
3901         valueOf(0x1000000000000000L), valueOf(0x12bf307ae81ffd59L),
3902         valueOf( 0xde0b6b3a7640000L), valueOf(0x4d28cb56c33fa539L),
3903         valueOf(0x1eca170c00000000L), valueOf(0x780c7372621bd74dL),
3904         valueOf(0x1e39a5057d810000L), valueOf(0x5b27ac993df97701L),
3905         valueOf(0x1000000000000000L), valueOf(0x27b95e997e21d9f1L),
3906         valueOf(0x5da0e1e53c5c8000L), valueOf( 0xb16a458ef403f19L),
3907         valueOf(0x16bcc41e90000000L), valueOf(0x2d04b7fdd9c0ef49L),
3908         valueOf(0x5658597bcaa24000L), valueOf( 0x6feb266931a75b7L),
3909         valueOf( 0xc29e98000000000L), valueOf(0x14adf4b7320334b9L),
3910         valueOf(0x226ed36478bfa000L), valueOf(0x383d9170b85ff80bL),
3911         valueOf(0x5a3c23e39c000000L), valueOf( 0x4e900abb53e6b71L),
3912         valueOf( 0x7600ec618141000L), valueOf( 0xaee5720ee830681L),
3913         valueOf(0x1000000000000000L), valueOf(0x172588ad4f5f0981L),
3914         valueOf(0x211e44f7d02c1000L), valueOf(0x2ee56725f06e5c71L),
3915         valueOf(0x41c21cb8e1000000L)};
3916 
3917     /*
3918      * These two arrays are the integer analogue of above.
3919      */
3920     private static int digitsPerInt[] = {0, 0, 30, 19, 15, 13, 11,
3921         11, 10, 9, 9, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6,
3922         6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5};
3923 
3924     private static int intRadix[] = {0, 0,
3925         0x40000000, 0x4546b3db, 0x40000000, 0x48c27395, 0x159fd800,
3926         0x75db9c97, 0x40000000, 0x17179149, 0x3b9aca00, 0xcc6db61,
3927         0x19a10000, 0x309f1021, 0x57f6c100, 0xa2f1b6f,  0x10000000,
3928         0x18754571, 0x247dbc80, 0x3547667b, 0x4c4b4000, 0x6b5a6e1d,
3929         0x6c20a40,  0x8d2d931,  0xb640000,  0xe8d4a51,  0x1269ae40,
3930         0x17179149, 0x1cb91000, 0x23744899, 0x2b73a840, 0x34e63b41,
3931         0x40000000, 0x4cfa3cc1, 0x5c13d840, 0x6d91b519, 0x39aa400
3932     };
3933 
3934     /**
3935      * These routines provide access to the two's complement representation
3936      * of BigIntegers.
3937      */
3938 
3939     /**
3940      * Returns the length of the two's complement representation in ints,
3941      * including space for at least one sign bit.
3942      */
3943     private int intLength() {
3944         return (bitLength() >>> 5) + 1;
3945     }
3946 
3947     /* Returns sign bit */
3948     private int signBit() {
3949         return signum < 0 ? 1 : 0;
3950     }
3951 
3952     /* Returns an int of sign bits */
3953     private int signInt() {
3954         return signum < 0 ? -1 : 0;
3955     }
3956 
3957     /**
3958      * Returns the specified int of the little-endian two's complement
3959      * representation (int 0 is the least significant).  The int number can
3960      * be arbitrarily high (values are logically preceded by infinitely many
3961      * sign ints).
3962      */
3963     private int getInt(int n) {
3964         if (n < 0)
3965             return 0;
3966         if (n >= mag.length)
3967             return signInt();
3968 
3969         int magInt = mag[mag.length-n-1];
3970 
3971         return (signum >= 0 ? magInt :
3972                 (n <= firstNonzeroIntNum() ? -magInt : ~magInt));
3973     }
3974 
3975     /**
3976      * Returns the index of the int that contains the first nonzero int in the
3977      * little-endian binary representation of the magnitude (int 0 is the
3978      * least significant). If the magnitude is zero, return value is undefined.
3979      */
3980     private int firstNonzeroIntNum() {
3981         int fn = firstNonzeroIntNum - 2;
3982         if (fn == -2) { // firstNonzeroIntNum not initialized yet
3983             fn = 0;
3984 
3985             // Search for the first nonzero int
3986             int i;
3987             int mlen = mag.length;
3988             for (i = mlen - 1; i >= 0 && mag[i] == 0; i--)
3989                 ;
3990             fn = mlen - i - 1;
3991             firstNonzeroIntNum = fn + 2; // offset by two to initialize
3992         }
3993         return fn;
3994     }
3995 
3996     /** use serialVersionUID from JDK 1.1. for interoperability */
3997     private static final long serialVersionUID = -8287574255936472291L;
3998 
3999     /**
4000      * Serializable fields for BigInteger.
4001      *
4002      * @serialField signum  int
4003      *              signum of this BigInteger.
4004      * @serialField magnitude int[]
4005      *              magnitude array of this BigInteger.
4006      * @serialField bitCount  int
4007      *              number of bits in this BigInteger
4008      * @serialField bitLength int
4009      *              the number of bits in the minimal two's-complement
4010      *              representation of this BigInteger
4011      * @serialField lowestSetBit int
4012      *              lowest set bit in the twos complement representation
4013      */
4014     private static final ObjectStreamField[] serialPersistentFields = {
4015         new ObjectStreamField("signum", Integer.TYPE),
4016         new ObjectStreamField("magnitude", byte[].class),
4017         new ObjectStreamField("bitCount", Integer.TYPE),
4018         new ObjectStreamField("bitLength", Integer.TYPE),
4019         new ObjectStreamField("firstNonzeroByteNum", Integer.TYPE),
4020         new ObjectStreamField("lowestSetBit", Integer.TYPE)
4021         };
4022 
4023     /**
4024      * Reconstitute the {@code BigInteger} instance from a stream (that is,
4025      * deserialize it). The magnitude is read in as an array of bytes
4026      * for historical reasons, but it is converted to an array of ints
4027      * and the byte array is discarded.
4028      * Note:
4029      * The current convention is to initialize the cache fields, bitCount,
4030      * bitLength and lowestSetBit, to 0 rather than some other marker value.
4031      * Therefore, no explicit action to set these fields needs to be taken in
4032      * readObject because those fields already have a 0 value be default since
4033      * defaultReadObject is not being used.
4034      */
4035     private void readObject(java.io.ObjectInputStream s)
4036         throws java.io.IOException, ClassNotFoundException {
4037         /*
4038          * In order to maintain compatibility with previous serialized forms,
4039          * the magnitude of a BigInteger is serialized as an array of bytes.
4040          * The magnitude field is used as a temporary store for the byte array
4041          * that is deserialized. The cached computation fields should be
4042          * transient but are serialized for compatibility reasons.
4043          */
4044 
4045         // prepare to read the alternate persistent fields
4046         ObjectInputStream.GetField fields = s.readFields();
4047 
4048         // Read the alternate persistent fields that we care about
4049         int sign = fields.get("signum", -2);
4050         byte[] magnitude = (byte[])fields.get("magnitude", null);
4051 
4052         // Validate signum
4053         if (sign < -1 || sign > 1) {
4054             String message = "BigInteger: Invalid signum value";
4055             if (fields.defaulted("signum"))
4056                 message = "BigInteger: Signum not present in stream";
4057             throw new java.io.StreamCorruptedException(message);
4058         }
4059         if ((magnitude.length == 0) != (sign == 0)) {
4060             String message = "BigInteger: signum-magnitude mismatch";
4061             if (fields.defaulted("magnitude"))
4062                 message = "BigInteger: Magnitude not present in stream";
4063             throw new java.io.StreamCorruptedException(message);
4064         }
4065 
4066         // Commit final fields via Unsafe
4067         UnsafeHolder.putSign(this, sign);
4068 
4069         // Calculate mag field from magnitude and discard magnitude
4070         UnsafeHolder.putMag(this, stripLeadingZeroBytes(magnitude));
4071     }
4072 
4073     // Support for resetting final fields while deserializing
4074     private static class UnsafeHolder {
4075         private static final sun.misc.Unsafe unsafe;
4076         private static final long signumOffset;
4077         private static final long magOffset;
4078         static {
4079             try {
4080                 unsafe = sun.misc.Unsafe.getUnsafe();
4081                 signumOffset = unsafe.objectFieldOffset
4082                     (BigInteger.class.getDeclaredField("signum"));
4083                 magOffset = unsafe.objectFieldOffset
4084                     (BigInteger.class.getDeclaredField("mag"));
4085             } catch (Exception ex) {
4086                 throw new ExceptionInInitializerError(ex);
4087             }
4088         }
4089 
4090         static void putSign(BigInteger bi, int sign) {
4091             unsafe.putIntVolatile(bi, signumOffset, sign);
4092         }
4093 
4094         static void putMag(BigInteger bi, int[] magnitude) {
4095             unsafe.putObjectVolatile(bi, magOffset, magnitude);
4096         }
4097     }
4098 
4099     /**
4100      * Save the {@code BigInteger} instance to a stream.
4101      * The magnitude of a BigInteger is serialized as a byte array for
4102      * historical reasons.
4103      *
4104      * @serialData two necessary fields are written as well as obsolete
4105      *             fields for compatibility with older versions.
4106      */
4107     private void writeObject(ObjectOutputStream s) throws IOException {
4108         // set the values of the Serializable fields
4109         ObjectOutputStream.PutField fields = s.putFields();
4110         fields.put("signum", signum);
4111         fields.put("magnitude", magSerializedForm());
4112         // The values written for cached fields are compatible with older
4113         // versions, but are ignored in readObject so don't otherwise matter.
4114         fields.put("bitCount", -1);
4115         fields.put("bitLength", -1);
4116         fields.put("lowestSetBit", -2);
4117         fields.put("firstNonzeroByteNum", -2);
4118 
4119         // save them
4120         s.writeFields();
4121 }
4122 
4123     /**
4124      * Returns the mag array as an array of bytes.
4125      */
4126     private byte[] magSerializedForm() {
4127         int len = mag.length;
4128 
4129         int bitLen = (len == 0 ? 0 : ((len - 1) << 5) + bitLengthForInt(mag[0]));
4130         int byteLen = (bitLen + 7) >>> 3;
4131         byte[] result = new byte[byteLen];
4132 
4133         for (int i = byteLen - 1, bytesCopied = 4, intIndex = len - 1, nextInt = 0;
4134              i>=0; i--) {
4135             if (bytesCopied == 4) {
4136                 nextInt = mag[intIndex--];
4137                 bytesCopied = 1;
4138             } else {
4139                 nextInt >>>= 8;
4140                 bytesCopied++;
4141             }
4142             result[i] = (byte)nextInt;
4143         }
4144         return result;
4145     }
4146 
4147     /**
4148      * Converts this {@code BigInteger} to a {@code long}, checking
4149      * for lost information.  If the value of this {@code BigInteger}
4150      * is out of the range of the {@code long} type, then an
4151      * {@code ArithmeticException} is thrown.
4152      *
4153      * @return this {@code BigInteger} converted to a {@code long}.
4154      * @throws ArithmeticException if the value of {@code this} will
4155      * not exactly fit in a {@code long}.
4156      * @see BigInteger#longValue
4157      * @since  1.8
4158      */
4159     public long longValueExact() {
4160         if (mag.length <= 2 && bitLength() <= 63)
4161             return longValue();
4162         else
4163             throw new ArithmeticException("BigInteger out of long range");
4164     }
4165 
4166     /**
4167      * Converts this {@code BigInteger} to an {@code int}, checking
4168      * for lost information.  If the value of this {@code BigInteger}
4169      * is out of the range of the {@code int} type, then an
4170      * {@code ArithmeticException} is thrown.
4171      *
4172      * @return this {@code BigInteger} converted to an {@code int}.
4173      * @throws ArithmeticException if the value of {@code this} will
4174      * not exactly fit in a {@code int}.
4175      * @see BigInteger#intValue
4176      * @since  1.8
4177      */
4178     public int intValueExact() {
4179         if (mag.length <= 1 && bitLength() <= 31)
4180             return intValue();
4181         else
4182             throw new ArithmeticException("BigInteger out of int range");
4183     }
4184 
4185     /**
4186      * Converts this {@code BigInteger} to a {@code short}, checking
4187      * for lost information.  If the value of this {@code BigInteger}
4188      * is out of the range of the {@code short} type, then an
4189      * {@code ArithmeticException} is thrown.
4190      *
4191      * @return this {@code BigInteger} converted to a {@code short}.
4192      * @throws ArithmeticException if the value of {@code this} will
4193      * not exactly fit in a {@code short}.
4194      * @see BigInteger#shortValue
4195      * @since  1.8
4196      */
4197     public short shortValueExact() {
4198         if (mag.length <= 1 && bitLength() <= 31) {
4199             int value = intValue();
4200             if (value >= Short.MIN_VALUE && value <= Short.MAX_VALUE)
4201                 return shortValue();
4202         }
4203         throw new ArithmeticException("BigInteger out of short range");
4204     }
4205 
4206     /**
4207      * Converts this {@code BigInteger} to a {@code byte}, checking
4208      * for lost information.  If the value of this {@code BigInteger}
4209      * is out of the range of the {@code byte} type, then an
4210      * {@code ArithmeticException} is thrown.
4211      *
4212      * @return this {@code BigInteger} converted to a {@code byte}.
4213      * @throws ArithmeticException if the value of {@code this} will
4214      * not exactly fit in a {@code byte}.
4215      * @see BigInteger#byteValue
4216      * @since  1.8
4217      */
4218     public byte byteValueExact() {
4219         if (mag.length <= 1 && bitLength() <= 31) {
4220             int value = intValue();
4221             if (value >= Byte.MIN_VALUE && value <= Byte.MAX_VALUE)
4222                 return byteValue();
4223         }
4224         throw new ArithmeticException("BigInteger out of byte range");
4225     }
4226 }
--- EOF ---