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





3468         }
3469         else
3470             retVal = cacheLine.get(exponent);
3471 
3472         return retVal;






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