1 /*
2  * Copyright (c) 1996, 2013, Oracle and/or its affiliates. All rights reserved.
3  * ORACLE PROPRIETARY/CONFIDENTIAL. Use is subject to license terms.
4  *
5  *
6  *
7  *
8  *
9  *
10  *
11  *
12  *
13  *
14  *
15  *
16  *
17  *
18  *
19  *
20  *
21  *
22  *
23  *
24  */

25
26 /*
27  * Portions Copyright (c) 1995  Colin Plumb.  All rights reserved.
28  */

29
30 package java.math;
31
32 import java.io.IOException;
33 import java.io.ObjectInputStream;
34 import java.io.ObjectOutputStream;
35 import java.io.ObjectStreamField;
36 import java.util.Arrays;
37 import java.util.Random;
38 import java.util.concurrent.ThreadLocalRandom;
39 import sun.misc.DoubleConsts;
40 import sun.misc.FloatConsts;
41
42 /**
43  * Immutable arbitrary-precision integers.  All operations behave as if
44  * BigIntegers were represented in two's-complement notation (like Java's
45  * primitive integer types).  BigInteger provides analogues to all of Java's
46  * primitive integer operators, and all relevant methods from java.lang.Math.
47  * Additionally, BigInteger provides operations for modular arithmetic, GCD
48  * calculation, primality testing, prime generation, bit manipulation,
49  * and a few other miscellaneous operations.
50  *
51  * <p>Semantics of arithmetic operations exactly mimic those of Java's integer
52  * arithmetic operators, as defined in <i>The Java Language Specification</i>.
53  * For example, division by zero throws an {@code ArithmeticException}, and
54  * division of a negative by a positive yields a negative (or zero) remainder.
55  * All of the details in the Spec concerning overflow are ignored, as
56  * BigIntegers are made as large as necessary to accommodate the results of an
57  * operation.
58  *
59  * <p>Semantics of shift operations extend those of Java's shift operators
60  * to allow for negative shift distances.  A right-shift with a negative
61  * shift distance results in a left shift, and vice-versa.  The unsigned
62  * right shift operator ({@code >>>}) is omitted, as this operation makes
63  * little sense in combination with the "infinite word size" abstraction
64  * provided by this class.
65  *
66  * <p>Semantics of bitwise logical operations exactly mimic those of Java's
67  * bitwise integer operators.  The binary operators ({@code and},
68  * {@code or}, {@code xor}) implicitly perform sign extension on the shorter
69  * of the two operands prior to performing the operation.
70  *
71  * <p>Comparison operations perform signed integer comparisons, analogous to
72  * those performed by Java's relational and equality operators.
73  *
74  * <p>Modular arithmetic operations are provided to compute residues, perform
75  * exponentiation, and compute multiplicative inverses.  These methods always
76  * return a non-negative result, between {@code 0} and {@code (modulus - 1)},
77  * inclusive.
78  *
79  * <p>Bit operations operate on a single bit of the two's-complement
80  * representation of their operand.  If necessary, the operand is sign-
81  * extended so that it contains the designated bit.  None of the single-bit
82  * operations can produce a BigInteger with a different sign from the
83  * BigInteger being operated on, as they affect only a single bit, and the
84  * "infinite word size" abstraction provided by this class ensures that there
85  * are infinitely many "virtual sign bits" preceding each BigInteger.
86  *
87  * <p>For the sake of brevity and clarity, pseudo-code is used throughout the
88  * descriptions of BigInteger methods.  The pseudo-code expression
89  * {@code (i + j)} is shorthand for "a BigInteger whose value is
90  * that of the BigInteger {@code i} plus that of the BigInteger {@code j}."
91  * The pseudo-code expression {@code (i == j)} is shorthand for
92  * "{@code trueif and only if the BigInteger {@code i} represents the same
93  * value as the BigInteger {@code j}."  Other pseudo-code expressions are
94  * interpreted similarly.
95  *
96  * <p>All methods and constructors in this class throw
97  * {@code NullPointerException} when passed
98  * a null object reference for any input parameter.
99  *
100  * BigInteger must support values in the range
101  * -2<sup>{@code Integer.MAX_VALUE}</sup> (exclusive) to
102  * +2<sup>{@code Integer.MAX_VALUE}</sup> (exclusive)
103  * and may support values outside of that range.
104  *
105  * The range of probable prime values is limited and may be less than
106  * the full supported positive range of {@code BigInteger}.
107  * The range must be at least 1 to 2<sup>500000000</sup>.
108  *
109  * @implNote
110  * BigInteger constructors and operations throw {@code ArithmeticException} when
111  * the result is out of the supported range of
112  * -2<sup>{@code Integer.MAX_VALUE}</sup> (exclusive) to
113  * +2<sup>{@code Integer.MAX_VALUE}</sup> (exclusive).
114  *
115  * @see     BigDecimal
116  * @author  Josh Bloch
117  * @author  Michael McCloskey
118  * @author  Alan Eliasen
119  * @author  Timothy Buktu
120  * @since JDK1.1
121  */

122
123 public class BigInteger extends Number implements Comparable<BigInteger> {
124     /**
125      * The signum of this BigInteger: -1 for negative, 0 for zero, or
126      * 1 for positive.  Note that the BigInteger zero <i>must</i> have
127      * a signum of 0.  This is necessary to ensures that there is exactly one
128      * representation for each BigInteger value.
129      *
130      * @serial
131      */

132     final int signum;
133
134     /**
135      * The magnitude of this BigInteger, in <i>big-endian</i> order: the
136      * zeroth element of this array is the most-significant int of the
137      * magnitude.  The magnitude must be "minimal" in that the most-significant
138      * int ({@code mag[0]}) must be non-zero.  This is necessary to
139      * ensure that there is exactly one representation for each BigInteger
140      * value.  Note that this implies that the BigInteger zero has a
141      * zero-length mag array.
142      */

143     final int[] mag;
144
145     // These "redundant fields" are initialized with recognizable nonsense
146     // values, and cached the first time they are needed (or never, if they
147     // aren't needed).
148
149      /**
150      * One plus the bitCount of this BigInteger. Zeros means unitialized.
151      *
152      * @serial
153      * @see #bitCount
154      * @deprecated Deprecated since logical value is offset from stored
155      * value and correction factor is applied in accessor method.
156      */

157     @Deprecated
158     private int bitCount;
159
160     /**
161      * One plus the bitLength of this BigInteger. Zeros means unitialized.
162      * (either value is acceptable).
163      *
164      * @serial
165      * @see #bitLength()
166      * @deprecated Deprecated since logical value is offset from stored
167      * value and correction factor is applied in accessor method.
168      */

169     @Deprecated
170     private int bitLength;
171
172     /**
173      * Two plus the lowest set bit of this BigInteger, as returned by
174      * getLowestSetBit().
175      *
176      * @serial
177      * @see #getLowestSetBit
178      * @deprecated Deprecated since logical value is offset from stored
179      * value and correction factor is applied in accessor method.
180      */

181     @Deprecated
182     private int lowestSetBit;
183
184     /**
185      * Two plus the index of the lowest-order int in the magnitude of this
186      * BigInteger that contains a nonzero int, or -2 (either value is acceptable).
187      * The least significant int has int-number 0, the next int in order of
188      * increasing significance has int-number 1, and so forth.
189      * @deprecated Deprecated since logical value is offset from stored
190      * value and correction factor is applied in accessor method.
191      */

192     @Deprecated
193     private int firstNonzeroIntNum;
194
195     /**
196      * This mask is used to obtain the value of an int as if it were unsigned.
197      */

198     final static long LONG_MASK = 0xffffffffL;
199
200     /**
201      * This constant limits {@code mag.length} of BigIntegers to the supported
202      * range.
203      */

204     private static final int MAX_MAG_LENGTH = Integer.MAX_VALUE / Integer.SIZE + 1; // (1 << 26)
205
206     /**
207      * Bit lengths larger than this constant can cause overflow in searchLen
208      * calculation and in BitSieve.singleSearch method.
209      */

210     private static final  int PRIME_SEARCH_BIT_LENGTH_LIMIT = 500000000;
211
212     /**
213      * The threshold value for using Karatsuba multiplication.  If the number
214      * of ints in both mag arrays are greater than this number, then
215      * Karatsuba multiplication will be used.   This value is found
216      * experimentally to work well.
217      */

218     private static final int KARATSUBA_THRESHOLD = 80;
219
220     /**
221      * The threshold value for using 3-way Toom-Cook multiplication.
222      * If the number of ints in each mag array is greater than the
223      * Karatsuba threshold, and the number of ints in at least one of
224      * the mag arrays is greater than this threshold, then Toom-Cook
225      * multiplication will be used.
226      */

227     private static final int TOOM_COOK_THRESHOLD = 240;
228
229     /**
230      * The threshold value for using Karatsuba squaring.  If the number
231      * of ints in the number are larger than this value,
232      * Karatsuba squaring will be used.   This value is found
233      * experimentally to work well.
234      */

235     private static final int KARATSUBA_SQUARE_THRESHOLD = 128;
236
237     /**
238      * The threshold value for using Toom-Cook squaring.  If the number
239      * of ints in the number are larger than this value,
240      * Toom-Cook squaring will be used.   This value is found
241      * experimentally to work well.
242      */

243     private static final int TOOM_COOK_SQUARE_THRESHOLD = 216;
244
245     /**
246      * The threshold value for using Burnikel-Ziegler division.  If the number
247      * of ints in the divisor are larger than this value, Burnikel-Ziegler
248      * division may be used.  This value is found experimentally to work well.
249      */

250     static final int BURNIKEL_ZIEGLER_THRESHOLD = 80;
251
252     /**
253      * The offset value for using Burnikel-Ziegler division.  If the number
254      * of ints in the divisor exceeds the Burnikel-Ziegler threshold, and the
255      * number of ints in the dividend is greater than the number of ints in the
256      * divisor plus this value, Burnikel-Ziegler division will be used.  This
257      * value is found experimentally to work well.
258      */

259     static final int BURNIKEL_ZIEGLER_OFFSET = 40;
260
261     /**
262      * The threshold value for using Schoenhage recursive base conversion. If
263      * the number of ints in the number are larger than this value,
264      * the Schoenhage algorithm will be used.  In practice, it appears that the
265      * Schoenhage routine is faster for any threshold down to 2, and is
266      * relatively flat for thresholds between 2-25, so this choice may be
267      * varied within this range for very small effect.
268      */

269     private static final int SCHOENHAGE_BASE_CONVERSION_THRESHOLD = 20;
270
271     /**
272      * The threshold value for using squaring code to perform multiplication
273      * of a {@code BigInteger} instance by itself.  If the number of ints in
274      * the number are larger than this value, {@code multiply(this)} will
275      * return {@code square()}.
276      */

277     private static final int MULTIPLY_SQUARE_THRESHOLD = 20;
278
279     /**
280      * The threshold for using an intrinsic version of
281      * implMontgomeryXXX to perform Montgomery multiplication.  If the
282      * number of ints in the number is more than this value we do not
283      * use the intrinsic.
284      */

285     private static final int MONTGOMERY_INTRINSIC_THRESHOLD = 512;
286
287
288     // Constructors
289
290     /**
291      * Translates a byte array containing the two's-complement binary
292      * representation of a BigInteger into a BigInteger.  The input array is
293      * assumed to be in <i>big-endian</i> byte-order: the most significant
294      * byte is in the zeroth element.
295      *
296      * @param  val big-endian two's-complement binary representation of
297      *         BigInteger.
298      * @throws NumberFormatException {@code val} is zero bytes long.
299      */

300     public BigInteger(byte[] val) {
301         if (val.length == 0)
302             throw new NumberFormatException("Zero length BigInteger");
303
304         if (val[0] < 0) {
305             mag = makePositive(val);
306             signum = -1;
307         } else {
308             mag = stripLeadingZeroBytes(val);
309             signum = (mag.length == 0 ? 0 : 1);
310         }
311         if (mag.length >= MAX_MAG_LENGTH) {
312             checkRange();
313         }
314     }
315
316     /**
317      * This private constructor translates an int array containing the
318      * two's-complement binary representation of a BigInteger into a
319      * BigInteger. The input array is assumed to be in <i>big-endian</i>
320      * int-order: the most significant int is in the zeroth element.
321      */

322     private BigInteger(int[] val) {
323         if (val.length == 0)
324             throw new NumberFormatException("Zero length BigInteger");
325
326         if (val[0] < 0) {
327             mag = makePositive(val);
328             signum = -1;
329         } else {
330             mag = trustedStripLeadingZeroInts(val);
331             signum = (mag.length == 0 ? 0 : 1);
332         }
333         if (mag.length >= MAX_MAG_LENGTH) {
334             checkRange();
335         }
336     }
337
338     /**
339      * Translates the sign-magnitude representation of a BigInteger into a
340      * BigInteger.  The sign is represented as an integer signum value: -1 for
341      * negative, 0 for zero, or 1 for positive.  The magnitude is a byte array
342      * in <i>big-endian</i> byte-order: the most significant byte is in the
343      * zeroth element.  A zero-length magnitude array is permissible, and will
344      * result in a BigInteger value of 0, whether signum is -1, 0 or 1.
345      *
346      * @param  signum signum of the number (-1 for negative, 0 for zero, 1
347      *         for positive).
348      * @param  magnitude big-endian binary representation of the magnitude of
349      *         the number.
350      * @throws NumberFormatException {@code signum} is not one of the three
351      *         legal values (-1, 0, and 1), or {@code signum} is 0 and
352      *         {@code magnitude} contains one or more non-zero bytes.
353      */

354     public BigInteger(int signum, byte[] magnitude) {
355         this.mag = stripLeadingZeroBytes(magnitude);
356
357         if (signum < -1 || signum > 1)
358             throw(new NumberFormatException("Invalid signum value"));
359
360         if (this.mag.length == 0) {
361             this.signum = 0;
362         } else {
363             if (signum == 0)
364                 throw(new NumberFormatException("signum-magnitude mismatch"));
365             this.signum = signum;
366         }
367         if (mag.length >= MAX_MAG_LENGTH) {
368             checkRange();
369         }
370     }
371
372     /**
373      * A constructor for internal use that translates the sign-magnitude
374      * representation of a BigInteger into a BigInteger. It checks the
375      * arguments and copies the magnitude so this constructor would be
376      * safe for external use.
377      */

378     private BigInteger(int signum, int[] magnitude) {
379         this.mag = stripLeadingZeroInts(magnitude);
380
381         if (signum < -1 || signum > 1)
382             throw(new NumberFormatException("Invalid signum value"));
383
384         if (this.mag.length == 0) {
385             this.signum = 0;
386         } else {
387             if (signum == 0)
388                 throw(new NumberFormatException("signum-magnitude mismatch"));
389             this.signum = signum;
390         }
391         if (mag.length >= MAX_MAG_LENGTH) {
392             checkRange();
393         }
394     }
395
396     /**
397      * Translates the String representation of a BigInteger in the
398      * specified radix into a BigInteger.  The String representation
399      * consists of an optional minus or plus sign followed by a
400      * sequence of one or more digits in the specified radix.  The
401      * character-to-digit mapping is provided by {@code
402      * Character.digit}.  The String may not contain any extraneous
403      * characters (whitespace, for example).
404      *
405      * @param val String representation of BigInteger.
406      * @param radix radix to be used in interpreting {@code val}.
407      * @throws NumberFormatException {@code val} is not a valid representation
408      *         of a BigInteger in the specified radix, or {@code radix} is
409      *         outside the range from {@link Character#MIN_RADIX} to
410      *         {@link Character#MAX_RADIX}, inclusive.
411      * @see    Character#digit
412      */

413     public BigInteger(String val, int radix) {
414         int cursor = 0, numDigits;
415         final int len = val.length();
416
417         if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
418             throw new NumberFormatException("Radix out of range");
419         if (len == 0)
420             throw new NumberFormatException("Zero length BigInteger");
421
422         // Check for at most one leading sign
423         int sign = 1;
424         int index1 = val.lastIndexOf('-');
425         int index2 = val.lastIndexOf('+');
426         if (index1 >= 0) {
427             if (index1 != 0 || index2 >= 0) {
428                 throw new NumberFormatException("Illegal embedded sign character");
429             }
430             sign = -1;
431             cursor = 1;
432         } else if (index2 >= 0) {
433             if (index2 != 0) {
434                 throw new NumberFormatException("Illegal embedded sign character");
435             }
436             cursor = 1;
437         }
438         if (cursor == len)
439             throw new NumberFormatException("Zero length BigInteger");
440
441         // Skip leading zeros and compute number of digits in magnitude
442         while (cursor < len &&
443                Character.digit(val.charAt(cursor), radix) == 0) {
444             cursor++;
445         }
446
447         if (cursor == len) {
448             signum = 0;
449             mag = ZERO.mag;
450             return;
451         }
452
453         numDigits = len - cursor;
454         signum = sign;
455
456         // Pre-allocate array of expected size. May be too large but can
457         // never be too small. Typically exact.
458         long numBits = ((numDigits * bitsPerDigit[radix]) >>> 10) + 1;
459         if (numBits + 31 >= (1L << 32)) {
460             reportOverflow();
461         }
462         int numWords = (int) (numBits + 31) >>> 5;
463         int[] magnitude = new int[numWords];
464
465         // Process first (potentially short) digit group
466         int firstGroupLen = numDigits % digitsPerInt[radix];
467         if (firstGroupLen == 0)
468             firstGroupLen = digitsPerInt[radix];
469         String group = val.substring(cursor, cursor += firstGroupLen);
470         magnitude[numWords - 1] = Integer.parseInt(group, radix);
471         if (magnitude[numWords - 1] < 0)
472             throw new NumberFormatException("Illegal digit");
473
474         // Process remaining digit groups
475         int superRadix = intRadix[radix];
476         int groupVal = 0;
477         while (cursor < len) {
478             group = val.substring(cursor, cursor += digitsPerInt[radix]);
479             groupVal = Integer.parseInt(group, radix);
480             if (groupVal < 0)
481                 throw new NumberFormatException("Illegal digit");
482             destructiveMulAdd(magnitude, superRadix, groupVal);
483         }
484         // Required for cases where the array was overallocated.
485         mag = trustedStripLeadingZeroInts(magnitude);
486         if (mag.length >= MAX_MAG_LENGTH) {
487             checkRange();
488         }
489     }
490
491     /*
492      * Constructs a new BigInteger using a char array with radix=10.
493      * Sign is precalculated outside and not allowed in the val.
494      */

495     BigInteger(char[] val, int sign, int len) {
496         int cursor = 0, numDigits;
497
498         // Skip leading zeros and compute number of digits in magnitude
499         while (cursor < len && Character.digit(val[cursor], 10) == 0) {
500             cursor++;
501         }
502         if (cursor == len) {
503             signum = 0;
504             mag = ZERO.mag;
505             return;
506         }
507
508         numDigits = len - cursor;
509         signum = sign;
510         // Pre-allocate array of expected size
511         int numWords;
512         if (len < 10) {
513             numWords = 1;
514         } else {
515             long numBits = ((numDigits * bitsPerDigit[10]) >>> 10) + 1;
516             if (numBits + 31 >= (1L << 32)) {
517                 reportOverflow();
518             }
519             numWords = (int) (numBits + 31) >>> 5;
520         }
521         int[] magnitude = new int[numWords];
522
523         // Process first (potentially short) digit group
524         int firstGroupLen = numDigits % digitsPerInt[10];
525         if (firstGroupLen == 0)
526             firstGroupLen = digitsPerInt[10];
527         magnitude[numWords - 1] = parseInt(val, cursor,  cursor += firstGroupLen);
528
529         // Process remaining digit groups
530         while (cursor < len) {
531             int groupVal = parseInt(val, cursor, cursor += digitsPerInt[10]);
532             destructiveMulAdd(magnitude, intRadix[10], groupVal);
533         }
534         mag = trustedStripLeadingZeroInts(magnitude);
535         if (mag.length >= MAX_MAG_LENGTH) {
536             checkRange();
537         }
538     }
539
540     // Create an integer with the digits between the two indexes
541     // Assumes start < end. The result may be negative, but it
542     // is to be treated as an unsigned value.
543     private int parseInt(char[] source, int start, int end) {
544         int result = Character.digit(source[start++], 10);
545         if (result == -1)
546             throw new NumberFormatException(new String(source));
547
548         for (int index = start; index < end; index++) {
549             int nextVal = Character.digit(source[index], 10);
550             if (nextVal == -1)
551                 throw new NumberFormatException(new String(source));
552             result = 10*result + nextVal;
553         }
554
555         return result;
556     }
557
558     // bitsPerDigit in the given radix times 1024
559     // Rounded up to avoid underallocation.
560     private static long bitsPerDigit[] = { 0, 0,
561         1024, 1624, 2048, 2378, 2648, 2875, 3072, 3247, 3402, 3543, 3672,
562         3790, 3899, 4001, 4096, 4186, 4271, 4350, 4426, 4498, 4567, 4633,
563         4696, 4756, 4814, 4870, 4923, 4975, 5025, 5074, 5120, 5166, 5210,
564                                            5253, 5295};
565
566     // Multiply x array times word y in place, and add word z
567     private static void destructiveMulAdd(int[] x, int y, int z) {
568         // Perform the multiplication word by word
569         long ylong = y & LONG_MASK;
570         long zlong = z & LONG_MASK;
571         int len = x.length;
572
573         long product = 0;
574         long carry = 0;
575         for (int i = len-1; i >= 0; i--) {
576             product = ylong * (x[i] & LONG_MASK) + carry;
577             x[i] = (int)product;
578             carry = product >>> 32;
579         }
580
581         // Perform the addition
582         long sum = (x[len-1] & LONG_MASK) + zlong;
583         x[len-1] = (int)sum;
584         carry = sum >>> 32;
585         for (int i = len-2; i >= 0; i--) {
586             sum = (x[i] & LONG_MASK) + carry;
587             x[i] = (int)sum;
588             carry = sum >>> 32;
589         }
590     }
591
592     /**
593      * Translates the decimal String representation of a BigInteger into a
594      * BigInteger.  The String representation consists of an optional minus
595      * sign followed by a sequence of one or more decimal digits.  The
596      * character-to-digit mapping is provided by {@code Character.digit}.
597      * The String may not contain any extraneous characters (whitespace, for
598      * example).
599      *
600      * @param val decimal String representation of BigInteger.
601      * @throws NumberFormatException {@code val} is not a valid representation
602      *         of a BigInteger.
603      * @see    Character#digit
604      */

605     public BigInteger(String val) {
606         this(val, 10);
607     }
608
609     /**
610      * Constructs a randomly generated BigInteger, uniformly distributed over
611      * the range 0 to (2<sup>{@code numBits}</sup> - 1), inclusive.
612      * The uniformity of the distribution assumes that a fair source of random
613      * bits is provided in {@code rnd}.  Note that this constructor always
614      * constructs a non-negative BigInteger.
615      *
616      * @param  numBits maximum bitLength of the new BigInteger.
617      * @param  rnd source of randomness to be used in computing the new
618      *         BigInteger.
619      * @throws IllegalArgumentException {@code numBits} is negative.
620      * @see #bitLength()
621      */

622     public BigInteger(int numBits, Random rnd) {
623         this(1, randomBits(numBits, rnd));
624     }
625
626     private static byte[] randomBits(int numBits, Random rnd) {
627         if (numBits < 0)
628             throw new IllegalArgumentException("numBits must be non-negative");
629         int numBytes = (int)(((long)numBits+7)/8); // avoid overflow
630         byte[] randomBits = new byte[numBytes];
631
632         // Generate random bytes and mask out any excess bits
633         if (numBytes > 0) {
634             rnd.nextBytes(randomBits);
635             int excessBits = 8*numBytes - numBits;
636             randomBits[0] &= (1 << (8-excessBits)) - 1;
637         }
638         return randomBits;
639     }
640
641     /**
642      * Constructs a randomly generated positive BigInteger that is probably
643      * prime, with the specified bitLength.
644      *
645      * <p>It is recommended that the {@link #probablePrime probablePrime}
646      * method be used in preference to this constructor unless there
647      * is a compelling need to specify a certainty.
648      *
649      * @param  bitLength bitLength of the returned BigInteger.
650      * @param  certainty a measure of the uncertainty that the caller is
651      *         willing to tolerate.  The probability that the new BigInteger
652      *         represents a prime number will exceed
653      *         (1 - 1/2<sup>{@code certainty}</sup>).  The execution time of
654      *         this constructor is proportional to the value of this parameter.
655      * @param  rnd source of random bits used to select candidates to be
656      *         tested for primality.
657      * @throws ArithmeticException {@code bitLength < 2} or {@code bitLength} is too large.
658      * @see    #bitLength()
659      */

660     public BigInteger(int bitLength, int certainty, Random rnd) {
661         BigInteger prime;
662
663         if (bitLength < 2)
664             throw new ArithmeticException("bitLength < 2");
665         prime = (bitLength < SMALL_PRIME_THRESHOLD
666                                 ? smallPrime(bitLength, certainty, rnd)
667                                 : largePrime(bitLength, certainty, rnd));
668         signum = 1;
669         mag = prime.mag;
670     }
671
672     // Minimum size in bits that the requested prime number has
673     // before we use the large prime number generating algorithms.
674     // The cutoff of 95 was chosen empirically for best performance.
675     private static final int SMALL_PRIME_THRESHOLD = 95;
676
677     // Certainty required to meet the spec of probablePrime
678     private static final int DEFAULT_PRIME_CERTAINTY = 100;
679
680     /**
681      * Returns a positive BigInteger that is probably prime, with the
682      * specified bitLength. The probability that a BigInteger returned
683      * by this method is composite does not exceed 2<sup>-100</sup>.
684      *
685      * @param  bitLength bitLength of the returned BigInteger.
686      * @param  rnd source of random bits used to select candidates to be
687      *         tested for primality.
688      * @return a BigInteger of {@code bitLength} bits that is probably prime
689      * @throws ArithmeticException {@code bitLength < 2} or {@code bitLength} is too large.
690      * @see    #bitLength()
691      * @since 1.4
692      */

693     public static BigInteger probablePrime(int bitLength, Random rnd) {
694         if (bitLength < 2)
695             throw new ArithmeticException("bitLength < 2");
696
697         return (bitLength < SMALL_PRIME_THRESHOLD ?
698                 smallPrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd) :
699                 largePrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd));
700     }
701
702     /**
703      * Find a random number of the specified bitLength that is probably prime.
704      * This method is used for smaller primes, its performance degrades on
705      * larger bitlengths.
706      *
707      * This method assumes bitLength > 1.
708      */

709     private static BigInteger smallPrime(int bitLength, int certainty, Random rnd) {
710         int magLen = (bitLength + 31) >>> 5;
711         int temp[] = new int[magLen];
712         int highBit = 1 << ((bitLength+31) & 0x1f);  // High bit of high int
713         int highMask = (highBit << 1) - 1;  // Bits to keep in high int
714
715         while (true) {
716             // Construct a candidate
717             for (int i=0; i < magLen; i++)
718                 temp[i] = rnd.nextInt();
719             temp[0] = (temp[0] & highMask) | highBit;  // Ensure exact length
720             if (bitLength > 2)
721                 temp[magLen-1] |= 1;  // Make odd if bitlen > 2
722
723             BigInteger p = new BigInteger(temp, 1);
724
725             // Do cheap "pre-test" if applicable
726             if (bitLength > 6) {
727                 long r = p.remainder(SMALL_PRIME_PRODUCT).longValue();
728                 if ((r%3==0)  || (r%5==0)  || (r%7==0)  || (r%11==0) ||
729                     (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) ||
730                     (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0))
731                     continue// Candidate is composite; try another
732             }
733
734             // All candidates of bitLength 2 and 3 are prime by this point
735             if (bitLength < 4)
736                 return p;
737
738             // Do expensive test if we survive pre-test (or it's inapplicable)
739             if (p.primeToCertainty(certainty, rnd))
740                 return p;
741         }
742     }
743
744     private static final BigInteger SMALL_PRIME_PRODUCT
745                        = valueOf(3L*5*7*11*13*17*19*23*29*31*37*41);
746
747     /**
748      * Find a random number of the specified bitLength that is probably prime.
749      * This method is more appropriate for larger bitlengths since it uses
750      * a sieve to eliminate most composites before using a more expensive
751      * test.
752      */

753     private static BigInteger largePrime(int bitLength, int certainty, Random rnd) {
754         BigInteger p;
755         p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
756         p.mag[p.mag.length-1] &= 0xfffffffe;
757
758         // Use a sieve length likely to contain the next prime number
759         int searchLen = getPrimeSearchLen(bitLength);
760         BitSieve searchSieve = new BitSieve(p, searchLen);
761         BigInteger candidate = searchSieve.retrieve(p, certainty, rnd);
762
763         while ((candidate == null) || (candidate.bitLength() != bitLength)) {
764             p = p.add(BigInteger.valueOf(2*searchLen));
765             if (p.bitLength() != bitLength)
766                 p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
767             p.mag[p.mag.length-1] &= 0xfffffffe;
768             searchSieve = new BitSieve(p, searchLen);
769             candidate = searchSieve.retrieve(p, certainty, rnd);
770         }
771         return candidate;
772     }
773
774    /**
775     * Returns the first integer greater than this {@code BigInteger} that
776     * is probably prime.  The probability that the number returned by this
777     * method is composite does not exceed 2<sup>-100</sup>. This method will
778     * never skip over a prime when searching: if it returns {@code p}, there
779     * is no prime {@code q} such that {@code this < q < p}.
780     *
781     * @return the first integer greater than this {@code BigInteger} that
782     *         is probably prime.
783     * @throws ArithmeticException {@code this < 0} or {@code this} is too large.
784     * @since 1.5
785     */

786     public BigInteger nextProbablePrime() {
787         if (this.signum < 0)
788             throw new ArithmeticException("start < 0: " + this);
789
790         // Handle trivial cases
791         if ((this.signum == 0) || this.equals(ONE))
792             return TWO;
793
794         BigInteger result = this.add(ONE);
795
796         // Fastpath for small numbers
797         if (result.bitLength() < SMALL_PRIME_THRESHOLD) {
798
799             // Ensure an odd number
800             if (!result.testBit(0))
801                 result = result.add(ONE);
802
803             while (true) {
804                 // Do cheap "pre-test" if applicable
805                 if (result.bitLength() > 6) {
806                     long r = result.remainder(SMALL_PRIME_PRODUCT).longValue();
807                     if ((r%3==0)  || (r%5==0)  || (r%7==0)  || (r%11==0) ||
808                         (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) ||
809                         (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0)) {
810                         result = result.add(TWO);
811                         continue// Candidate is composite; try another
812                     }
813                 }
814
815                 // All candidates of bitLength 2 and 3 are prime by this point
816                 if (result.bitLength() < 4)
817                     return result;
818
819                 // The expensive test
820                 if (result.primeToCertainty(DEFAULT_PRIME_CERTAINTY, null))
821                     return result;
822
823                 result = result.add(TWO);
824             }
825         }
826
827         // Start at previous even number
828         if (result.testBit(0))
829             result = result.subtract(ONE);
830
831         // Looking for the next large prime
832         int searchLen = getPrimeSearchLen(result.bitLength());
833
834         while (true) {
835            BitSieve searchSieve = new BitSieve(result, searchLen);
836            BigInteger candidate = searchSieve.retrieve(result,
837                                                  DEFAULT_PRIME_CERTAINTY, null);
838            if (candidate != null)
839                return candidate;
840            result = result.add(BigInteger.valueOf(2 * searchLen));
841         }
842     }
843
844     private static int getPrimeSearchLen(int bitLength) {
845         if (bitLength > PRIME_SEARCH_BIT_LENGTH_LIMIT + 1) {
846             throw new ArithmeticException("Prime search implementation restriction on bitLength");
847         }
848         return bitLength / 20 * 64;
849     }
850
851     /**
852      * Returns {@code trueif this BigInteger is probably prime,
853      * {@code falseif it's definitely composite.
854      *
855      * This method assumes bitLength > 2.
856      *
857      * @param  certainty a measure of the uncertainty that the caller is
858      *         willing to tolerate: if the call returns {@code true}
859      *         the probability that this BigInteger is prime exceeds
860      *         {@code (1 - 1/2<sup>certainty</sup>)}.  The execution time of
861      *         this method is proportional to the value of this parameter.
862      * @return {@code trueif this BigInteger is probably prime,
863      *         {@code falseif it's definitely composite.
864      */

865     boolean primeToCertainty(int certainty, Random random) {
866         int rounds = 0;
867         int n = (Math.min(certainty, Integer.MAX_VALUE-1)+1)/2;
868
869         // The relationship between the certainty and the number of rounds
870         // we perform is given in the draft standard ANSI X9.80, "PRIME
871         // NUMBER GENERATION, PRIMALITY TESTING, AND PRIMALITY CERTIFICATES".
872         int sizeInBits = this.bitLength();
873         if (sizeInBits < 100) {
874             rounds = 50;
875             rounds = n < rounds ? n : rounds;
876             return passesMillerRabin(rounds, random);
877         }
878
879         if (sizeInBits < 256) {
880             rounds = 27;
881         } else if (sizeInBits < 512) {
882             rounds = 15;
883         } else if (sizeInBits < 768) {
884             rounds = 8;
885         } else if (sizeInBits < 1024) {
886             rounds = 4;
887         } else {
888             rounds = 2;
889         }
890         rounds = n < rounds ? n : rounds;
891
892         return passesMillerRabin(rounds, random) && passesLucasLehmer();
893     }
894
895     /**
896      * Returns true iff this BigInteger is a Lucas-Lehmer probable prime.
897      *
898      * The following assumptions are made:
899      * This BigInteger is a positive, odd number.
900      */

901     private boolean passesLucasLehmer() {
902         BigInteger thisPlusOne = this.add(ONE);
903
904         // Step 1
905         int d = 5;
906         while (jacobiSymbol(d, this) != -1) {
907             // 5, -7, 9, -11, ...
908             d = (d < 0) ? Math.abs(d)+2 : -(d+2);
909         }
910
911         // Step 2
912         BigInteger u = lucasLehmerSequence(d, thisPlusOne, this);
913
914         // Step 3
915         return u.mod(this).equals(ZERO);
916     }
917
918     /**
919      * Computes Jacobi(p,n).
920      * Assumes n positive, odd, n>=3.
921      */

922     private static int jacobiSymbol(int p, BigInteger n) {
923         if (p == 0)
924             return 0;
925
926         // Algorithm and comments adapted from Colin Plumb's C library.
927         int j = 1;
928         int u = n.mag[n.mag.length-1];
929
930         // Make p positive
931         if (p < 0) {
932             p = -p;
933             int n8 = u & 7;
934             if ((n8 == 3) || (n8 == 7))
935                 j = -j; // 3 (011) or 7 (111) mod 8
936         }
937
938         // Get rid of factors of 2 in p
939         while ((p & 3) == 0)
940             p >>= 2;
941         if ((p & 1) == 0) {
942             p >>= 1;
943             if (((u ^ (u>>1)) & 2) != 0)
944                 j = -j; // 3 (011) or 5 (101) mod 8
945         }
946         if (p == 1)
947             return j;
948         // Then, apply quadratic reciprocity
949         if ((p & u & 2) != 0)   // p = u = 3 (mod 4)?
950             j = -j;
951         // And reduce u mod p
952         u = n.mod(BigInteger.valueOf(p)).intValue();
953
954         // Now compute Jacobi(u,p), u < p
955         while (u != 0) {
956             while ((u & 3) == 0)
957                 u >>= 2;
958             if ((u & 1) == 0) {
959                 u >>= 1;
960                 if (((p ^ (p>>1)) & 2) != 0)
961                     j = -j;     // 3 (011) or 5 (101) mod 8
962             }
963             if (u == 1)
964                 return j;
965             // Now both u and p are odd, so use quadratic reciprocity
966             assert (u < p);
967             int t = u; u = p; p = t;
968             if ((u & p & 2) != 0) // u = p = 3 (mod 4)?
969                 j = -j;
970             // Now u >= p, so it can be reduced
971             u %= p;
972         }
973         return 0;
974     }
975
976     private static BigInteger lucasLehmerSequence(int z, BigInteger k, BigInteger n) {
977         BigInteger d = BigInteger.valueOf(z);
978         BigInteger u = ONE; BigInteger u2;
979         BigInteger v = ONE; BigInteger v2;
980
981         for (int i=k.bitLength()-2; i >= 0; i--) {
982             u2 = u.multiply(v).mod(n);
983
984             v2 = v.square().add(d.multiply(u.square())).mod(n);
985             if (v2.testBit(0))
986                 v2 = v2.subtract(n);
987
988             v2 = v2.shiftRight(1);
989
990             u = u2; v = v2;
991             if (k.testBit(i)) {
992                 u2 = u.add(v).mod(n);
993                 if (u2.testBit(0))
994                     u2 = u2.subtract(n);
995
996                 u2 = u2.shiftRight(1);
997                 v2 = v.add(d.multiply(u)).mod(n);
998                 if (v2.testBit(0))
999                     v2 = v2.subtract(n);
1000                 v2 = v2.shiftRight(1);
1001
1002                 u = u2; v = v2;
1003             }
1004         }
1005         return u;
1006     }
1007
1008     /**
1009      * Returns true iff this BigInteger passes the specified number of
1010      * Miller-Rabin tests. This test is taken from the DSA spec (NIST FIPS
1011      * 186-2).
1012      *
1013      * The following assumptions are made:
1014      * This BigInteger is a positive, odd number greater than 2.
1015      * iterations<=50.
1016      */

1017     private boolean passesMillerRabin(int iterations, Random rnd) {
1018         // Find a and m such that m is odd and this == 1 + 2**a * m
1019         BigInteger thisMinusOne = this.subtract(ONE);
1020         BigInteger m = thisMinusOne;
1021         int a = m.getLowestSetBit();
1022         m = m.shiftRight(a);
1023
1024         // Do the tests
1025         if (rnd == null) {
1026             rnd = ThreadLocalRandom.current();
1027         }
1028         for (int i=0; i < iterations; i++) {
1029             // Generate a uniform random on (1, this)
1030             BigInteger b;
1031             do {
1032                 b = new BigInteger(this.bitLength(), rnd);
1033             } while (b.compareTo(ONE) <= 0 || b.compareTo(this) >= 0);
1034
1035             int j = 0;
1036             BigInteger z = b.modPow(m, this);
1037             while (!((j == 0 && z.equals(ONE)) || z.equals(thisMinusOne))) {
1038                 if (j > 0 && z.equals(ONE) || ++j == a)
1039                     return false;
1040                 z = z.modPow(TWO, this);
1041             }
1042         }
1043         return true;
1044     }
1045
1046     /**
1047      * This internal constructor differs from its public cousin
1048      * with the arguments reversed in two ways: it assumes that its
1049      * arguments are correct, and it doesn't copy the magnitude array.
1050      */

1051     BigInteger(int[] magnitude, int signum) {
1052         this.signum = (magnitude.length == 0 ? 0 : signum);
1053         this.mag = magnitude;
1054         if (mag.length >= MAX_MAG_LENGTH) {
1055             checkRange();
1056         }
1057     }
1058
1059     /**
1060      * This private constructor is for internal use and assumes that its
1061      * arguments are correct.
1062      */

1063     private BigInteger(byte[] magnitude, int signum) {
1064         this.signum = (magnitude.length == 0 ? 0 : signum);
1065         this.mag = stripLeadingZeroBytes(magnitude);
1066         if (mag.length >= MAX_MAG_LENGTH) {
1067             checkRange();
1068         }
1069     }
1070
1071     /**
1072      * Throws an {@code ArithmeticException} if the {@code BigInteger} would be
1073      * out of the supported range.
1074      *
1075      * @throws ArithmeticException if {@code this} exceeds the supported range.
1076      */

1077     private void checkRange() {
1078         if (mag.length > MAX_MAG_LENGTH || mag.length == MAX_MAG_LENGTH && mag[0] < 0) {
1079             reportOverflow();
1080         }
1081     }
1082
1083     private static void reportOverflow() {
1084         throw new ArithmeticException("BigInteger would overflow supported range");
1085     }
1086
1087     //Static Factory Methods
1088
1089     /**
1090      * Returns a BigInteger whose value is equal to that of the
1091      * specified {@code long}.  This "static factory method" is
1092      * provided in preference to a ({@code long}) constructor
1093      * because it allows for reuse of frequently used BigIntegers.
1094      *
1095      * @param  val value of the BigInteger to return.
1096      * @return a BigInteger with the specified value.
1097      */

1098     public static BigInteger valueOf(long val) {
1099         // If -MAX_CONSTANT < val < MAX_CONSTANT, return stashed constant
1100         if (val == 0)
1101             return ZERO;
1102         if (val > 0 && val <= MAX_CONSTANT)
1103             return posConst[(int) val];
1104         else if (val < 0 && val >= -MAX_CONSTANT)
1105             return negConst[(int) -val];
1106
1107         return new BigInteger(val);
1108     }
1109
1110     /**
1111      * Constructs a BigInteger with the specified value, which may not be zero.
1112      */

1113     private BigInteger(long val) {
1114         if (val < 0) {
1115             val = -val;
1116             signum = -1;
1117         } else {
1118             signum = 1;
1119         }
1120
1121         int highWord = (int)(val >>> 32);
1122         if (highWord == 0) {
1123             mag = new int[1];
1124             mag[0] = (int)val;
1125         } else {
1126             mag = new int[2];
1127             mag[0] = highWord;
1128             mag[1] = (int)val;
1129         }
1130     }
1131
1132     /**
1133      * Returns a BigInteger with the given two's complement representation.
1134      * Assumes that the input array will not be modified (the returned
1135      * BigInteger will reference the input array if feasible).
1136      */

1137     private static BigInteger valueOf(int val[]) {
1138         return (val[0] > 0 ? new BigInteger(val, 1) : new BigInteger(val));
1139     }
1140
1141     // Constants
1142
1143     /**
1144      * Initialize static constant array when class is loaded.
1145      */

1146     private final static int MAX_CONSTANT = 16;
1147     private static BigInteger posConst[] = new BigInteger[MAX_CONSTANT+1];
1148     private static BigInteger negConst[] = new BigInteger[MAX_CONSTANT+1];
1149
1150     /**
1151      * The cache of powers of each radix.  This allows us to not have to
1152      * recalculate powers of radix^(2^n) more than once.  This speeds
1153      * Schoenhage recursive base conversion significantly.
1154      */

1155     private static volatile BigInteger[][] powerCache;
1156
1157     /** The cache of logarithms of radices for base conversion. */
1158     private static final double[] logCache;
1159
1160     /** The natural log of 2.  This is used in computing cache indices. */
1161     private static final double LOG_TWO = Math.log(2.0);
1162
1163     static {
1164         for (int i = 1; i <= MAX_CONSTANT; i++) {
1165             int[] magnitude = new int[1];
1166             magnitude[0] = i;
1167             posConst[i] = new BigInteger(magnitude,  1);
1168             negConst[i] = new BigInteger(magnitude, -1);
1169         }
1170
1171         /*
1172          * Initialize the cache of radix^(2^x) values used for base conversion
1173          * with just the very first value.  Additional values will be created
1174          * on demand.
1175          */

1176         powerCache = new BigInteger[Character.MAX_RADIX+1][];
1177         logCache = new double[Character.MAX_RADIX+1];
1178
1179         for (int i=Character.MIN_RADIX; i <= Character.MAX_RADIX; i++) {
1180             powerCache[i] = new BigInteger[] { BigInteger.valueOf(i) };
1181             logCache[i] = Math.log(i);
1182         }
1183     }
1184
1185     /**
1186      * The BigInteger constant zero.
1187      *
1188      * @since   1.2
1189      */

1190     public static final BigInteger ZERO = new BigInteger(new int[0], 0);
1191
1192     /**
1193      * The BigInteger constant one.
1194      *
1195      * @since   1.2
1196      */

1197     public static final BigInteger ONE = valueOf(1);
1198
1199     /**
1200      * The BigInteger constant two.  (Not exported.)
1201      */

1202     private static final BigInteger TWO = valueOf(2);
1203
1204     /**
1205      * The BigInteger constant -1.  (Not exported.)
1206      */

1207     private static final BigInteger NEGATIVE_ONE = valueOf(-1);
1208
1209     /**
1210      * The BigInteger constant ten.
1211      *
1212      * @since   1.5
1213      */

1214     public static final BigInteger TEN = valueOf(10);
1215
1216     // Arithmetic Operations
1217
1218     /**
1219      * Returns a BigInteger whose value is {@code (this + val)}.
1220      *
1221      * @param  val value to be added to this BigInteger.
1222      * @return {@code this + val}
1223      */

1224     public BigInteger add(BigInteger val) {
1225         if (val.signum == 0)
1226             return this;
1227         if (signum == 0)
1228             return val;
1229         if (val.signum == signum)
1230             return new BigInteger(add(mag, val.mag), signum);
1231
1232         int cmp = compareMagnitude(val);
1233         if (cmp == 0)
1234             return ZERO;
1235         int[] resultMag = (cmp > 0 ? subtract(mag, val.mag)
1236                            : subtract(val.mag, mag));
1237         resultMag = trustedStripLeadingZeroInts(resultMag);
1238
1239         return new BigInteger(resultMag, cmp == signum ? 1 : -1);
1240     }
1241
1242     /**
1243      * Package private methods used by BigDecimal code to add a BigInteger
1244      * with a long. Assumes val is not equal to INFLATED.
1245      */

1246     BigInteger add(long val) {
1247         if (val == 0)
1248             return this;
1249         if (signum == 0)
1250             return valueOf(val);
1251         if (Long.signum(val) == signum)
1252             return new BigInteger(add(mag, Math.abs(val)), signum);
1253         int cmp = compareMagnitude(val);
1254         if (cmp == 0)
1255             return ZERO;
1256         int[] resultMag = (cmp > 0 ? subtract(mag, Math.abs(val)) : subtract(Math.abs(val), mag));
1257         resultMag = trustedStripLeadingZeroInts(resultMag);
1258         return new BigInteger(resultMag, cmp == signum ? 1 : -1);
1259     }
1260
1261     /**
1262      * Adds the contents of the int array x and long value val. This
1263      * method allocates a new int array to hold the answer and returns
1264      * a reference to that array.  Assumes x.length &gt; 0 and val is
1265      * non-negative
1266      */

1267     private static int[] add(int[] x, long val) {
1268         int[] y;
1269         long sum = 0;
1270         int xIndex = x.length;
1271         int[] result;
1272         int highWord = (int)(val >>> 32);
1273         if (highWord == 0) {
1274             result = new int[xIndex];
1275             sum = (x[--xIndex] & LONG_MASK) + val;
1276             result[xIndex] = (int)sum;
1277         } else {
1278             if (xIndex == 1) {
1279                 result = new int[2];
1280                 sum = val  + (x[0] & LONG_MASK);
1281                 result[1] = (int)sum;
1282                 result[0] = (int)(sum >>> 32);
1283                 return result;
1284             } else {
1285                 result = new int[xIndex];
1286                 sum = (x[--xIndex] & LONG_MASK) + (val & LONG_MASK);
1287                 result[xIndex] = (int)sum;
1288                 sum = (x[--xIndex] & LONG_MASK) + (highWord & LONG_MASK) + (sum >>> 32);
1289                 result[xIndex] = (int)sum;
1290             }
1291         }
1292         // Copy remainder of longer number while carry propagation is required
1293         boolean carry = (sum >>> 32 != 0);
1294         while (xIndex > 0 && carry)
1295             carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
1296         // Copy remainder of longer number
1297         while (xIndex > 0)
1298             result[--xIndex] = x[xIndex];
1299         // Grow result if necessary
1300         if (carry) {
1301             int bigger[] = new int[result.length + 1];
1302             System.arraycopy(result, 0, bigger, 1, result.length);
1303             bigger[0] = 0x01;
1304             return bigger;
1305         }
1306         return result;
1307     }
1308
1309     /**
1310      * Adds the contents of the int arrays x and y. This method allocates
1311      * a new int array to hold the answer and returns a reference to that
1312      * array.
1313      */

1314     private static int[] add(int[] x, int[] y) {
1315         // If x is shorter, swap the two arrays
1316         if (x.length < y.length) {
1317             int[] tmp = x;
1318             x = y;
1319             y = tmp;
1320         }
1321
1322         int xIndex = x.length;
1323         int yIndex = y.length;
1324         int result[] = new int[xIndex];
1325         long sum = 0;
1326         if (yIndex == 1) {
1327             sum = (x[--xIndex] & LONG_MASK) + (y[0] & LONG_MASK) ;
1328             result[xIndex] = (int)sum;
1329         } else {
1330             // Add common parts of both numbers
1331             while (yIndex > 0) {
1332                 sum = (x[--xIndex] & LONG_MASK) +
1333                       (y[--yIndex] & LONG_MASK) + (sum >>> 32);
1334                 result[xIndex] = (int)sum;
1335             }
1336         }
1337         // Copy remainder of longer number while carry propagation is required
1338         boolean carry = (sum >>> 32 != 0);
1339         while (xIndex > 0 && carry)
1340             carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
1341
1342         // Copy remainder of longer number
1343         while (xIndex > 0)
1344             result[--xIndex] = x[xIndex];
1345
1346         // Grow result if necessary
1347         if (carry) {
1348             int bigger[] = new int[result.length + 1];
1349             System.arraycopy(result, 0, bigger, 1, result.length);
1350             bigger[0] = 0x01;
1351             return bigger;
1352         }
1353         return result;
1354     }
1355
1356     private static int[] subtract(long val, int[] little) {
1357         int highWord = (int)(val >>> 32);
1358         if (highWord == 0) {
1359             int result[] = new int[1];
1360             result[0] = (int)(val - (little[0] & LONG_MASK));
1361             return result;
1362         } else {
1363             int result[] = new int[2];
1364             if (little.length == 1) {
1365                 long difference = ((int)val & LONG_MASK) - (little[0] & LONG_MASK);
1366                 result[1] = (int)difference;
1367                 // Subtract remainder of longer number while borrow propagates
1368                 boolean borrow = (difference >> 32 != 0);
1369                 if (borrow) {
1370                     result[0] = highWord - 1;
1371                 } else {        // Copy remainder of longer number
1372                     result[0] = highWord;
1373                 }
1374                 return result;
1375             } else { // little.length == 2
1376                 long difference = ((int)val & LONG_MASK) - (little[1] & LONG_MASK);
1377                 result[1] = (int)difference;
1378                 difference = (highWord & LONG_MASK) - (little[0] & LONG_MASK) + (difference >> 32);
1379                 result[0] = (int)difference;
1380                 return result;
1381             }
1382         }
1383     }
1384
1385     /**
1386      * Subtracts the contents of the second argument (val) from the
1387      * first (big).  The first int array (big) must represent a larger number
1388      * than the second.  This method allocates the space necessary to hold the
1389      * answer.
1390      * assumes val &gt;= 0
1391      */

1392     private static int[] subtract(int[] big, long val) {
1393         int highWord = (int)(val >>> 32);
1394         int bigIndex = big.length;
1395         int result[] = new int[bigIndex];
1396         long difference = 0;
1397
1398         if (highWord == 0) {
1399             difference = (big[--bigIndex] & LONG_MASK) - val;
1400             result[bigIndex] = (int)difference;
1401         } else {
1402             difference = (big[--bigIndex] & LONG_MASK) - (val & LONG_MASK);
1403             result[bigIndex] = (int)difference;
1404             difference = (big[--bigIndex] & LONG_MASK) - (highWord & LONG_MASK) + (difference >> 32);
1405             result[bigIndex] = (int)difference;
1406         }
1407
1408         // Subtract remainder of longer number while borrow propagates
1409         boolean borrow = (difference >> 32 != 0);
1410         while (bigIndex > 0 && borrow)
1411             borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
1412
1413         // Copy remainder of longer number
1414         while (bigIndex > 0)
1415             result[--bigIndex] = big[bigIndex];
1416
1417         return result;
1418     }
1419
1420     /**
1421      * Returns a BigInteger whose value is {@code (this - val)}.
1422      *
1423      * @param  val value to be subtracted from this BigInteger.
1424      * @return {@code this - val}
1425      */

1426     public BigInteger subtract(BigInteger val) {
1427         if (val.signum == 0)
1428             return this;
1429         if (signum == 0)
1430             return val.negate();
1431         if (val.signum != signum)
1432             return new BigInteger(add(mag, val.mag), signum);
1433
1434         int cmp = compareMagnitude(val);
1435         if (cmp == 0)
1436             return ZERO;
1437         int[] resultMag = (cmp > 0 ? subtract(mag, val.mag)
1438                            : subtract(val.mag, mag));
1439         resultMag = trustedStripLeadingZeroInts(resultMag);
1440         return new BigInteger(resultMag, cmp == signum ? 1 : -1);
1441     }
1442
1443     /**
1444      * Subtracts the contents of the second int arrays (little) from the
1445      * first (big).  The first int array (big) must represent a larger number
1446      * than the second.  This method allocates the space necessary to hold the
1447      * answer.
1448      */

1449     private static int[] subtract(int[] big, int[] little) {
1450         int bigIndex = big.length;
1451         int result[] = new int[bigIndex];
1452         int littleIndex = little.length;
1453         long difference = 0;
1454
1455         // Subtract common parts of both numbers
1456         while (littleIndex > 0) {
1457             difference = (big[--bigIndex] & LONG_MASK) -
1458                          (little[--littleIndex] & LONG_MASK) +
1459                          (difference >> 32);
1460             result[bigIndex] = (int)difference;
1461         }
1462
1463         // Subtract remainder of longer number while borrow propagates
1464         boolean borrow = (difference >> 32 != 0);
1465         while (bigIndex > 0 && borrow)
1466             borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
1467
1468         // Copy remainder of longer number
1469         while (bigIndex > 0)
1470             result[--bigIndex] = big[bigIndex];
1471
1472         return result;
1473     }
1474
1475     /**
1476      * Returns a BigInteger whose value is {@code (this * val)}.
1477      *
1478      * @implNote An implementation may offer better algorithmic
1479      * performance when {@code val == this}.
1480      *
1481      * @param  val value to be multiplied by this BigInteger.
1482      * @return {@code this * val}
1483      */

1484     public BigInteger multiply(BigInteger val) {
1485         if (val.signum == 0 || signum == 0)
1486             return ZERO;
1487
1488         int xlen = mag.length;
1489
1490         if (val == this && xlen > MULTIPLY_SQUARE_THRESHOLD) {
1491             return square();
1492         }
1493
1494         int ylen = val.mag.length;
1495
1496         if ((xlen < KARATSUBA_THRESHOLD) || (ylen < KARATSUBA_THRESHOLD)) {
1497             int resultSign = signum == val.signum ? 1 : -1;
1498             if (val.mag.length == 1) {
1499                 return multiplyByInt(mag,val.mag[0], resultSign);
1500             }
1501             if (mag.length == 1) {
1502                 return multiplyByInt(val.mag,mag[0], resultSign);
1503             }
1504             int[] result = multiplyToLen(mag, xlen,
1505                                          val.mag, ylen, null);
1506             result = trustedStripLeadingZeroInts(result);
1507             return new BigInteger(result, resultSign);
1508         } else {
1509             if ((xlen < TOOM_COOK_THRESHOLD) && (ylen < TOOM_COOK_THRESHOLD)) {
1510                 return multiplyKaratsuba(this, val);
1511             } else {
1512                 return multiplyToomCook3(this, val);
1513             }
1514         }
1515     }
1516
1517     private static BigInteger multiplyByInt(int[] x, int y, int sign) {
1518         if (Integer.bitCount(y) == 1) {
1519             return new BigInteger(shiftLeft(x,Integer.numberOfTrailingZeros(y)), sign);
1520         }
1521         int xlen = x.length;
1522         int[] rmag =  new int[xlen + 1];
1523         long carry = 0;
1524         long yl = y & LONG_MASK;
1525         int rstart = rmag.length - 1;
1526         for (int i = xlen - 1; i >= 0; i--) {
1527             long product = (x[i] & LONG_MASK) * yl + carry;
1528             rmag[rstart--] = (int)product;
1529             carry = product >>> 32;
1530         }
1531         if (carry == 0L) {
1532             rmag = java.util.Arrays.copyOfRange(rmag, 1, rmag.length);
1533         } else {
1534             rmag[rstart] = (int)carry;
1535         }
1536         return new BigInteger(rmag, sign);
1537     }
1538
1539     /**
1540      * Package private methods used by BigDecimal code to multiply a BigInteger
1541      * with a long. Assumes v is not equal to INFLATED.
1542      */

1543     BigInteger multiply(long v) {
1544         if (v == 0 || signum == 0)
1545           return ZERO;
1546         if (v == BigDecimal.INFLATED)
1547             return multiply(BigInteger.valueOf(v));
1548         int rsign = (v > 0 ? signum : -signum);
1549         if (v < 0)
1550             v = -v;
1551         long dh = v >>> 32;      // higher order bits
1552         long dl = v & LONG_MASK; // lower order bits
1553
1554         int xlen = mag.length;
1555         int[] value = mag;
1556         int[] rmag = (dh == 0L) ? (new int[xlen + 1]) : (new int[xlen + 2]);
1557         long carry = 0;
1558         int rstart = rmag.length - 1;
1559         for (int i = xlen - 1; i >= 0; i--) {
1560             long product = (value[i] & LONG_MASK) * dl + carry;
1561             rmag[rstart--] = (int)product;
1562             carry = product >>> 32;
1563         }
1564         rmag[rstart] = (int)carry;
1565         if (dh != 0L) {
1566             carry = 0;
1567             rstart = rmag.length - 2;
1568             for (int i = xlen - 1; i >= 0; i--) {
1569                 long product = (value[i] & LONG_MASK) * dh +
1570                     (rmag[rstart] & LONG_MASK) + carry;
1571                 rmag[rstart--] = (int)product;
1572                 carry = product >>> 32;
1573             }
1574             rmag[0] = (int)carry;
1575         }
1576         if (carry == 0L)
1577             rmag = java.util.Arrays.copyOfRange(rmag, 1, rmag.length);
1578         return new BigInteger(rmag, rsign);
1579     }
1580
1581     /**
1582      * Multiplies int arrays x and y to the specified lengths and places
1583      * the result into z. There will be no leading zeros in the resultant array.
1584      */

1585     private static int[] multiplyToLen(int[] x, int xlen, int[] y, int ylen, int[] z) {
1586         int xstart = xlen - 1;
1587         int ystart = ylen - 1;
1588
1589         if (z == null || z.length < (xlen+ ylen))
1590             z = new int[xlen+ylen];
1591
1592         long carry = 0;
1593         for (int j=ystart, k=ystart+1+xstart; j >= 0; j--, k--) {
1594             long product = (y[j] & LONG_MASK) *
1595                            (x[xstart] & LONG_MASK) + carry;
1596             z[k] = (int)product;
1597             carry = product >>> 32;
1598         }
1599         z[xstart] = (int)carry;
1600
1601         for (int i = xstart-1; i >= 0; i--) {
1602             carry = 0;
1603             for (int j=ystart, k=ystart+1+i; j >= 0; j--, k--) {
1604                 long product = (y[j] & LONG_MASK) *
1605                                (x[i] & LONG_MASK) +
1606                                (z[k] & LONG_MASK) + carry;
1607                 z[k] = (int)product;
1608                 carry = product >>> 32;
1609             }
1610             z[i] = (int)carry;
1611         }
1612         return z;
1613     }
1614
1615     /**
1616      * Multiplies two BigIntegers using the Karatsuba multiplication
1617      * algorithm.  This is a recursive divide-and-conquer algorithm which is
1618      * more efficient for large numbers than what is commonly called the
1619      * "grade-school" algorithm used in multiplyToLen.  If the numbers to be
1620      * multiplied have length n, the "grade-school" algorithm has an
1621      * asymptotic complexity of O(n^2).  In contrast, the Karatsuba algorithm
1622      * has complexity of O(n^(log2(3))), or O(n^1.585).  It achieves this
1623      * increased performance by doing 3 multiplies instead of 4 when
1624      * evaluating the product.  As it has some overhead, should be used when
1625      * both numbers are larger than a certain threshold (found
1626      * experimentally).
1627      *
1628      * See:  http://en.wikipedia.org/wiki/Karatsuba_algorithm
1629      */

1630     private static BigInteger multiplyKaratsuba(BigInteger x, BigInteger y) {
1631         int xlen = x.mag.length;
1632         int ylen = y.mag.length;
1633
1634         // The number of ints in each half of the number.
1635         int half = (Math.max(xlen, ylen)+1) / 2;
1636
1637         // xl and yl are the lower halves of x and y respectively,
1638         // xh and yh are the upper halves.
1639         BigInteger xl = x.getLower(half);
1640         BigInteger xh = x.getUpper(half);
1641         BigInteger yl = y.getLower(half);
1642         BigInteger yh = y.getUpper(half);
1643
1644         BigInteger p1 = xh.multiply(yh);  // p1 = xh*yh
1645         BigInteger p2 = xl.multiply(yl);  // p2 = xl*yl
1646
1647         // p3=(xh+xl)*(yh+yl)
1648         BigInteger p3 = xh.add(xl).multiply(yh.add(yl));
1649
1650         // result = p1 * 2^(32*2*half) + (p3 - p1 - p2) * 2^(32*half) + p2
1651         BigInteger result = p1.shiftLeft(32*half).add(p3.subtract(p1).subtract(p2)).shiftLeft(32*half).add(p2);
1652
1653         if (x.signum != y.signum) {
1654             return result.negate();
1655         } else {
1656             return result;
1657         }
1658     }
1659
1660     /**
1661      * Multiplies two BigIntegers using a 3-way Toom-Cook multiplication
1662      * algorithm.  This is a recursive divide-and-conquer algorithm which is
1663      * more efficient for large numbers than what is commonly called the
1664      * "grade-school" algorithm used in multiplyToLen.  If the numbers to be
1665      * multiplied have length n, the "grade-school" algorithm has an
1666      * asymptotic complexity of O(n^2).  In contrast, 3-way Toom-Cook has a
1667      * complexity of about O(n^1.465).  It achieves this increased asymptotic
1668      * performance by breaking each number into three parts and by doing 5
1669      * multiplies instead of 9 when evaluating the product.  Due to overhead
1670      * (additions, shifts, and one division) in the Toom-Cook algorithm, it
1671      * should only be used when both numbers are larger than a certain
1672      * threshold (found experimentally).  This threshold is generally larger
1673      * than that for Karatsuba multiplication, so this algorithm is generally
1674      * only used when numbers become significantly larger.
1675      *
1676      * The algorithm used is the "optimal" 3-way Toom-Cook algorithm outlined
1677      * by Marco Bodrato.
1678      *
1679      *  See: http://bodrato.it/toom-cook/
1680      *       http://bodrato.it/papers/#WAIFI2007
1681      *
1682      * "Towards Optimal Toom-Cook Multiplication for Univariate and
1683      * Multivariate Polynomials in Characteristic 2 and 0." by Marco BODRATO;
1684      * In C.Carlet and B.Sunar, Eds., "WAIFI'07 proceedings", p. 116-133,
1685      * LNCS #4547. Springer, Madrid, Spain, June 21-22, 2007.
1686      *
1687      */

1688     private static BigInteger multiplyToomCook3(BigInteger a, BigInteger b) {
1689         int alen = a.mag.length;
1690         int blen = b.mag.length;
1691
1692         int largest = Math.max(alen, blen);
1693
1694         // k is the size (in ints) of the lower-order slices.
1695         int k = (largest+2)/3;   // Equal to ceil(largest/3)
1696
1697         // r is the size (in ints) of the highest-order slice.
1698         int r = largest - 2*k;
1699
1700         // Obtain slices of the numbers. a2 and b2 are the most significant
1701         // bits of the numbers a and b, and a0 and b0 the least significant.
1702         BigInteger a0, a1, a2, b0, b1, b2;
1703         a2 = a.getToomSlice(k, r, 0, largest);
1704         a1 = a.getToomSlice(k, r, 1, largest);
1705         a0 = a.getToomSlice(k, r, 2, largest);
1706         b2 = b.getToomSlice(k, r, 0, largest);
1707         b1 = b.getToomSlice(k, r, 1, largest);
1708         b0 = b.getToomSlice(k, r, 2, largest);
1709
1710         BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1, db1;
1711
1712         v0 = a0.multiply(b0);
1713         da1 = a2.add(a0);
1714         db1 = b2.add(b0);
1715         vm1 = da1.subtract(a1).multiply(db1.subtract(b1));
1716         da1 = da1.add(a1);
1717         db1 = db1.add(b1);
1718         v1 = da1.multiply(db1);
1719         v2 = da1.add(a2).shiftLeft(1).subtract(a0).multiply(
1720              db1.add(b2).shiftLeft(1).subtract(b0));
1721         vinf = a2.multiply(b2);
1722
1723         // The algorithm requires two divisions by 2 and one by 3.
1724         // All divisions are known to be exact, that is, they do not produce
1725         // remainders, and all results are positive.  The divisions by 2 are
1726         // implemented as right shifts which are relatively efficient, leaving
1727         // only an exact division by 3, which is done by a specialized
1728         // linear-time algorithm.
1729         t2 = v2.subtract(vm1).exactDivideBy3();
1730         tm1 = v1.subtract(vm1).shiftRight(1);
1731         t1 = v1.subtract(v0);
1732         t2 = t2.subtract(t1).shiftRight(1);
1733         t1 = t1.subtract(tm1).subtract(vinf);
1734         t2 = t2.subtract(vinf.shiftLeft(1));
1735         tm1 = tm1.subtract(t2);
1736
1737         // Number of bits to shift left.
1738         int ss = k*32;
1739
1740         BigInteger result = vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
1741
1742         if (a.signum != b.signum) {
1743             return result.negate();
1744         } else {
1745             return result;
1746         }
1747     }
1748
1749
1750     /**
1751      * Returns a slice of a BigInteger for use in Toom-Cook multiplication.
1752      *
1753      * @param lowerSize The size of the lower-order bit slices.
1754      * @param upperSize The size of the higher-order bit slices.
1755      * @param slice The index of which slice is requested, which must be a
1756      * number from 0 to size-1. Slice 0 is the highest-order bits, and slice
1757      * size-1 are the lowest-order bits. Slice 0 may be of different size than
1758      * the other slices.
1759      * @param fullsize The size of the larger integer array, used to align
1760      * slices to the appropriate position when multiplying different-sized
1761      * numbers.
1762      */

1763     private BigInteger getToomSlice(int lowerSize, int upperSize, int slice,
1764                                     int fullsize) {
1765         int start, end, sliceSize, len, offset;
1766
1767         len = mag.length;
1768         offset = fullsize - len;
1769
1770         if (slice == 0) {
1771             start = 0 - offset;
1772             end = upperSize - 1 - offset;
1773         } else {
1774             start = upperSize + (slice-1)*lowerSize - offset;
1775             end = start + lowerSize - 1;
1776         }
1777
1778         if (start < 0) {
1779             start = 0;
1780         }
1781         if (end < 0) {
1782            return ZERO;
1783         }
1784
1785         sliceSize = (end-start) + 1;
1786
1787         if (sliceSize <= 0) {
1788             return ZERO;
1789         }
1790
1791         // While performing Toom-Cook, all slices are positive and
1792         // the sign is adjusted when the final number is composed.
1793         if (start == 0 && sliceSize >= len) {
1794             return this.abs();
1795         }
1796
1797         int intSlice[] = new int[sliceSize];
1798         System.arraycopy(mag, start, intSlice, 0, sliceSize);
1799
1800         return new BigInteger(trustedStripLeadingZeroInts(intSlice), 1);
1801     }
1802
1803     /**
1804      * Does an exact division (that is, the remainder is known to be zero)
1805      * of the specified number by 3.  This is used in Toom-Cook
1806      * multiplication.  This is an efficient algorithm that runs in linear
1807      * time.  If the argument is not exactly divisible by 3, results are
1808      * undefined.  Note that this is expected to be called with positive
1809      * arguments only.
1810      */

1811     private BigInteger exactDivideBy3() {
1812         int len = mag.length;
1813         int[] result = new int[len];
1814         long x, w, q, borrow;
1815         borrow = 0L;
1816         for (int i=len-1; i >= 0; i--) {
1817             x = (mag[i] & LONG_MASK);
1818             w = x - borrow;
1819             if (borrow > x) {      // Did we make the number go negative?
1820                 borrow = 1L;
1821             } else {
1822                 borrow = 0L;
1823             }
1824
1825             // 0xAAAAAAAB is the modular inverse of 3 (mod 2^32).  Thus,
1826             // the effect of this is to divide by 3 (mod 2^32).
1827             // This is much faster than division on most architectures.
1828             q = (w * 0xAAAAAAABL) & LONG_MASK;
1829             result[i] = (int) q;
1830
1831             // Now check the borrow. The second check can of course be
1832             // eliminated if the first fails.
1833             if (q >= 0x55555556L) {
1834                 borrow++;
1835                 if (q >= 0xAAAAAAABL)
1836                     borrow++;
1837             }
1838         }
1839         result = trustedStripLeadingZeroInts(result);
1840         return new BigInteger(result, signum);
1841     }
1842
1843     /**
1844      * Returns a new BigInteger representing n lower ints of the number.
1845      * This is used by Karatsuba multiplication and Karatsuba squaring.
1846      */

1847     private BigInteger getLower(int n) {
1848         int len = mag.length;
1849
1850         if (len <= n) {
1851             return abs();
1852         }
1853
1854         int lowerInts[] = new int[n];
1855         System.arraycopy(mag, len-n, lowerInts, 0, n);
1856
1857         return new BigInteger(trustedStripLeadingZeroInts(lowerInts), 1);
1858     }
1859
1860     /**
1861      * Returns a new BigInteger representing mag.length-n upper
1862      * ints of the number.  This is used by Karatsuba multiplication and
1863      * Karatsuba squaring.
1864      */

1865     private BigInteger getUpper(int n) {
1866         int len = mag.length;
1867
1868         if (len <= n) {
1869             return ZERO;
1870         }
1871
1872         int upperLen = len - n;
1873         int upperInts[] = new int[upperLen];
1874         System.arraycopy(mag, 0, upperInts, 0, upperLen);
1875
1876         return new BigInteger(trustedStripLeadingZeroInts(upperInts), 1);
1877     }
1878
1879     // Squaring
1880
1881     /**
1882      * Returns a BigInteger whose value is {@code (this<sup>2</sup>)}.
1883      *
1884      * @return {@code this<sup>2</sup>}
1885      */

1886     private BigInteger square() {
1887         if (signum == 0) {
1888             return ZERO;
1889         }
1890         int len = mag.length;
1891
1892         if (len < KARATSUBA_SQUARE_THRESHOLD) {
1893             int[] z = squareToLen(mag, len, null);
1894             return new BigInteger(trustedStripLeadingZeroInts(z), 1);
1895         } else {
1896             if (len < TOOM_COOK_SQUARE_THRESHOLD) {
1897                 return squareKaratsuba();
1898             } else {
1899                 return squareToomCook3();
1900             }
1901         }
1902     }
1903
1904     /**
1905      * Squares the contents of the int array x. The result is placed into the
1906      * int array z.  The contents of x are not changed.
1907      */

1908     private static final int[] squareToLen(int[] x, int len, int[] z) {
1909          int zlen = len << 1;
1910          if (z == null || z.length < zlen)
1911              z = new int[zlen];
1912
1913          // Execute checks before calling intrinsified method.
1914          implSquareToLenChecks(x, len, z, zlen);
1915          return implSquareToLen(x, len, z, zlen);
1916      }
1917
1918      /**
1919       * Parameters validation.
1920       */

1921      private static void implSquareToLenChecks(int[] x, int len, int[] z, int zlen) throws RuntimeException {
1922          if (len < 1) {
1923              throw new IllegalArgumentException("invalid input length: " + len);
1924          }
1925          if (len > x.length) {
1926              throw new IllegalArgumentException("input length out of bound: " +
1927                                         len + " > " + x.length);
1928          }
1929          if (len * 2 > z.length) {
1930              throw new IllegalArgumentException("input length out of bound: " +
1931                                         (len * 2) + " > " + z.length);
1932          }
1933          if (zlen < 1) {
1934              throw new IllegalArgumentException("invalid input length: " + zlen);
1935          }
1936          if (zlen > z.length) {
1937              throw new IllegalArgumentException("input length out of bound: " +
1938                                         len + " > " + z.length);
1939          }
1940      }
1941
1942      /**
1943       * Java Runtime may use intrinsic for this method.
1944       */

1945      private static final int[] implSquareToLen(int[] x, int len, int[] z, int zlen) {
1946         /*
1947          * The algorithm used here is adapted from Colin Plumb's C library.
1948          * Technique: Consider the partial products in the multiplication
1949          * of "abcde" by itself:
1950          *
1951          *               a  b  c  d  e
1952          *            *  a  b  c  d  e
1953          *          ==================
1954          *              ae be ce de ee
1955          *           ad bd cd dd de
1956          *        ac bc cc cd ce
1957          *     ab bb bc bd be
1958          *  aa ab ac ad ae
1959          *
1960          * Note that everything above the main diagonal:
1961          *              ae be ce de = (abcd) * e
1962          *           ad bd cd       = (abc) * d
1963          *        ac bc             = (ab) * c
1964          *     ab                   = (a) * b
1965          *
1966          * is a copy of everything below the main diagonal:
1967          *                       de
1968          *                 cd ce
1969          *           bc bd be
1970          *     ab ac ad ae
1971          *
1972          * Thus, the sum is 2 * (off the diagonal) + diagonal.
1973          *
1974          * This is accumulated beginning with the diagonal (which
1975          * consist of the squares of the digits of the input), which is then
1976          * divided by two, the off-diagonal added, and multiplied by two
1977          * again.  The low bit is simply a copy of the low bit of the
1978          * input, so it doesn't need special care.
1979          */

1980
1981         // Store the squares, right shifted one bit (i.e., divided by 2)
1982         int lastProductLowWord = 0;
1983         for (int j=0, i=0; j < len; j++) {
1984             long piece = (x[j] & LONG_MASK);
1985             long product = piece * piece;
1986             z[i++] = (lastProductLowWord << 31) | (int)(product >>> 33);
1987             z[i++] = (int)(product >>> 1);
1988             lastProductLowWord = (int)product;
1989         }
1990
1991         // Add in off-diagonal sums
1992         for (int i=len, offset=1; i > 0; i--, offset+=2) {
1993             int t = x[i-1];
1994             t = mulAdd(z, x, offset, i-1, t);
1995             addOne(z, offset-1, i, t);
1996         }
1997
1998         // Shift back up and set low bit
1999         primitiveLeftShift(z, zlen, 1);
2000         z[zlen-1] |= x[len-1] & 1;
2001
2002         return z;
2003     }
2004
2005     /**
2006      * Squares a BigInteger using the Karatsuba squaring algorithm.  It should
2007      * be used when both numbers are larger than a certain threshold (found
2008      * experimentally).  It is a recursive divide-and-conquer algorithm that
2009      * has better asymptotic performance than the algorithm used in
2010      * squareToLen.
2011      */

2012     private BigInteger squareKaratsuba() {
2013         int half = (mag.length+1) / 2;
2014
2015         BigInteger xl = getLower(half);
2016         BigInteger xh = getUpper(half);
2017
2018         BigInteger xhs = xh.square();  // xhs = xh^2
2019         BigInteger xls = xl.square();  // xls = xl^2
2020
2021         // xh^2 << 64  +  (((xl+xh)^2 - (xh^2 + xl^2)) << 32) + xl^2
2022         return xhs.shiftLeft(half*32).add(xl.add(xh).square().subtract(xhs.add(xls))).shiftLeft(half*32).add(xls);
2023     }
2024
2025     /**
2026      * Squares a BigInteger using the 3-way Toom-Cook squaring algorithm.  It
2027      * should be used when both numbers are larger than a certain threshold
2028      * (found experimentally).  It is a recursive divide-and-conquer algorithm
2029      * that has better asymptotic performance than the algorithm used in
2030      * squareToLen or squareKaratsuba.
2031      */

2032     private BigInteger squareToomCook3() {
2033         int len = mag.length;
2034
2035         // k is the size (in ints) of the lower-order slices.
2036         int k = (len+2)/3;   // Equal to ceil(largest/3)
2037
2038         // r is the size (in ints) of the highest-order slice.
2039         int r = len - 2*k;
2040
2041         // Obtain slices of the numbers. a2 is the most significant
2042         // bits of the number, and a0 the least significant.
2043         BigInteger a0, a1, a2;
2044         a2 = getToomSlice(k, r, 0, len);
2045         a1 = getToomSlice(k, r, 1, len);
2046         a0 = getToomSlice(k, r, 2, len);
2047         BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1;
2048
2049         v0 = a0.square();
2050         da1 = a2.add(a0);
2051         vm1 = da1.subtract(a1).square();
2052         da1 = da1.add(a1);
2053         v1 = da1.square();
2054         vinf = a2.square();
2055         v2 = da1.add(a2).shiftLeft(1).subtract(a0).square();
2056
2057         // The algorithm requires two divisions by 2 and one by 3.
2058         // All divisions are known to be exact, that is, they do not produce
2059         // remainders, and all results are positive.  The divisions by 2 are
2060         // implemented as right shifts which are relatively efficient, leaving
2061         // only a division by 3.
2062         // The division by 3 is done by an optimized algorithm for this case.
2063         t2 = v2.subtract(vm1).exactDivideBy3();
2064         tm1 = v1.subtract(vm1).shiftRight(1);
2065         t1 = v1.subtract(v0);
2066         t2 = t2.subtract(t1).shiftRight(1);
2067         t1 = t1.subtract(tm1).subtract(vinf);
2068         t2 = t2.subtract(vinf.shiftLeft(1));
2069         tm1 = tm1.subtract(t2);
2070
2071         // Number of bits to shift left.
2072         int ss = k*32;
2073
2074         return vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
2075     }
2076
2077     // Division
2078
2079     /**
2080      * Returns a BigInteger whose value is {@code (this / val)}.
2081      *
2082      * @param  val value by which this BigInteger is to be divided.
2083      * @return {@code this / val}
2084      * @throws ArithmeticException if {@code val} is zero.
2085      */

2086     public BigInteger divide(BigInteger val) {
2087         if (val.mag.length < BURNIKEL_ZIEGLER_THRESHOLD ||
2088                 mag.length - val.mag.length < BURNIKEL_ZIEGLER_OFFSET) {
2089             return divideKnuth(val);
2090         } else {
2091             return divideBurnikelZiegler(val);
2092         }
2093     }
2094
2095     /**
2096      * Returns a BigInteger whose value is {@code (this / val)} using an O(n^2) algorithm from Knuth.
2097      *
2098      * @param  val value by which this BigInteger is to be divided.
2099      * @return {@code this / val}
2100      * @throws ArithmeticException if {@code val} is zero.
2101      * @see MutableBigInteger#divideKnuth(MutableBigInteger, MutableBigInteger, boolean)
2102      */

2103     private BigInteger divideKnuth(BigInteger val) {
2104         MutableBigInteger q = new MutableBigInteger(),
2105                           a = new MutableBigInteger(this.mag),
2106                           b = new MutableBigInteger(val.mag);
2107
2108         a.divideKnuth(b, q, false);
2109         return q.toBigInteger(this.signum * val.signum);
2110     }
2111
2112     /**
2113      * Returns an array of two BigIntegers containing {@code (this / val)}
2114      * followed by {@code (this % val)}.
2115      *
2116      * @param  val value by which this BigInteger is to be divided, and the
2117      *         remainder computed.
2118      * @return an array of two BigIntegers: the quotient {@code (this / val)}
2119      *         is the initial element, and the remainder {@code (this % val)}
2120      *         is the final element.
2121      * @throws ArithmeticException if {@code val} is zero.
2122      */

2123     public BigInteger[] divideAndRemainder(BigInteger val) {
2124         if (val.mag.length < BURNIKEL_ZIEGLER_THRESHOLD ||
2125                 mag.length - val.mag.length < BURNIKEL_ZIEGLER_OFFSET) {
2126             return divideAndRemainderKnuth(val);
2127         } else {
2128             return divideAndRemainderBurnikelZiegler(val);
2129         }
2130     }
2131
2132     /** Long division */
2133     private BigInteger[] divideAndRemainderKnuth(BigInteger val) {
2134         BigInteger[] result = new BigInteger[2];
2135         MutableBigInteger q = new MutableBigInteger(),
2136                           a = new MutableBigInteger(this.mag),
2137                           b = new MutableBigInteger(val.mag);
2138         MutableBigInteger r = a.divideKnuth(b, q);
2139         result[0] = q.toBigInteger(this.signum == val.signum ? 1 : -1);
2140         result[1] = r.toBigInteger(this.signum);
2141         return result;
2142     }
2143
2144     /**
2145      * Returns a BigInteger whose value is {@code (this % val)}.
2146      *
2147      * @param  val value by which this BigInteger is to be divided, and the
2148      *         remainder computed.
2149      * @return {@code this % val}
2150      * @throws ArithmeticException if {@code val} is zero.
2151      */

2152     public BigInteger remainder(BigInteger val) {
2153         if (val.mag.length < BURNIKEL_ZIEGLER_THRESHOLD ||
2154                 mag.length - val.mag.length < BURNIKEL_ZIEGLER_OFFSET) {
2155             return remainderKnuth(val);
2156         } else {
2157             return remainderBurnikelZiegler(val);
2158         }
2159     }
2160
2161     /** Long division */
2162     private BigInteger remainderKnuth(BigInteger val) {
2163         MutableBigInteger q = new MutableBigInteger(),
2164                           a = new MutableBigInteger(this.mag),
2165                           b = new MutableBigInteger(val.mag);
2166
2167         return a.divideKnuth(b, q).toBigInteger(this.signum);
2168     }
2169
2170     /**
2171      * Calculates {@code this / val} using the Burnikel-Ziegler algorithm.
2172      * @param  val the divisor
2173      * @return {@code this / val}
2174      */

2175     private BigInteger divideBurnikelZiegler(BigInteger val) {
2176         return divideAndRemainderBurnikelZiegler(val)[0];
2177     }
2178
2179     /**
2180      * Calculates {@code this % val} using the Burnikel-Ziegler algorithm.
2181      * @param val the divisor
2182      * @return {@code this % val}
2183      */

2184     private BigInteger remainderBurnikelZiegler(BigInteger val) {
2185         return divideAndRemainderBurnikelZiegler(val)[1];
2186     }
2187
2188     /**
2189      * Computes {@code this / val} and {@code this % val} using the
2190      * Burnikel-Ziegler algorithm.
2191      * @param val the divisor
2192      * @return an array containing the quotient and remainder
2193      */

2194     private BigInteger[] divideAndRemainderBurnikelZiegler(BigInteger val) {
2195         MutableBigInteger q = new MutableBigInteger();
2196         MutableBigInteger r = new MutableBigInteger(this).divideAndRemainderBurnikelZiegler(new MutableBigInteger(val), q);
2197         BigInteger qBigInt = q.isZero() ? ZERO : q.toBigInteger(signum*val.signum);
2198         BigInteger rBigInt = r.isZero() ? ZERO : r.toBigInteger(signum);
2199         return new BigInteger[] {qBigInt, rBigInt};
2200     }
2201
2202     /**
2203      * Returns a BigInteger whose value is <tt>(this<sup>exponent</sup>)</tt>.
2204      * Note that {@code exponent} is an integer rather than a BigInteger.
2205      *
2206      * @param  exponent exponent to which this BigInteger is to be raised.
2207      * @return <tt>this<sup>exponent</sup></tt>
2208      * @throws ArithmeticException {@code exponent} is negative.  (This would
2209      *         cause the operation to yield a non-integer value.)
2210      */

2211     public BigInteger pow(int exponent) {
2212         if (exponent < 0) {
2213             throw new ArithmeticException("Negative exponent");
2214         }
2215         if (signum == 0) {
2216             return (exponent == 0 ? ONE : this);
2217         }
2218
2219         BigInteger partToSquare = this.abs();
2220
2221         // Factor out powers of two from the base, as the exponentiation of
2222         // these can be done by left shifts only.
2223         // The remaining part can then be exponentiated faster.  The
2224         // powers of two will be multiplied back at the end.
2225         int powersOfTwo = partToSquare.getLowestSetBit();
2226         long bitsToShift = (long)powersOfTwo * exponent;
2227         if (bitsToShift > Integer.MAX_VALUE) {
2228             reportOverflow();
2229         }
2230
2231         int remainingBits;
2232
2233         // Factor the powers of two out quickly by shifting right, if needed.
2234         if (powersOfTwo > 0) {
2235             partToSquare = partToSquare.shiftRight(powersOfTwo);
2236             remainingBits = partToSquare.bitLength();
2237             if (remainingBits == 1) {  // Nothing left but +/- 1?
2238                 if (signum < 0 && (exponent&1) == 1) {
2239                     return NEGATIVE_ONE.shiftLeft(powersOfTwo*exponent);
2240                 } else {
2241                     return ONE.shiftLeft(powersOfTwo*exponent);
2242                 }
2243             }
2244         } else {
2245             remainingBits = partToSquare.bitLength();
2246             if (remainingBits == 1) { // Nothing left but +/- 1?
2247                 if (signum < 0  && (exponent&1) == 1) {
2248                     return NEGATIVE_ONE;
2249                 } else {
2250                     return ONE;
2251                 }
2252             }
2253         }
2254
2255         // This is a quick way to approximate the size of the result,
2256         // similar to doing log2[n] * exponent.  This will give an upper bound
2257         // of how big the result can be, and which algorithm to use.
2258         long scaleFactor = (long)remainingBits * exponent;
2259
2260         // Use slightly different algorithms for small and large operands.
2261         // See if the result will safely fit into a long. (Largest 2^63-1)
2262         if (partToSquare.mag.length == 1 && scaleFactor <= 62) {
2263             // Small number algorithm.  Everything fits into a long.
2264             int newSign = (signum <0  && (exponent&1) == 1 ? -1 : 1);
2265             long result = 1;
2266             long baseToPow2 = partToSquare.mag[0] & LONG_MASK;
2267
2268             int workingExponent = exponent;
2269
2270             // Perform exponentiation using repeated squaring trick
2271             while (workingExponent != 0) {
2272                 if ((workingExponent & 1) == 1) {
2273                     result = result * baseToPow2;
2274                 }
2275
2276                 if ((workingExponent >>>= 1) != 0) {
2277                     baseToPow2 = baseToPow2 * baseToPow2;
2278                 }
2279             }
2280
2281             // Multiply back the powers of two (quickly, by shifting left)
2282             if (powersOfTwo > 0) {
2283                 if (bitsToShift + scaleFactor <= 62) { // Fits in long?
2284                     return valueOf((result << bitsToShift) * newSign);
2285                 } else {
2286                     return valueOf(result*newSign).shiftLeft((int) bitsToShift);
2287                 }
2288             }
2289             else {
2290                 return valueOf(result*newSign);
2291             }
2292         } else {
2293             // Large number algorithm.  This is basically identical to
2294             // the algorithm above, but calls multiply() and square()
2295             // which may use more efficient algorithms for large numbers.
2296             BigInteger answer = ONE;
2297
2298             int workingExponent = exponent;
2299             // Perform exponentiation using repeated squaring trick
2300             while (workingExponent != 0) {
2301                 if ((workingExponent & 1) == 1) {
2302                     answer = answer.multiply(partToSquare);
2303                 }
2304
2305                 if ((workingExponent >>>= 1) != 0) {
2306                     partToSquare = partToSquare.square();
2307                 }
2308             }
2309             // Multiply back the (exponentiated) powers of two (quickly,
2310             // by shifting left)
2311             if (powersOfTwo > 0) {
2312                 answer = answer.shiftLeft(powersOfTwo*exponent);
2313             }
2314
2315             if (signum < 0 && (exponent&1) == 1) {
2316                 return answer.negate();
2317             } else {
2318                 return answer;
2319             }
2320         }
2321     }
2322
2323     /**
2324      * Returns a BigInteger whose value is the greatest common divisor of
2325      * {@code abs(this)} and {@code abs(val)}.  Returns 0 if
2326      * {@code this == 0 && val == 0}.
2327      *
2328      * @param  val value with which the GCD is to be computed.
2329      * @return {@code GCD(abs(this), abs(val))}
2330      */

2331     public BigInteger gcd(BigInteger val) {
2332         if (val.signum == 0)
2333             return this.abs();
2334         else if (this.signum == 0)
2335             return val.abs();
2336
2337         MutableBigInteger a = new MutableBigInteger(this);
2338         MutableBigInteger b = new MutableBigInteger(val);
2339
2340         MutableBigInteger result = a.hybridGCD(b);
2341
2342         return result.toBigInteger(1);
2343     }
2344
2345     /**
2346      * Package private method to return bit length for an integer.
2347      */

2348     static int bitLengthForInt(int n) {
2349         return 32 - Integer.numberOfLeadingZeros(n);
2350     }
2351
2352     /**
2353      * Left shift int array a up to len by n bits. Returns the array that
2354      * results from the shift since space may have to be reallocated.
2355      */

2356     private static int[] leftShift(int[] a, int len, int n) {
2357         int nInts = n >>> 5;
2358         int nBits = n&0x1F;
2359         int bitsInHighWord = bitLengthForInt(a[0]);
2360
2361         // If shift can be done without recopy, do so
2362         if (n <= (32-bitsInHighWord)) {
2363             primitiveLeftShift(a, len, nBits);
2364             return a;
2365         } else { // Array must be resized
2366             if (nBits <= (32-bitsInHighWord)) {
2367                 int result[] = new int[nInts+len];
2368                 System.arraycopy(a, 0, result, 0, len);
2369                 primitiveLeftShift(result, result.length, nBits);
2370                 return result;
2371             } else {
2372                 int result[] = new int[nInts+len+1];
2373                 System.arraycopy(a, 0, result, 0, len);
2374                 primitiveRightShift(result, result.length, 32 - nBits);
2375                 return result;
2376             }
2377         }
2378     }
2379
2380     // shifts a up to len right n bits assumes no leading zeros, 0<n<32
2381     static void primitiveRightShift(int[] a, int len, int n) {
2382         int n2 = 32 - n;
2383         for (int i=len-1, c=a[i]; i > 0; i--) {
2384             int b = c;
2385             c = a[i-1];
2386             a[i] = (c << n2) | (b >>> n);
2387         }
2388         a[0] >>>= n;
2389     }
2390
2391     // shifts a up to len left n bits assumes no leading zeros, 0<=n<32
2392     static void primitiveLeftShift(int[] a, int len, int n) {
2393         if (len == 0 || n == 0)
2394             return;
2395
2396         int n2 = 32 - n;
2397         for (int i=0, c=a[i], m=i+len-1; i < m; i++) {
2398             int b = c;
2399             c = a[i+1];
2400             a[i] = (b << n) | (c >>> n2);
2401         }
2402         a[len-1] <<= n;
2403     }
2404
2405     /**
2406      * Calculate bitlength of contents of the first len elements an int array,
2407      * assuming there are no leading zero ints.
2408      */

2409     private static int bitLength(int[] val, int len) {
2410         if (len == 0)
2411             return 0;
2412         return ((len - 1) << 5) + bitLengthForInt(val[0]);
2413     }
2414
2415     /**
2416      * Returns a BigInteger whose value is the absolute value of this
2417      * BigInteger.
2418      *
2419      * @return {@code abs(this)}
2420      */

2421     public BigInteger abs() {
2422         return (signum >= 0 ? this : this.negate());
2423     }
2424
2425     /**
2426      * Returns a BigInteger whose value is {@code (-this)}.
2427      *
2428      * @return {@code -this}
2429      */

2430     public BigInteger negate() {
2431         return new BigInteger(this.mag, -this.signum);
2432     }
2433
2434     /**
2435      * Returns the signum function of this BigInteger.
2436      *
2437      * @return -1, 0 or 1 as the value of this BigInteger is negative, zero or
2438      *         positive.
2439      */

2440     public int signum() {
2441         return this.signum;
2442     }
2443
2444     // Modular Arithmetic Operations
2445
2446     /**
2447      * Returns a BigInteger whose value is {@code (this mod m}).  This method
2448      * differs from {@code remainder} in that it always returns a
2449      * <i>non-negative</i> BigInteger.
2450      *
2451      * @param  m the modulus.
2452      * @return {@code this mod m}
2453      * @throws ArithmeticException {@code m} &le; 0
2454      * @see    #remainder
2455      */

2456     public BigInteger mod(BigInteger m) {
2457         if (m.signum <= 0)
2458             throw new ArithmeticException("BigInteger: modulus not positive");
2459
2460         BigInteger result = this.remainder(m);
2461         return (result.signum >= 0 ? result : result.add(m));
2462     }
2463
2464     /**
2465      * Returns a BigInteger whose value is
2466      * <tt>(this<sup>exponent</sup> mod m)</tt>.  (Unlike {@code pow}, this
2467      * method permits negative exponents.)
2468      *
2469      * @param  exponent the exponent.
2470      * @param  m the modulus.
2471      * @return <tt>this<sup>exponent</sup> mod m</tt>
2472      * @throws ArithmeticException {@code m} &le; 0 or the exponent is
2473      *         negative and this BigInteger is not <i>relatively
2474      *         prime</i> to {@code m}.
2475      * @see    #modInverse
2476      */

2477     public BigInteger modPow(BigInteger exponent, BigInteger m) {
2478         if (m.signum <= 0)
2479             throw new ArithmeticException("BigInteger: modulus not positive");
2480
2481         // Trivial cases
2482         if (exponent.signum == 0)
2483             return (m.equals(ONE) ? ZERO : ONE);
2484
2485         if (this.equals(ONE))
2486             return (m.equals(ONE) ? ZERO : ONE);
2487
2488         if (this.equals(ZERO) && exponent.signum >= 0)
2489             return ZERO;
2490
2491         if (this.equals(negConst[1]) && (!exponent.testBit(0)))
2492             return (m.equals(ONE) ? ZERO : ONE);
2493
2494         boolean invertResult;
2495         if ((invertResult = (exponent.signum < 0)))
2496             exponent = exponent.negate();
2497
2498         BigInteger base = (this.signum < 0 || this.compareTo(m) >= 0
2499                            ? this.mod(m) : this);
2500         BigInteger result;
2501         if (m.testBit(0)) { // odd modulus
2502             result = base.oddModPow(exponent, m);
2503         } else {
2504             /*
2505              * Even modulus.  Tear it into an "odd part" (m1) and power of two
2506              * (m2), exponentiate mod m1, manually exponentiate mod m2, and
2507              * use Chinese Remainder Theorem to combine results.
2508              */

2509
2510             // Tear m apart into odd part (m1) and power of 2 (m2)
2511             int p = m.getLowestSetBit();   // Max pow of 2 that divides m
2512
2513             BigInteger m1 = m.shiftRight(p);  // m/2**p
2514             BigInteger m2 = ONE.shiftLeft(p); // 2**p
2515
2516             // Calculate new base from m1
2517             BigInteger base2 = (this.signum < 0 || this.compareTo(m1) >= 0
2518                                 ? this.mod(m1) : this);
2519
2520             // Caculate (base ** exponent) mod m1.
2521             BigInteger a1 = (m1.equals(ONE) ? ZERO :
2522                              base2.oddModPow(exponent, m1));
2523
2524             // Calculate (this ** exponent) mod m2
2525             BigInteger a2 = base.modPow2(exponent, p);
2526
2527             // Combine results using Chinese Remainder Theorem
2528             BigInteger y1 = m2.modInverse(m1);
2529             BigInteger y2 = m1.modInverse(m2);
2530
2531             if (m.mag.length < MAX_MAG_LENGTH / 2) {
2532                 result = a1.multiply(m2).multiply(y1).add(a2.multiply(m1).multiply(y2)).mod(m);
2533             } else {
2534                 MutableBigInteger t1 = new MutableBigInteger();
2535                 new MutableBigInteger(a1.multiply(m2)).multiply(new MutableBigInteger(y1), t1);
2536                 MutableBigInteger t2 = new MutableBigInteger();
2537                 new MutableBigInteger(a2.multiply(m1)).multiply(new MutableBigInteger(y2), t2);
2538                 t1.add(t2);
2539                 MutableBigInteger q = new MutableBigInteger();
2540                 result = t1.divide(new MutableBigInteger(m), q).toBigInteger();
2541             }
2542         }
2543
2544         return (invertResult ? result.modInverse(m) : result);
2545     }
2546
2547     // Montgomery multiplication.  These are wrappers for
2548     // implMontgomeryXX routines which are expected to be replaced by
2549     // virtual machine intrinsics.  We don't use the intrinsics for
2550     // very large operands: MONTGOMERY_INTRINSIC_THRESHOLD should be
2551     // larger than any reasonable crypto key.
2552     private static int[] montgomeryMultiply(int[] a, int[] b, int[] n, int len, long inv,
2553                                             int[] product) {
2554         implMontgomeryMultiplyChecks(a, b, n, len, product);
2555         if (len > MONTGOMERY_INTRINSIC_THRESHOLD) {
2556             // Very long argument: do not use an intrinsic
2557             product = multiplyToLen(a, len, b, len, product);
2558             return montReduce(product, n, len, (int)inv);
2559         } else {
2560             return implMontgomeryMultiply(a, b, n, len, inv, materialize(product, len));
2561         }
2562     }
2563     private static int[] montgomerySquare(int[] a, int[] n, int len, long inv,
2564                                           int[] product) {
2565         implMontgomeryMultiplyChecks(a, a, n, len, product);
2566         if (len > MONTGOMERY_INTRINSIC_THRESHOLD) {
2567             // Very long argument: do not use an intrinsic
2568             product = squareToLen(a, len, product);
2569             return montReduce(product, n, len, (int)inv);
2570         } else {
2571             return implMontgomerySquare(a, n, len, inv, materialize(product, len));
2572         }
2573     }
2574
2575     // Range-check everything.
2576     private static void implMontgomeryMultiplyChecks
2577         (int[] a, int[] b, int[] n, int len, int[] product) throws RuntimeException {
2578         if (len % 2 != 0) {
2579             throw new IllegalArgumentException("input array length must be even: " + len);
2580         }
2581
2582         if (len < 1) {
2583             throw new IllegalArgumentException("invalid input length: " + len);
2584         }
2585
2586         if (len > a.length ||
2587             len > b.length ||
2588             len > n.length ||
2589             (product != null && len > product.length)) {
2590             throw new IllegalArgumentException("input array length out of bound: " + len);
2591         }
2592     }
2593
2594     // Make sure that the int array z (which is expected to contain
2595     // the result of a Montgomery multiplication) is present and
2596     // sufficiently large.
2597     private static int[] materialize(int[] z, int len) {
2598          if (z == null || z.length < len)
2599              z = new int[len];
2600          return z;
2601     }
2602
2603     // These methods are intended to be be replaced by virtual machine
2604     // intrinsics.
2605     private static int[] implMontgomeryMultiply(int[] a, int[] b, int[] n, int len,
2606                                          long inv, int[] product) {
2607         product = multiplyToLen(a, len, b, len, product);
2608         return montReduce(product, n, len, (int)inv);
2609     }
2610     private static int[] implMontgomerySquare(int[] a, int[] n, int len,
2611                                        long inv, int[] product) {
2612         product = squareToLen(a, len, product);
2613         return montReduce(product, n, len, (int)inv);
2614     }
2615
2616     static int[] bnExpModThreshTable = {7, 25, 81, 241, 673, 1793,
2617                                                 Integer.MAX_VALUE}; // Sentinel
2618
2619     /**
2620      * Returns a BigInteger whose value is x to the power of y mod z.
2621      * Assumes: z is odd && x < z.
2622      */

2623     private BigInteger oddModPow(BigInteger y, BigInteger z) {
2624     /*
2625      * The algorithm is adapted from Colin Plumb's C library.
2626      *
2627      * The window algorithm:
2628      * The idea is to keep a running product of b1 = n^(high-order bits of exp)
2629      * and then keep appending exponent bits to it.  The following patterns
2630      * apply to a 3-bit window (k = 3):
2631      * To append   0: square
2632      * To append   1: square, multiply by n^1
2633      * To append  10: square, multiply by n^1, square
2634      * To append  11: square, square, multiply by n^3
2635      * To append 100: square, multiply by n^1, square, square
2636      * To append 101: square, square, square, multiply by n^5
2637      * To append 110: square, square, multiply by n^3, square
2638      * To append 111: square, square, square, multiply by n^7
2639      *
2640      * Since each pattern involves only one multiply, the longer the pattern
2641      * the better, except that a 0 (no multiplies) can be appended directly.
2642      * We precompute a table of odd powers of n, up to 2^k, and can then
2643      * multiply k bits of exponent at a time.  Actually, assuming random
2644      * exponents, there is on average one zero bit between needs to
2645      * multiply (1/2 of the time there's none, 1/4 of the time there's 1,
2646      * 1/8 of the time, there's 2, 1/32 of the time, there's 3, etc.), so
2647      * you have to do one multiply per k+1 bits of exponent.
2648      *
2649      * The loop walks down the exponent, squaring the result buffer as
2650      * it goes.  There is a wbits+1 bit lookahead buffer, buf, that is
2651      * filled with the upcoming exponent bits.  (What is read after the
2652      * end of the exponent is unimportant, but it is filled with zero here.)
2653      * When the most-significant bit of this buffer becomes set, i.e.
2654      * (buf & tblmask) != 0, we have to decide what pattern to multiply
2655      * by, and when to do it.  We decide, remember to do it in future
2656      * after a suitable number of squarings have passed (e.g. a pattern
2657      * of "100" in the buffer requires that we multiply by n^1 immediately;
2658      * a pattern of "110" calls for multiplying by n^3 after one more
2659      * squaring), clear the buffer, and continue.
2660      *
2661      * When we start, there is one more optimization: the result buffer
2662      * is implcitly one, so squaring it or multiplying by it can be
2663      * optimized away.  Further, if we start with a pattern like "100"
2664      * in the lookahead window, rather than placing n into the buffer
2665      * and then starting to square it, we have already computed n^2
2666      * to compute the odd-powers table, so we can place that into
2667      * the buffer and save a squaring.
2668      *
2669      * This means that if you have a k-bit window, to compute n^z,
2670      * where z is the high k bits of the exponent, 1/2 of the time
2671      * it requires no squarings.  1/4 of the time, it requires 1
2672      * squaring, ... 1/2^(k-1) of the time, it reqires k-2 squarings.
2673      * And the remaining 1/2^(k-1) of the time, the top k bits are a
2674      * 1 followed by k-1 0 bits, so it again only requires k-2
2675      * squarings, not k-1.  The average of these is 1.  Add that
2676      * to the one squaring we have to do to compute the table,
2677      * and you'll see that a k-bit window saves k-2 squarings
2678      * as well as reducing the multiplies.  (It actually doesn't
2679      * hurt in the case k = 1, either.)
2680      */

2681         // Special case for exponent of one
2682         if (y.equals(ONE))
2683             return this;
2684
2685         // Special case for base of zero
2686         if (signum == 0)
2687             return ZERO;
2688
2689         int[] base = mag.clone();
2690         int[] exp = y.mag;
2691         int[] mod = z.mag;
2692         int modLen = mod.length;
2693
2694         // Make modLen even. It is conventional to use a cryptographic
2695         // modulus that is 512, 768, 1024, or 2048 bits, so this code
2696         // will not normally be executed. However, it is necessary for
2697         // the correct functioning of the HotSpot intrinsics.
2698         if ((modLen & 1) != 0) {
2699             int[] x = new int[modLen + 1];
2700             System.arraycopy(mod, 0, x, 1, modLen);
2701             mod = x;
2702             modLen++;
2703         }
2704
2705         // Select an appropriate window size
2706         int wbits = 0;
2707         int ebits = bitLength(exp, exp.length);
2708         // if exponent is 65537 (0x10001), use minimum window size
2709         if ((ebits != 17) || (exp[0] != 65537)) {
2710             while (ebits > bnExpModThreshTable[wbits]) {
2711                 wbits++;
2712             }
2713         }
2714
2715         // Calculate appropriate table size
2716         int tblmask = 1 << wbits;
2717
2718         // Allocate table for precomputed odd powers of base in Montgomery form
2719         int[][] table = new int[tblmask][];
2720         for (int i=0; i < tblmask; i++)
2721             table[i] = new int[modLen];
2722
2723         // Compute the modular inverse of the least significant 64-bit
2724         // digit of the modulus
2725         long n0 = (mod[modLen-1] & LONG_MASK) + ((mod[modLen-2] & LONG_MASK) << 32);
2726         long inv = -MutableBigInteger.inverseMod64(n0);
2727
2728         // Convert base to Montgomery form
2729         int[] a = leftShift(base, base.length, modLen << 5);
2730
2731         MutableBigInteger q = new MutableBigInteger(),
2732                           a2 = new MutableBigInteger(a),
2733                           b2 = new MutableBigInteger(mod);
2734         b2.normalize(); // MutableBigInteger.divide() assumes that its
2735                         // divisor is in normal form.
2736
2737         MutableBigInteger r= a2.divide(b2, q);
2738         table[0] = r.toIntArray();
2739
2740         // Pad table[0] with leading zeros so its length is at least modLen
2741         if (table[0].length < modLen) {
2742            int offset = modLen - table[0].length;
2743            int[] t2 = new int[modLen];
2744            System.arraycopy(table[0], 0, t2, offset, table[0].length);
2745            table[0] = t2;
2746         }
2747
2748         // Set b to the square of the base
2749         int[] b = montgomerySquare(table[0], mod, modLen, inv, null);
2750
2751         // Set t to high half of b
2752         int[] t = Arrays.copyOf(b, modLen);
2753
2754         // Fill in the table with odd powers of the base
2755         for (int i=1; i < tblmask; i++) {
2756             table[i] = montgomeryMultiply(t, table[i-1], mod, modLen, inv, null);
2757         }
2758
2759         // Pre load the window that slides over the exponent
2760         int bitpos = 1 << ((ebits-1) & (32-1));
2761
2762         int buf = 0;
2763         int elen = exp.length;
2764         int eIndex = 0;
2765         for (int i = 0; i <= wbits; i++) {
2766             buf = (buf << 1) | (((exp[eIndex] & bitpos) != 0)?1:0);
2767             bitpos >>>= 1;
2768             if (bitpos == 0) {
2769                 eIndex++;
2770                 bitpos = 1 << (32-1);
2771                 elen--;
2772             }
2773         }
2774
2775         int multpos = ebits;
2776
2777         // The first iteration, which is hoisted out of the main loop
2778         ebits--;
2779         boolean isone = true;
2780
2781         multpos = ebits - wbits;
2782         while ((buf & 1) == 0) {
2783             buf >>>= 1;
2784             multpos++;
2785         }
2786
2787         int[] mult = table[buf >>> 1];
2788
2789         buf = 0;
2790         if (multpos == ebits)
2791             isone = false;
2792
2793         // The main loop
2794         while (true) {
2795             ebits--;
2796             // Advance the window
2797             buf <<= 1;
2798
2799             if (elen != 0) {
2800                 buf |= ((exp[eIndex] & bitpos) != 0) ? 1 : 0;
2801                 bitpos >>>= 1;
2802                 if (bitpos == 0) {
2803                     eIndex++;
2804                     bitpos = 1 << (32-1);
2805                     elen--;
2806                 }
2807             }
2808
2809             // Examine the window for pending multiplies
2810             if ((buf & tblmask) != 0) {
2811                 multpos = ebits - wbits;
2812                 while ((buf & 1) == 0) {
2813                     buf >>>= 1;
2814                     multpos++;
2815                 }
2816                 mult = table[buf >>> 1];
2817                 buf = 0;
2818             }
2819
2820             // Perform multiply
2821             if (ebits == multpos) {
2822                 if (isone) {
2823                     b = mult.clone();
2824                     isone = false;
2825                 } else {
2826                     t = b;
2827                     a = montgomeryMultiply(t, mult, mod, modLen, inv, a);
2828                     t = a; a = b; b = t;
2829                 }
2830             }
2831
2832             // Check if done
2833             if (ebits == 0)
2834                 break;
2835
2836             // Square the input
2837             if (!isone) {
2838                 t = b;
2839                 a = montgomerySquare(t, mod, modLen, inv, a);
2840                 t = a; a = b; b = t;
2841             }
2842         }
2843
2844         // Convert result out of Montgomery form and return
2845         int[] t2 = new int[2*modLen];
2846         System.arraycopy(b, 0, t2, modLen, modLen);
2847
2848         b = montReduce(t2, mod, modLen, (int)inv);
2849
2850         t2 = Arrays.copyOf(b, modLen);
2851
2852         return new BigInteger(1, t2);
2853     }
2854
2855     /**
2856      * Montgomery reduce n, modulo mod.  This reduces modulo mod and divides
2857      * by 2^(32*mlen). Adapted from Colin Plumb's C library.
2858      */

2859     private static int[] montReduce(int[] n, int[] mod, int mlen, int inv) {
2860         int c=0;
2861         int len = mlen;
2862         int offset=0;
2863
2864         do {
2865             int nEnd = n[n.length-1-offset];
2866             int carry = mulAdd(n, mod, offset, mlen, inv * nEnd);
2867             c += addOne(n, offset, mlen, carry);
2868             offset++;
2869         } while (--len > 0);
2870
2871         while (c > 0)
2872             c += subN(n, mod, mlen);
2873
2874         while (intArrayCmpToLen(n, mod, mlen) >= 0)
2875             subN(n, mod, mlen);
2876
2877         return n;
2878     }
2879
2880
2881     /*
2882      * Returns -1, 0 or +1 as big-endian unsigned int array arg1 is less than,
2883      * equal to, or greater than arg2 up to length len.
2884      */

2885     private static int intArrayCmpToLen(int[] arg1, int[] arg2, int len) {
2886         for (int i=0; i < len; i++) {
2887             long b1 = arg1[i] & LONG_MASK;
2888             long b2 = arg2[i] & LONG_MASK;
2889             if (b1 < b2)
2890                 return -1;
2891             if (b1 > b2)
2892                 return 1;
2893         }
2894         return 0;
2895     }
2896
2897     /**
2898      * Subtracts two numbers of same length, returning borrow.
2899      */

2900     private static int subN(int[] a, int[] b, int len) {
2901         long sum = 0;
2902
2903         while (--len >= 0) {
2904             sum = (a[len] & LONG_MASK) -
2905                  (b[len] & LONG_MASK) + (sum >> 32);
2906             a[len] = (int)sum;
2907         }
2908
2909         return (int)(sum >> 32);
2910     }
2911
2912     /**
2913      * Multiply an array by one word k and add to result, return the carry
2914      */

2915     static int mulAdd(int[] out, int[] in, int offset, int len, int k) {
2916         implMulAddCheck(out, in, offset, len, k);
2917         return implMulAdd(out, in, offset, len, k);
2918     }
2919
2920     /**
2921      * Parameters validation.
2922      */

2923     private static void implMulAddCheck(int[] out, int[] in, int offset, int len, int k) {
2924         if (len > in.length) {
2925             throw new IllegalArgumentException("input length is out of bound: " + len + " > " + in.length);
2926         }
2927         if (offset < 0) {
2928             throw new IllegalArgumentException("input offset is invalid: " + offset);
2929         }
2930         if (offset > (out.length - 1)) {
2931             throw new IllegalArgumentException("input offset is out of bound: " + offset + " > " + (out.length - 1));
2932         }
2933         if (len > (out.length - offset)) {
2934             throw new IllegalArgumentException("input len is out of bound: " + len + " > " + (out.length - offset));
2935         }
2936     }
2937
2938     /**
2939      * Java Runtime may use intrinsic for this method.
2940      */

2941     private static int implMulAdd(int[] out, int[] in, int offset, int len, int k) {
2942         long kLong = k & LONG_MASK;
2943         long carry = 0;
2944
2945         offset = out.length-offset - 1;
2946         for (int j=len-1; j >= 0; j--) {
2947             long product = (in[j] & LONG_MASK) * kLong +
2948                            (out[offset] & LONG_MASK) + carry;
2949             out[offset--] = (int)product;
2950             carry = product >>> 32;
2951         }
2952         return (int)carry;
2953     }
2954
2955     /**
2956      * Add one word to the number a mlen words into a. Return the resulting
2957      * carry.
2958      */

2959     static int addOne(int[] a, int offset, int mlen, int carry) {
2960         offset = a.length-1-mlen-offset;
2961         long t = (a[offset] & LONG_MASK) + (carry & LONG_MASK);
2962
2963         a[offset] = (int)t;
2964         if ((t >>> 32) == 0)
2965             return 0;
2966         while (--mlen >= 0) {
2967             if (--offset < 0) { // Carry out of number
2968                 return 1;
2969             } else {
2970                 a[offset]++;
2971                 if (a[offset] != 0)
2972                     return 0;
2973             }
2974         }
2975         return 1;
2976     }
2977
2978     /**
2979      * Returns a BigInteger whose value is (this ** exponent) mod (2**p)
2980      */

2981     private BigInteger modPow2(BigInteger exponent, int p) {
2982         /*
2983          * Perform exponentiation using repeated squaring trick, chopping off
2984          * high order bits as indicated by modulus.
2985          */

2986         BigInteger result = ONE;
2987         BigInteger baseToPow2 = this.mod2(p);
2988         int expOffset = 0;
2989
2990         int limit = exponent.bitLength();
2991
2992         if (this.testBit(0))
2993            limit = (p-1) < limit ? (p-1) : limit;
2994
2995         while (expOffset < limit) {
2996             if (exponent.testBit(expOffset))
2997                 result = result.multiply(baseToPow2).mod2(p);
2998             expOffset++;
2999             if (expOffset < limit)
3000                 baseToPow2 = baseToPow2.square().mod2(p);
3001         }
3002
3003         return result;
3004     }
3005
3006     /**
3007      * Returns a BigInteger whose value is this mod(2**p).
3008      * Assumes that this {@code BigInteger >= 0} and {@code p > 0}.
3009      */

3010     private BigInteger mod2(int p) {
3011         if (bitLength() <= p)
3012             return this;
3013
3014         // Copy remaining ints of mag
3015         int numInts = (p + 31) >>> 5;
3016         int[] mag = new int[numInts];
3017         System.arraycopy(this.mag, (this.mag.length - numInts), mag, 0, numInts);
3018
3019         // Mask out any excess bits
3020         int excessBits = (numInts << 5) - p;
3021         mag[0] &= (1L << (32-excessBits)) - 1;
3022
3023         return (mag[0] == 0 ? new BigInteger(1, mag) : new BigInteger(mag, 1));
3024     }
3025
3026     /**
3027      * Returns a BigInteger whose value is {@code (this}<sup>-1</sup> {@code mod m)}.
3028      *
3029      * @param  m the modulus.
3030      * @return {@code this}<sup>-1</sup> {@code mod m}.
3031      * @throws ArithmeticException {@code  m} &le; 0, or this BigInteger
3032      *         has no multiplicative inverse mod m (that is, this BigInteger
3033      *         is not <i>relatively prime</i> to m).
3034      */

3035     public BigInteger modInverse(BigInteger m) {
3036         if (m.signum != 1)
3037             throw new ArithmeticException("BigInteger: modulus not positive");
3038
3039         if (m.equals(ONE))
3040             return ZERO;
3041
3042         // Calculate (this mod m)
3043         BigInteger modVal = this;
3044         if (signum < 0 || (this.compareMagnitude(m) >= 0))
3045             modVal = this.mod(m);
3046
3047         if (modVal.equals(ONE))
3048             return ONE;
3049
3050         MutableBigInteger a = new MutableBigInteger(modVal);
3051         MutableBigInteger b = new MutableBigInteger(m);
3052
3053         MutableBigInteger result = a.mutableModInverse(b);
3054         return result.toBigInteger(1);
3055     }
3056
3057     // Shift Operations
3058
3059     /**
3060      * Returns a BigInteger whose value is {@code (this << n)}.
3061      * The shift distance, {@code n}, may be negative, in which case
3062      * this method performs a right shift.
3063      * (Computes <tt>floor(this * 2<sup>n</sup>)</tt>.)
3064      *
3065      * @param  n shift distance, in bits.
3066      * @return {@code this << n}
3067      * @see #shiftRight
3068      */

3069     public BigInteger shiftLeft(int n) {
3070         if (signum == 0)
3071             return ZERO;
3072         if (n > 0) {
3073             return new BigInteger(shiftLeft(mag, n), signum);
3074         } else if (n == 0) {
3075             return this;
3076         } else {
3077             // Possible int overflow in (-n) is not a trouble,
3078             // because shiftRightImpl considers its argument unsigned
3079             return shiftRightImpl(-n);
3080         }
3081     }
3082
3083     /**
3084      * Returns a magnitude array whose value is {@code (mag << n)}.
3085      * The shift distance, {@code n}, is considered unnsigned.
3086      * (Computes <tt>this * 2<sup>n</sup></tt>.)
3087      *
3088      * @param mag magnitude, the most-significant int ({@code mag[0]}) must be non-zero.
3089      * @param  n unsigned shift distance, in bits.
3090      * @return {@code mag << n}
3091      */

3092     private static int[] shiftLeft(int[] mag, int n) {
3093         int nInts = n >>> 5;
3094         int nBits = n & 0x1f;
3095         int magLen = mag.length;
3096         int newMag[] = null;
3097
3098         if (nBits == 0) {
3099             newMag = new int[magLen + nInts];
3100             System.arraycopy(mag, 0, newMag, 0, magLen);
3101         } else {
3102             int i = 0;
3103             int nBits2 = 32 - nBits;
3104             int highBits = mag[0] >>> nBits2;
3105             if (highBits != 0) {
3106                 newMag = new int[magLen + nInts + 1];
3107                 newMag[i++] = highBits;
3108             } else {
3109                 newMag = new int[magLen + nInts];
3110             }
3111             int j=0;
3112             while (j < magLen-1)
3113                 newMag[i++] = mag[j++] << nBits | mag[j] >>> nBits2;
3114             newMag[i] = mag[j] << nBits;
3115         }
3116         return newMag;
3117     }
3118
3119     /**
3120      * Returns a BigInteger whose value is {@code (this >> n)}.  Sign
3121      * extension is performed.  The shift distance, {@code n}, may be
3122      * negative, in which case this method performs a left shift.
3123      * (Computes <tt>floor(this / 2<sup>n</sup>)</tt>.)
3124      *
3125      * @param  n shift distance, in bits.
3126      * @return {@code this >> n}
3127      * @see #shiftLeft
3128      */

3129     public BigInteger shiftRight(int n) {
3130         if (signum == 0)
3131             return ZERO;
3132         if (n > 0) {
3133             return shiftRightImpl(n);
3134         } else if (n == 0) {
3135             return this;
3136         } else {
3137             // Possible int overflow in {@code -n} is not a trouble,
3138             // because shiftLeft considers its argument unsigned
3139             return new BigInteger(shiftLeft(mag, -n), signum);
3140         }
3141     }
3142
3143     /**
3144      * Returns a BigInteger whose value is {@code (this >> n)}. The shift
3145      * distance, {@code n}, is considered unsigned.
3146      * (Computes <tt>floor(this * 2<sup>-n</sup>)</tt>.)
3147      *
3148      * @param  n unsigned shift distance, in bits.
3149      * @return {@code this >> n}
3150      */

3151     private BigInteger shiftRightImpl(int n) {
3152         int nInts = n >>> 5;
3153         int nBits = n & 0x1f;
3154         int magLen = mag.length;
3155         int newMag[] = null;
3156
3157         // Special case: entire contents shifted off the end
3158         if (nInts >= magLen)
3159             return (signum >= 0 ? ZERO : negConst[1]);
3160
3161         if (nBits == 0) {
3162             int newMagLen = magLen - nInts;
3163             newMag = Arrays.copyOf(mag, newMagLen);
3164         } else {
3165             int i = 0;
3166             int highBits = mag[0] >>> nBits;
3167             if (highBits != 0) {
3168                 newMag = new int[magLen - nInts];
3169                 newMag[i++] = highBits;
3170             } else {
3171                 newMag = new int[magLen - nInts -1];
3172             }
3173
3174             int nBits2 = 32 - nBits;
3175             int j=0;
3176             while (j < magLen - nInts - 1)
3177                 newMag[i++] = (mag[j++] << nBits2) | (mag[j] >>> nBits);
3178         }
3179
3180         if (signum < 0) {
3181             // Find out whether any one-bits were shifted off the end.
3182             boolean onesLost = false;
3183             for (int i=magLen-1, j=magLen-nInts; i >= j && !onesLost; i--)
3184                 onesLost = (mag[i] != 0);
3185             if (!onesLost && nBits != 0)
3186                 onesLost = (mag[magLen - nInts - 1] << (32 - nBits) != 0);
3187
3188             if (onesLost)
3189                 newMag = javaIncrement(newMag);
3190         }
3191
3192         return new BigInteger(newMag, signum);
3193     }
3194
3195     int[] javaIncrement(int[] val) {
3196         int lastSum = 0;
3197         for (int i=val.length-1;  i >= 0 && lastSum == 0; i--)
3198             lastSum = (val[i] += 1);
3199         if (lastSum == 0) {
3200             val = new int[val.length+1];
3201             val[0] = 1;
3202         }
3203         return val;
3204     }
3205
3206     // Bitwise Operations
3207
3208     /**
3209      * Returns a BigInteger whose value is {@code (this & val)}.  (This
3210      * method returns a negative BigInteger if and only if this and val are
3211      * both negative.)
3212      *
3213      * @param val value to be AND'ed with this BigInteger.
3214      * @return {@code this & val}
3215      */

3216     public BigInteger and(BigInteger val) {
3217         int[] result = new int[Math.max(intLength(), val.intLength())];
3218         for (int i=0; i < result.length; i++)
3219             result[i] = (getInt(result.length-i-1)
3220                          & val.getInt(result.length-i-1));
3221
3222         return valueOf(result);
3223     }
3224
3225     /**
3226      * Returns a BigInteger whose value is {@code (this | val)}.  (This method
3227      * returns a negative BigInteger if and only if either this or val is
3228      * negative.)
3229      *
3230      * @param val value to be OR'ed with this BigInteger.
3231      * @return {@code this | val}
3232      */

3233     public BigInteger or(BigInteger val) {
3234         int[] result = new int[Math.max(intLength(), val.intLength())];
3235         for (int i=0; i < result.length; i++)
3236             result[i] = (getInt(result.length-i-1)
3237                          | val.getInt(result.length-i-1));
3238
3239         return valueOf(result);
3240     }
3241
3242     /**
3243      * Returns a BigInteger whose value is {@code (this ^ val)}.  (This method
3244      * returns a negative BigInteger if and only if exactly one of this and
3245      * val are negative.)
3246      *
3247      * @param val value to be XOR'ed with this BigInteger.
3248      * @return {@code this ^ val}
3249      */

3250     public BigInteger xor(BigInteger val) {
3251         int[] result = new int[Math.max(intLength(), val.intLength())];
3252         for (int i=0; i < result.length; i++)
3253             result[i] = (getInt(result.length-i-1)
3254                          ^ val.getInt(result.length-i-1));
3255
3256         return valueOf(result);
3257     }
3258
3259     /**
3260      * Returns a BigInteger whose value is {@code (~this)}.  (This method
3261      * returns a negative value if and only if this BigInteger is
3262      * non-negative.)
3263      *
3264      * @return {@code ~this}
3265      */

3266     public BigInteger not() {
3267         int[] result = new int[intLength()];
3268         for (int i=0; i < result.length; i++)
3269             result[i] = ~getInt(result.length-i-1);
3270
3271         return valueOf(result);
3272     }
3273
3274     /**
3275      * Returns a BigInteger whose value is {@code (this & ~val)}.  This
3276      * method, which is equivalent to {@code and(val.not())}, is provided as
3277      * a convenience for masking operations.  (This method returns a negative
3278      * BigInteger if and only if {@code this} is negative and {@code val} is
3279      * positive.)
3280      *
3281      * @param val value to be complemented and AND'ed with this BigInteger.
3282      * @return {@code this & ~val}
3283      */

3284     public BigInteger andNot(BigInteger val) {
3285         int[] result = new int[Math.max(intLength(), val.intLength())];
3286         for (int i=0; i < result.length; i++)
3287             result[i] = (getInt(result.length-i-1)
3288                          & ~val.getInt(result.length-i-1));
3289
3290         return valueOf(result);
3291     }
3292
3293
3294     // Single Bit Operations
3295
3296     /**
3297      * Returns {@code trueif and only if the designated bit is set.
3298      * (Computes {@code ((this & (1<<n)) != 0)}.)
3299      *
3300      * @param  n index of bit to test.
3301      * @return {@code trueif and only if the designated bit is set.
3302      * @throws ArithmeticException {@code n} is negative.
3303      */

3304     public boolean testBit(int n) {
3305         if (n < 0)
3306             throw new ArithmeticException("Negative bit address");
3307
3308         return (getInt(n >>> 5) & (1 << (n & 31))) != 0;
3309     }
3310
3311     /**
3312      * Returns a BigInteger whose value is equivalent to this BigInteger
3313      * with the designated bit set.  (Computes {@code (this | (1<<n))}.)
3314      *
3315      * @param  n index of bit to set.
3316      * @return {@code this | (1<<n)}
3317      * @throws ArithmeticException {@code n} is negative.
3318      */

3319     public BigInteger setBit(int n) {
3320         if (n < 0)
3321             throw new ArithmeticException("Negative bit address");
3322
3323         int intNum = n >>> 5;
3324         int[] result = new int[Math.max(intLength(), intNum+2)];
3325
3326         for (int i=0; i < result.length; i++)
3327             result[result.length-i-1] = getInt(i);
3328
3329         result[result.length-intNum-1] |= (1 << (n & 31));
3330
3331         return valueOf(result);
3332     }
3333
3334     /**
3335      * Returns a BigInteger whose value is equivalent to this BigInteger
3336      * with the designated bit cleared.
3337      * (Computes {@code (this & ~(1<<n))}.)
3338      *
3339      * @param  n index of bit to clear.
3340      * @return {@code this & ~(1<<n)}
3341      * @throws ArithmeticException {@code n} is negative.
3342      */

3343     public BigInteger clearBit(int n) {
3344         if (n < 0)
3345             throw new ArithmeticException("Negative bit address");
3346
3347         int intNum = n >>> 5;
3348         int[] result = new int[Math.max(intLength(), ((n + 1) >>> 5) + 1)];
3349
3350         for (int i=0; i < result.length; i++)
3351             result[result.length-i-1] = getInt(i);
3352
3353         result[result.length-intNum-1] &= ~(1 << (n & 31));
3354
3355         return valueOf(result);
3356     }
3357
3358     /**
3359      * Returns a BigInteger whose value is equivalent to this BigInteger
3360      * with the designated bit flipped.
3361      * (Computes {@code (this ^ (1<<n))}.)
3362      *
3363      * @param  n index of bit to flip.
3364      * @return {@code this ^ (1<<n)}
3365      * @throws ArithmeticException {@code n} is negative.
3366      */

3367     public BigInteger flipBit(int n) {
3368         if (n < 0)
3369             throw new ArithmeticException("Negative bit address");
3370
3371         int intNum = n >>> 5;
3372         int[] result = new int[Math.max(intLength(), intNum+2)];
3373
3374         for (int i=0; i < result.length; i++)
3375             result[result.length-i-1] = getInt(i);
3376
3377         result[result.length-intNum-1] ^= (1 << (n & 31));
3378
3379         return valueOf(result);
3380     }
3381
3382     /**
3383      * Returns the index of the rightmost (lowest-order) one bit in this
3384      * BigInteger (the number of zero bits to the right of the rightmost
3385      * one bit).  Returns -1 if this BigInteger contains no one bits.
3386      * (Computes {@code (this == 0? -1 : log2(this & -this))}.)
3387      *
3388      * @return index of the rightmost one bit in this BigInteger.
3389      */

3390     public int getLowestSetBit() {
3391         @SuppressWarnings("deprecation"int lsb = lowestSetBit - 2;
3392         if (lsb == -2) {  // lowestSetBit not initialized yet
3393             lsb = 0;
3394             if (signum == 0) {
3395                 lsb -= 1;
3396             } else {
3397                 // Search for lowest order nonzero int
3398                 int i,b;
3399                 for (i=0; (b = getInt(i)) == 0; i++)
3400                     ;
3401                 lsb += (i << 5) + Integer.numberOfTrailingZeros(b);
3402             }
3403             lowestSetBit = lsb + 2;
3404         }
3405         return lsb;
3406     }
3407
3408
3409     // Miscellaneous Bit Operations
3410
3411     /**
3412      * Returns the number of bits in the minimal two's-complement
3413      * representation of this BigInteger, <i>excluding</i> a sign bit.
3414      * For positive BigIntegers, this is equivalent to the number of bits in
3415      * the ordinary binary representation.  (Computes
3416      * {@code (ceil(log2(this < 0 ? -this : this+1)))}.)
3417      *
3418      * @return number of bits in the minimal two's-complement
3419      *         representation of this BigInteger, <i>excluding</i> a sign bit.
3420      */

3421     public int bitLength() {
3422         @SuppressWarnings("deprecation"int n = bitLength - 1;
3423         if (n == -1) { // bitLength not initialized yet
3424             int[] m = mag;
3425             int len = m.length;
3426             if (len == 0) {
3427                 n = 0; // offset by one to initialize
3428             }  else {
3429                 // Calculate the bit length of the magnitude
3430                 int magBitLength = ((len - 1) << 5) + bitLengthForInt(mag[0]);
3431                  if (signum < 0) {
3432                      // Check if magnitude is a power of two
3433                      boolean pow2 = (Integer.bitCount(mag[0]) == 1);
3434                      for (int i=1; i< len && pow2; i++)
3435                          pow2 = (mag[i] == 0);
3436
3437                      n = (pow2 ? magBitLength -1 : magBitLength);
3438                  } else {
3439                      n = magBitLength;
3440                  }
3441             }
3442             bitLength = n + 1;
3443         }
3444         return n;
3445     }
3446
3447     /**
3448      * Returns the number of bits in the two's complement representation
3449      * of this BigInteger that differ from its sign bit.  This method is
3450      * useful when implementing bit-vector style sets atop BigIntegers.
3451      *
3452      * @return number of bits in the two's complement representation
3453      *         of this BigInteger that differ from its sign bit.
3454      */

3455     public int bitCount() {
3456         @SuppressWarnings("deprecation"int bc = bitCount - 1;
3457         if (bc == -1) {  // bitCount not initialized yet
3458             bc = 0;      // offset by one to initialize
3459             // Count the bits in the magnitude
3460             for (int i=0; i < mag.length; i++)
3461                 bc += Integer.bitCount(mag[i]);
3462             if (signum < 0) {
3463                 // Count the trailing zeros in the magnitude
3464                 int magTrailingZeroCount = 0, j;
3465                 for (j=mag.length-1; mag[j] == 0; j--)
3466                     magTrailingZeroCount += 32;
3467                 magTrailingZeroCount += Integer.numberOfTrailingZeros(mag[j]);
3468                 bc += magTrailingZeroCount - 1;
3469             }
3470             bitCount = bc + 1;
3471         }
3472         return bc;
3473     }
3474
3475     // Primality Testing
3476
3477     /**
3478      * Returns {@code trueif this BigInteger is probably prime,
3479      * {@code falseif it's definitely composite.  If
3480      * {@code certainty} is &le; 0, {@code true} is
3481      * returned.
3482      *
3483      * @param  certainty a measure of the uncertainty that the caller is
3484      *         willing to tolerate: if the call returns {@code true}
3485      *         the probability that this BigInteger is prime exceeds
3486      *         (1 - 1/2<sup>{@code certainty}</sup>).  The execution time of
3487      *         this method is proportional to the value of this parameter.
3488      * @return {@code trueif this BigInteger is probably prime,
3489      *         {@code falseif it's definitely composite.
3490      */

3491     public boolean isProbablePrime(int certainty) {
3492         if (certainty <= 0)
3493             return true;
3494         BigInteger w = this.abs();
3495         if (w.equals(TWO))
3496             return true;
3497         if (!w.testBit(0) || w.equals(ONE))
3498             return false;
3499
3500         return w.primeToCertainty(certainty, null);
3501     }
3502
3503     // Comparison Operations
3504
3505     /**
3506      * Compares this BigInteger with the specified BigInteger.  This
3507      * method is provided in preference to individual methods for each
3508      * of the six boolean comparison operators ({@literal <}, ==,
3509      * {@literal >}, {@literal >=}, !=, {@literal <=}).  The suggested
3510      * idiom for performing these comparisons is: {@code
3511      * (x.compareTo(y)} &lt;<i>op</i>&gt; {@code 0)}, where
3512      * &lt;<i>op</i>&gt; is one of the six comparison operators.
3513      *
3514      * @param  val BigInteger to which this BigInteger is to be compared.
3515      * @return -1, 0 or 1 as this BigInteger is numerically less than, equal
3516      *         to, or greater than {@code val}.
3517      */

3518     public int compareTo(BigInteger val) {
3519         if (signum == val.signum) {
3520             switch (signum) {
3521             case 1:
3522                 return compareMagnitude(val);
3523             case -1:
3524                 return val.compareMagnitude(this);
3525             default:
3526                 return 0;
3527             }
3528         }
3529         return signum > val.signum ? 1 : -1;
3530     }
3531
3532     /**
3533      * Compares the magnitude array of this BigInteger with the specified
3534      * BigInteger's. This is the version of compareTo ignoring sign.
3535      *
3536      * @param val BigInteger whose magnitude array to be compared.
3537      * @return -1, 0 or 1 as this magnitude array is less than, equal to or
3538      *         greater than the magnitude aray for the specified BigInteger's.
3539      */

3540     final int compareMagnitude(BigInteger val) {
3541         int[] m1 = mag;
3542         int len1 = m1.length;
3543         int[] m2 = val.mag;
3544         int len2 = m2.length;
3545         if (len1 < len2)
3546             return -1;
3547         if (len1 > len2)
3548             return 1;
3549         for (int i = 0; i < len1; i++) {
3550             int a = m1[i];
3551             int b = m2[i];
3552             if (a != b)
3553                 return ((a & LONG_MASK) < (b & LONG_MASK)) ? -1 : 1;
3554         }
3555         return 0;
3556     }
3557
3558     /**
3559      * Version of compareMagnitude that compares magnitude with long value.
3560      * val can't be Long.MIN_VALUE.
3561      */

3562     final int compareMagnitude(long val) {
3563         assert val != Long.MIN_VALUE;
3564         int[] m1 = mag;
3565         int len = m1.length;
3566         if (len > 2) {
3567             return 1;
3568         }
3569         if (val < 0) {
3570             val = -val;
3571         }
3572         int highWord = (int)(val >>> 32);
3573         if (highWord == 0) {
3574             if (len < 1)
3575                 return -1;
3576             if (len > 1)
3577                 return 1;
3578             int a = m1[0];
3579             int b = (int)val;
3580             if (a != b) {
3581                 return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
3582             }
3583             return 0;
3584         } else {
3585             if (len < 2)
3586                 return -1;
3587             int a = m1[0];
3588             int b = highWord;
3589             if (a != b) {
3590                 return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
3591             }
3592             a = m1[1];
3593             b = (int)val;
3594             if (a != b) {
3595                 return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
3596             }
3597             return 0;
3598         }
3599     }
3600
3601     /**
3602      * Compares this BigInteger with the specified Object for equality.
3603      *
3604      * @param  x Object to which this BigInteger is to be compared.
3605      * @return {@code trueif and only if the specified Object is a
3606      *         BigInteger whose value is numerically equal to this BigInteger.
3607      */

3608     public boolean equals(Object x) {
3609         // This test is just an optimization, which may or may not help
3610         if (x == this)
3611             return true;
3612
3613         if (!(x instanceof BigInteger))
3614             return false;
3615
3616         BigInteger xInt = (BigInteger) x;
3617         if (xInt.signum != signum)
3618             return false;
3619
3620         int[] m = mag;
3621         int len = m.length;
3622         int[] xm = xInt.mag;
3623         if (len != xm.length)
3624             return false;
3625
3626         for (int i = 0; i < len; i++)
3627             if (xm[i] != m[i])
3628                 return false;
3629
3630         return true;
3631     }
3632
3633     /**
3634      * Returns the minimum of this BigInteger and {@code val}.
3635      *
3636      * @param  val value with which the minimum is to be computed.
3637      * @return the BigInteger whose value is the lesser of this BigInteger and
3638      *         {@code val}.  If they are equal, either may be returned.
3639      */

3640     public BigInteger min(BigInteger val) {
3641         return (compareTo(val) < 0 ? this : val);
3642     }
3643
3644     /**
3645      * Returns the maximum of this BigInteger and {@code val}.
3646      *
3647      * @param  val value with which the maximum is to be computed.
3648      * @return the BigInteger whose value is the greater of this and
3649      *         {@code val}.  If they are equal, either may be returned.
3650      */

3651     public BigInteger max(BigInteger val) {
3652         return (compareTo(val) > 0 ? this : val);
3653     }
3654
3655
3656     // Hash Function
3657
3658     /**
3659      * Returns the hash code for this BigInteger.
3660      *
3661      * @return hash code for this BigInteger.
3662      */

3663     public int hashCode() {
3664         int hashCode = 0;
3665
3666         for (int i=0; i < mag.length; i++)
3667             hashCode = (int)(31*hashCode + (mag[i] & LONG_MASK));
3668
3669         return hashCode * signum;
3670     }
3671
3672     /**
3673      * Returns the String representation of this BigInteger in the
3674      * given radix.  If the radix is outside the range from {@link
3675      * Character#MIN_RADIX} to {@link Character#MAX_RADIX} inclusive,
3676      * it will default to 10 (as is the case for
3677      * {@code Integer.toString}).  The digit-to-character mapping
3678      * provided by {@code Character.forDigit} is used, and a minus
3679      * sign is prepended if appropriate.  (This representation is
3680      * compatible with the {@link #BigInteger(String, int) (String,
3681      * int)} constructor.)
3682      *
3683      * @param  radix  radix of the String representation.
3684      * @return String representation of this BigInteger in the given radix.
3685      * @see    Integer#toString
3686      * @see    Character#forDigit
3687      * @see    #BigInteger(java.lang.String, int)
3688      */

3689     public String toString(int radix) {
3690         if (signum == 0)
3691             return "0";
3692         if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
3693             radix = 10;
3694
3695         // If it's small enough, use smallToString.
3696         if (mag.length <= SCHOENHAGE_BASE_CONVERSION_THRESHOLD)
3697            return smallToString(radix);
3698
3699         // Otherwise use recursive toString, which requires positive arguments.
3700         // The results will be concatenated into this StringBuilder
3701         StringBuilder sb = new StringBuilder();
3702         if (signum < 0) {
3703             toString(this.negate(), sb, radix, 0);
3704             sb.insert(0, '-');
3705         }
3706         else
3707             toString(this, sb, radix, 0);
3708
3709         return sb.toString();
3710     }
3711
3712     /** This method is used to perform toString when arguments are small. */
3713     private String smallToString(int radix) {
3714         if (signum == 0) {
3715             return "0";
3716         }
3717
3718         // Compute upper bound on number of digit groups and allocate space
3719         int maxNumDigitGroups = (4*mag.length + 6)/7;
3720         String digitGroup[] = new String[maxNumDigitGroups];
3721
3722         // Translate number to string, a digit group at a time
3723         BigInteger tmp = this.abs();
3724         int numGroups = 0;
3725         while (tmp.signum != 0) {
3726             BigInteger d = longRadix[radix];
3727
3728             MutableBigInteger q = new MutableBigInteger(),
3729                               a = new MutableBigInteger(tmp.mag),
3730                               b = new MutableBigInteger(d.mag);
3731             MutableBigInteger r = a.divide(b, q);
3732             BigInteger q2 = q.toBigInteger(tmp.signum * d.signum);
3733             BigInteger r2 = r.toBigInteger(tmp.signum * d.signum);
3734
3735             digitGroup[numGroups++] = Long.toString(r2.longValue(), radix);
3736             tmp = q2;
3737         }
3738
3739         // Put sign (if any) and first digit group into result buffer
3740         StringBuilder buf = new StringBuilder(numGroups*digitsPerLong[radix]+1);
3741         if (signum < 0) {
3742             buf.append('-');
3743         }
3744         buf.append(digitGroup[numGroups-1]);
3745
3746         // Append remaining digit groups padded with leading zeros
3747         for (int i=numGroups-2; i >= 0; i--) {
3748             // Prepend (any) leading zeros for this digit group
3749             int numLeadingZeros = digitsPerLong[radix]-digitGroup[i].length();
3750             if (numLeadingZeros != 0) {
3751                 buf.append(zeros[numLeadingZeros]);
3752             }
3753             buf.append(digitGroup[i]);
3754         }
3755         return buf.toString();
3756     }
3757
3758     /**
3759      * Converts the specified BigInteger to a string and appends to
3760      * {@code sb}.  This implements the recursive Schoenhage algorithm
3761      * for base conversions.
3762      * <p/>
3763      * See Knuth, Donald,  _The Art of Computer Programming_, Vol. 2,
3764      * Answers to Exercises (4.4) Question 14.
3765      *
3766      * @param u      The number to convert to a string.
3767      * @param sb     The StringBuilder that will be appended to in place.
3768      * @param radix  The base to convert to.
3769      * @param digits The minimum number of digits to pad to.
3770      */

3771     private static void toString(BigInteger u, StringBuilder sb, int radix,
3772                                  int digits) {
3773         /* If we're smaller than a certain threshold, use the smallToString
3774            method, padding with leading zeroes when necessary. */

3775         if (u.mag.length <= SCHOENHAGE_BASE_CONVERSION_THRESHOLD) {
3776             String s = u.smallToString(radix);
3777
3778             // Pad with internal zeros if necessary.
3779             // Don't pad if we're at the beginning of the string.
3780             if ((s.length() < digits) && (sb.length() > 0)) {
3781                 for (int i=s.length(); i < digits; i++) { // May be a faster way to
3782                     sb.append('0');                    // do this?
3783                 }
3784             }
3785
3786             sb.append(s);
3787             return;
3788         }
3789
3790         int b, n;
3791         b = u.bitLength();
3792
3793         // Calculate a value for n in the equation radix^(2^n) = u
3794         // and subtract 1 from that value.  This is used to find the
3795         // cache index that contains the best value to divide u.
3796         n = (int) Math.round(Math.log(b * LOG_TWO / logCache[radix]) / LOG_TWO - 1.0);
3797         BigInteger v = getRadixConversionCache(radix, n);
3798         BigInteger[] results;
3799         results = u.divideAndRemainder(v);
3800
3801         int expectedDigits = 1 << n;
3802
3803         // Now recursively build the two halves of each number.
3804         toString(results[0], sb, radix, digits-expectedDigits);
3805         toString(results[1], sb, radix, expectedDigits);
3806     }
3807
3808     /**
3809      * Returns the value radix^(2^exponent) from the cache.
3810      * If this value doesn't already exist in the cache, it is added.
3811      * <p/>
3812      * This could be changed to a more complicated caching method using
3813      * {@code Future}.
3814      */

3815     private static BigInteger getRadixConversionCache(int radix, int exponent) {
3816         BigInteger[] cacheLine = powerCache[radix]; // volatile read
3817         if (exponent < cacheLine.length) {
3818             return cacheLine[exponent];
3819         }
3820
3821         int oldLength = cacheLine.length;
3822         cacheLine = Arrays.copyOf(cacheLine, exponent + 1);
3823         for (int i = oldLength; i <= exponent; i++) {
3824             cacheLine[i] = cacheLine[i - 1].pow(2);
3825         }
3826
3827         BigInteger[][] pc = powerCache; // volatile read again
3828         if (exponent >= pc[radix].length) {
3829             pc = pc.clone();
3830             pc[radix] = cacheLine;
3831             powerCache = pc; // volatile write, publish
3832         }
3833         return cacheLine[exponent];
3834     }
3835
3836     /* zero[i] is a string of i consecutive zeros. */
3837     private static String zeros[] = new String[64];
3838     static {
3839         zeros[63] =
3840             "000000000000000000000000000000000000000000000000000000000000000";
3841         for (int i=0; i < 63; i++)
3842             zeros[i] = zeros[63].substring(0, i);
3843     }
3844
3845     /**
3846      * Returns the decimal String representation of this BigInteger.
3847      * The digit-to-character mapping provided by
3848      * {@code Character.forDigit} is used, and a minus sign is
3849      * prepended if appropriate.  (This representation is compatible
3850      * with the {@link #BigInteger(String) (String)} constructor, and
3851      * allows for String concatenation with Java's + operator.)
3852      *
3853      * @return decimal String representation of this BigInteger.
3854      * @see    Character#forDigit
3855      * @see    #BigInteger(java.lang.String)
3856      */

3857     public String toString() {
3858         return toString(10);
3859     }
3860
3861     /**
3862      * Returns a byte array containing the two's-complement
3863      * representation of this BigInteger.  The byte array will be in
3864      * <i>big-endian</i> byte-order: the most significant byte is in
3865      * the zeroth element.  The array will contain the minimum number
3866      * of bytes required to represent this BigInteger, including at
3867      * least one sign bit, which is {@code (ceil((this.bitLength() +
3868      * 1)/8))}.  (This representation is compatible with the
3869      * {@link #BigInteger(byte[]) (byte[])} constructor.)
3870      *
3871      * @return a byte array containing the two's-complement representation of
3872      *         this BigInteger.
3873      * @see    #BigInteger(byte[])
3874      */

3875     public byte[] toByteArray() {
3876         int byteLen = bitLength()/8 + 1;
3877         byte[] byteArray = new byte[byteLen];
3878
3879         for (int i=byteLen-1, bytesCopied=4, nextInt=0, intIndex=0; i >= 0; i--) {
3880             if (bytesCopied == 4) {
3881                 nextInt = getInt(intIndex++);
3882                 bytesCopied = 1;
3883             } else {
3884                 nextInt >>>= 8;
3885                 bytesCopied++;
3886             }
3887             byteArray[i] = (byte)nextInt;
3888         }
3889         return byteArray;
3890     }
3891
3892     /**
3893      * Converts this BigInteger to an {@code int}.  This
3894      * conversion is analogous to a
3895      * <i>narrowing primitive conversion</i> from {@code long} to
3896      * {@code int} as defined in section 5.1.3 of
3897      * <cite>The Java&trade; Language Specification</cite>:
3898      * if this BigInteger is too big to fit in an
3899      * {@code int}, only the low-order 32 bits are returned.
3900      * Note that this conversion can lose information about the
3901      * overall magnitude of the BigInteger value as well as return a
3902      * result with the opposite sign.
3903      *
3904      * @return this BigInteger converted to an {@code int}.
3905      * @see #intValueExact()
3906      */

3907     public int intValue() {
3908         int result = 0;
3909         result = getInt(0);
3910         return result;
3911     }
3912
3913     /**
3914      * Converts this BigInteger to a {@code long}.  This
3915      * conversion is analogous to a
3916      * <i>narrowing primitive conversion</i> from {@code long} to
3917      * {@code int} as defined in section 5.1.3 of
3918      * <cite>The Java&trade; Language Specification</cite>:
3919      * if this BigInteger is too big to fit in a
3920      * {@code long}, only the low-order 64 bits are returned.
3921      * Note that this conversion can lose information about the
3922      * overall magnitude of the BigInteger value as well as return a
3923      * result with the opposite sign.
3924      *
3925      * @return this BigInteger converted to a {@code long}.
3926      * @see #longValueExact()
3927      */

3928     public long longValue() {
3929         long result = 0;
3930
3931         for (int i=1; i >= 0; i--)
3932             result = (result << 32) + (getInt(i) & LONG_MASK);
3933         return result;
3934     }
3935
3936     /**
3937      * Converts this BigInteger to a {@code float}.  This
3938      * conversion is similar to the
3939      * <i>narrowing primitive conversion</i> from {@code double} to
3940      * {@code float} as defined in section 5.1.3 of
3941      * <cite>The Java&trade; Language Specification</cite>:
3942      * if this BigInteger has too great a magnitude
3943      * to represent as a {@code float}, it will be converted to
3944      * {@link Float#NEGATIVE_INFINITY} or {@link
3945      * Float#POSITIVE_INFINITY} as appropriate.  Note that even when
3946      * the return value is finite, this conversion can lose
3947      * information about the precision of the BigInteger value.
3948      *
3949      * @return this BigInteger converted to a {@code float}.
3950      */

3951     public float floatValue() {
3952         if (signum == 0) {
3953             return 0.0f;
3954         }
3955
3956         int exponent = ((mag.length - 1) << 5) + bitLengthForInt(mag[0]) - 1;
3957
3958         // exponent == floor(log2(abs(this)))
3959         if (exponent < Long.SIZE - 1) {
3960             return longValue();
3961         } else if (exponent > Float.MAX_EXPONENT) {
3962             return signum > 0 ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3963         }
3964
3965         /*
3966          * We need the top SIGNIFICAND_WIDTH bits, including the "implicit"
3967          * one bit. To make rounding easier, we pick out the top
3968          * SIGNIFICAND_WIDTH + 1 bits, so we have one to help us round up or
3969          * down. twiceSignifFloor will contain the top SIGNIFICAND_WIDTH + 1
3970          * bits, and signifFloor the top SIGNIFICAND_WIDTH.
3971          *
3972          * It helps to consider the real number signif = abs(this) *
3973          * 2^(SIGNIFICAND_WIDTH - 1 - exponent).
3974          */

3975         int shift = exponent - FloatConsts.SIGNIFICAND_WIDTH;
3976
3977         int twiceSignifFloor;
3978         // twiceSignifFloor will be == abs().shiftRight(shift).intValue()
3979         // We do the shift into an int directly to improve performance.
3980
3981         int nBits = shift & 0x1f;
3982         int nBits2 = 32 - nBits;
3983
3984         if (nBits == 0) {
3985             twiceSignifFloor = mag[0];
3986         } else {
3987             twiceSignifFloor = mag[0] >>> nBits;
3988             if (twiceSignifFloor == 0) {
3989                 twiceSignifFloor = (mag[0] << nBits2) | (mag[1] >>> nBits);
3990             }
3991         }
3992
3993         int signifFloor = twiceSignifFloor >> 1;
3994         signifFloor &= FloatConsts.SIGNIF_BIT_MASK; // remove the implied bit
3995
3996         /*
3997          * We round up if either the fractional part of signif is strictly
3998          * greater than 0.5 (which is true if the 0.5 bit is set and any lower
3999          * bit is set), or if the fractional part of signif is >= 0.5 and
4000          * signifFloor is odd (which is true if both the 0.5 bit and the 1 bit
4001          * are set). This is equivalent to the desired HALF_EVEN rounding.
4002          */

4003         boolean increment = (twiceSignifFloor & 1) != 0
4004                 && ((signifFloor & 1) != 0 || abs().getLowestSetBit() < shift);
4005         int signifRounded = increment ? signifFloor + 1 : signifFloor;
4006         int bits = ((exponent + FloatConsts.EXP_BIAS))
4007                 << (FloatConsts.SIGNIFICAND_WIDTH - 1);
4008         bits += signifRounded;
4009         /*
4010          * If signifRounded == 2^24, we'd need to set all of the significand
4011          * bits to zero and add 1 to the exponent. This is exactly the behavior
4012          * we get from just adding signifRounded to bits directly. If the
4013          * exponent is Float.MAX_EXPONENT, we round up (correctly) to
4014          * Float.POSITIVE_INFINITY.
4015          */

4016         bits |= signum & FloatConsts.SIGN_BIT_MASK;
4017         return Float.intBitsToFloat(bits);
4018     }
4019
4020     /**
4021      * Converts this BigInteger to a {@code double}.  This
4022      * conversion is similar to the
4023      * <i>narrowing primitive conversion</i> from {@code double} to
4024      * {@code float} as defined in section 5.1.3 of
4025      * <cite>The Java&trade; Language Specification</cite>:
4026      * if this BigInteger has too great a magnitude
4027      * to represent as a {@code double}, it will be converted to
4028      * {@link Double#NEGATIVE_INFINITY} or {@link
4029      * Double#POSITIVE_INFINITY} as appropriate.  Note that even when
4030      * the return value is finite, this conversion can lose
4031      * information about the precision of the BigInteger value.
4032      *
4033      * @return this BigInteger converted to a {@code double}.
4034      */

4035     public double doubleValue() {
4036         if (signum == 0) {
4037             return 0.0;
4038         }
4039
4040         int exponent = ((mag.length - 1) << 5) + bitLengthForInt(mag[0]) - 1;
4041
4042         // exponent == floor(log2(abs(this))Double)
4043         if (exponent < Long.SIZE - 1) {
4044             return longValue();
4045         } else if (exponent > Double.MAX_EXPONENT) {
4046             return signum > 0 ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
4047         }
4048
4049         /*
4050          * We need the top SIGNIFICAND_WIDTH bits, including the "implicit"
4051          * one bit. To make rounding easier, we pick out the top
4052          * SIGNIFICAND_WIDTH + 1 bits, so we have one to help us round up or
4053          * down. twiceSignifFloor will contain the top SIGNIFICAND_WIDTH + 1
4054          * bits, and signifFloor the top SIGNIFICAND_WIDTH.
4055          *
4056          * It helps to consider the real number signif = abs(this) *
4057          * 2^(SIGNIFICAND_WIDTH - 1 - exponent).
4058          */

4059         int shift = exponent - DoubleConsts.SIGNIFICAND_WIDTH;
4060
4061         long twiceSignifFloor;
4062         // twiceSignifFloor will be == abs().shiftRight(shift).longValue()
4063         // We do the shift into a long directly to improve performance.
4064
4065         int nBits = shift & 0x1f;
4066         int nBits2 = 32 - nBits;
4067
4068         int highBits;
4069         int lowBits;
4070         if (nBits == 0) {
4071             highBits = mag[0];
4072             lowBits = mag[1];
4073         } else {
4074             highBits = mag[0] >>> nBits;
4075             lowBits = (mag[0] << nBits2) | (mag[1] >>> nBits);
4076             if (highBits == 0) {
4077                 highBits = lowBits;
4078                 lowBits = (mag[1] << nBits2) | (mag[2] >>> nBits);
4079             }
4080         }
4081
4082         twiceSignifFloor = ((highBits & LONG_MASK) << 32)
4083                 | (lowBits & LONG_MASK);
4084
4085         long signifFloor = twiceSignifFloor >> 1;
4086         signifFloor &= DoubleConsts.SIGNIF_BIT_MASK; // remove the implied bit
4087
4088         /*
4089          * We round up if either the fractional part of signif is strictly
4090          * greater than 0.5 (which is true if the 0.5 bit is set and any lower
4091          * bit is set), or if the fractional part of signif is >= 0.5 and
4092          * signifFloor is odd (which is true if both the 0.5 bit and the 1 bit
4093          * are set). This is equivalent to the desired HALF_EVEN rounding.
4094          */

4095         boolean increment = (twiceSignifFloor & 1) != 0
4096                 && ((signifFloor & 1) != 0 || abs().getLowestSetBit() < shift);
4097         long signifRounded = increment ? signifFloor + 1 : signifFloor;
4098         long bits = (long) ((exponent + DoubleConsts.EXP_BIAS))
4099                 << (DoubleConsts.SIGNIFICAND_WIDTH - 1);
4100         bits += signifRounded;
4101         /*
4102          * If signifRounded == 2^53, we'd need to set all of the significand
4103          * bits to zero and add 1 to the exponent. This is exactly the behavior
4104          * we get from just adding signifRounded to bits directly. If the
4105          * exponent is Double.MAX_EXPONENT, we round up (correctly) to
4106          * Double.POSITIVE_INFINITY.
4107          */

4108         bits |= signum & DoubleConsts.SIGN_BIT_MASK;
4109         return Double.longBitsToDouble(bits);
4110     }
4111
4112     /**
4113      * Returns a copy of the input array stripped of any leading zero bytes.
4114      */

4115     private static int[] stripLeadingZeroInts(int val[]) {
4116         int vlen = val.length;
4117         int keep;
4118
4119         // Find first nonzero byte
4120         for (keep = 0; keep < vlen && val[keep] == 0; keep++)
4121             ;
4122         return java.util.Arrays.copyOfRange(val, keep, vlen);
4123     }
4124
4125     /**
4126      * Returns the input array stripped of any leading zero bytes.
4127      * Since the source is trusted the copying may be skipped.
4128      */

4129     private static int[] trustedStripLeadingZeroInts(int val[]) {
4130         int vlen = val.length;
4131         int keep;
4132
4133         // Find first nonzero byte
4134         for (keep = 0; keep < vlen && val[keep] == 0; keep++)
4135             ;
4136         return keep == 0 ? val : java.util.Arrays.copyOfRange(val, keep, vlen);
4137     }
4138
4139     /**
4140      * Returns a copy of the input array stripped of any leading zero bytes.
4141      */

4142     private static int[] stripLeadingZeroBytes(byte a[]) {
4143         int byteLength = a.length;
4144         int keep;
4145
4146         // Find first nonzero byte
4147         for (keep = 0; keep < byteLength && a[keep] == 0; keep++)
4148             ;
4149
4150         // Allocate new array and copy relevant part of input array
4151         int intLength = ((byteLength - keep) + 3) >>> 2;
4152         int[] result = new int[intLength];
4153         int b = byteLength - 1;
4154         for (int i = intLength-1; i >= 0; i--) {
4155             result[i] = a[b--] & 0xff;
4156             int bytesRemaining = b - keep + 1;
4157             int bytesToTransfer = Math.min(3, bytesRemaining);
4158             for (int j=8; j <= (bytesToTransfer << 3); j += 8)
4159                 result[i] |= ((a[b--] & 0xff) << j);
4160         }
4161         return result;
4162     }
4163
4164     /**
4165      * Takes an array a representing a negative 2's-complement number and
4166      * returns the minimal (no leading zero bytes) unsigned whose value is -a.
4167      */

4168     private static int[] makePositive(byte a[]) {
4169         int keep, k;
4170         int byteLength = a.length;
4171
4172         // Find first non-sign (0xff) byte of input
4173         for (keep=0; keep < byteLength && a[keep] == -1; keep++)
4174             ;
4175
4176
4177         /* Allocate output array.  If all non-sign bytes are 0x00, we must
4178          * allocate space for one extra output byte. */

4179         for (k=keep; k < byteLength && a[k] == 0; k++)
4180             ;
4181
4182         int extraByte = (k == byteLength) ? 1 : 0;
4183         int intLength = ((byteLength - keep + extraByte) + 3) >>> 2;
4184         int result[] = new int[intLength];
4185
4186         /* Copy one's complement of input into output, leaving extra
4187          * byte (if it exists) == 0x00 */

4188         int b = byteLength - 1;
4189         for (int i = intLength-1; i >= 0; i--) {
4190             result[i] = a[b--] & 0xff;
4191             int numBytesToTransfer = Math.min(3, b-keep+1);
4192             if (numBytesToTransfer < 0)
4193                 numBytesToTransfer = 0;
4194             for (int j=8; j <= 8*numBytesToTransfer; j += 8)
4195                 result[i] |= ((a[b--] & 0xff) << j);
4196
4197             // Mask indicates which bits must be complemented
4198             int mask = -1 >>> (8*(3-numBytesToTransfer));
4199             result[i] = ~result[i] & mask;
4200         }
4201
4202         // Add one to one's complement to generate two's complement
4203         for (int i=result.length-1; i >= 0; i--) {
4204             result[i] = (int)((result[i] & LONG_MASK) + 1);
4205             if (result[i] != 0)
4206                 break;
4207         }
4208
4209         return result;
4210     }
4211
4212     /**
4213      * Takes an array a representing a negative 2's-complement number and
4214      * returns the minimal (no leading zero ints) unsigned whose value is -a.
4215      */

4216     private static int[] makePositive(int a[]) {
4217         int keep, j;
4218
4219         // Find first non-sign (0xffffffff) int of input
4220         for (keep=0; keep < a.length && a[keep] == -1; keep++)
4221             ;
4222
4223         /* Allocate output array.  If all non-sign ints are 0x00, we must
4224          * allocate space for one extra output int. */

4225         for (j=keep; j < a.length && a[j] == 0; j++)
4226             ;
4227         int extraInt = (j == a.length ? 1 : 0);
4228         int result[] = new int[a.length - keep + extraInt];
4229
4230         /* Copy one's complement of input into output, leaving extra
4231          * int (if it exists) == 0x00 */

4232         for (int i = keep; i < a.length; i++)
4233             result[i - keep + extraInt] = ~a[i];
4234
4235         // Add one to one's complement to generate two's complement
4236         for (int i=result.length-1; ++result[i] == 0; i--)
4237             ;
4238
4239         return result;
4240     }
4241
4242     /*
4243      * The following two arrays are used for fast String conversions.  Both
4244      * are indexed by radix.  The first is the number of digits of the given
4245      * radix that can fit in a Java long without "going negative", i.e., the
4246      * highest integer n such that radix**n < 2**63.  The second is the
4247      * "long radix" that tears each number into "long digits", each of which
4248      * consists of the number of digits in the corresponding element in
4249      * digitsPerLong (longRadix[i] = i**digitPerLong[i]).  Both arrays have
4250      * nonsense values in their 0 and 1 elements, as radixes 0 and 1 are not
4251      * used.
4252      */

4253     private static int digitsPerLong[] = {0, 0,
4254         62, 39, 31, 27, 24, 22, 20, 19, 18, 18, 17, 17, 16, 16, 15, 15, 15, 14,
4255         14, 14, 14, 13, 13, 13, 13, 13, 13, 12, 12, 12, 12, 12, 12, 12, 12};
4256
4257     private static BigInteger longRadix[] = {nullnull,
4258         valueOf(0x4000000000000000L), valueOf(0x383d9170b85ff80bL),
4259         valueOf(0x4000000000000000L), valueOf(0x6765c793fa10079dL),
4260         valueOf(0x41c21cb8e1000000L), valueOf(0x3642798750226111L),
4261         valueOf(0x1000000000000000L), valueOf(0x12bf307ae81ffd59L),
4262         valueOf( 0xde0b6b3a7640000L), valueOf(0x4d28cb56c33fa539L),
4263         valueOf(0x1eca170c00000000L), valueOf(0x780c7372621bd74dL),
4264         valueOf(0x1e39a5057d810000L), valueOf(0x5b27ac993df97701L),
4265         valueOf(0x1000000000000000L), valueOf(0x27b95e997e21d9f1L),
4266         valueOf(0x5da0e1e53c5c8000L), valueOf( 0xb16a458ef403f19L),
4267         valueOf(0x16bcc41e90000000L), valueOf(0x2d04b7fdd9c0ef49L),
4268         valueOf(0x5658597bcaa24000L), valueOf( 0x6feb266931a75b7L),
4269         valueOf( 0xc29e98000000000L), valueOf(0x14adf4b7320334b9L),
4270         valueOf(0x226ed36478bfa000L), valueOf(0x383d9170b85ff80bL),
4271         valueOf(0x5a3c23e39c000000L), valueOf( 0x4e900abb53e6b71L),
4272         valueOf( 0x7600ec618141000L), valueOf( 0xaee5720ee830681L),
4273         valueOf(0x1000000000000000L), valueOf(0x172588ad4f5f0981L),
4274         valueOf(0x211e44f7d02c1000L), valueOf(0x2ee56725f06e5c71L),
4275         valueOf(0x41c21cb8e1000000L)};
4276
4277     /*
4278      * These two arrays are the integer analogue of above.
4279      */

4280     private static int digitsPerInt[] = {0, 0, 30, 19, 15, 13, 11,
4281         11, 10, 9, 9, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6,
4282         6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5};
4283
4284     private static int intRadix[] = {0, 0,
4285         0x40000000, 0x4546b3db, 0x40000000, 0x48c27395, 0x159fd800,
4286         0x75db9c97, 0x40000000, 0x17179149, 0x3b9aca00, 0xcc6db61,
4287         0x19a10000, 0x309f1021, 0x57f6c100, 0xa2f1b6f,  0x10000000,
4288         0x18754571, 0x247dbc80, 0x3547667b, 0x4c4b4000, 0x6b5a6e1d,
4289         0x6c20a40,  0x8d2d931,  0xb640000,  0xe8d4a51,  0x1269ae40,
4290         0x17179149, 0x1cb91000, 0x23744899, 0x2b73a840, 0x34e63b41,
4291         0x40000000, 0x4cfa3cc1, 0x5c13d840, 0x6d91b519, 0x39aa400
4292     };
4293
4294     /**
4295      * These routines provide access to the two's complement representation
4296      * of BigIntegers.
4297      */

4298
4299     /**
4300      * Returns the length of the two's complement representation in ints,
4301      * including space for at least one sign bit.
4302      */

4303     private int intLength() {
4304         return (bitLength() >>> 5) + 1;
4305     }
4306
4307     /* Returns sign bit */
4308     private int signBit() {
4309         return signum < 0 ? 1 : 0;
4310     }
4311
4312     /* Returns an int of sign bits */
4313     private int signInt() {
4314         return signum < 0 ? -1 : 0;
4315     }
4316
4317     /**
4318      * Returns the specified int of the little-endian two's complement
4319      * representation (int 0 is the least significant).  The int number can
4320      * be arbitrarily high (values are logically preceded by infinitely many
4321      * sign ints).
4322      */

4323     private int getInt(int n) {
4324         if (n < 0)
4325             return 0;
4326         if (n >= mag.length)
4327             return signInt();
4328
4329         int magInt = mag[mag.length-n-1];
4330
4331         return (signum >= 0 ? magInt :
4332                 (n <= firstNonzeroIntNum() ? -magInt : ~magInt));
4333     }
4334
4335     /**
4336      * Returns the index of the int that contains the first nonzero int in the
4337      * little-endian binary representation of the magnitude (int 0 is the
4338      * least significant). If the magnitude is zero, return value is undefined.
4339      */

4340     private int firstNonzeroIntNum() {
4341         int fn = firstNonzeroIntNum - 2;
4342         if (fn == -2) { // firstNonzeroIntNum not initialized yet
4343             fn = 0;
4344
4345             // Search for the first nonzero int
4346             int i;
4347             int mlen = mag.length;
4348             for (i = mlen - 1; i >= 0 && mag[i] == 0; i--)
4349                 ;
4350             fn = mlen - i - 1;
4351             firstNonzeroIntNum = fn + 2; // offset by two to initialize
4352         }
4353         return fn;
4354     }
4355
4356     /** use serialVersionUID from JDK 1.1. for interoperability */
4357     private static final long serialVersionUID = -8287574255936472291L;
4358
4359     /**
4360      * Serializable fields for BigInteger.
4361      *
4362      * @serialField signum  int
4363      *              signum of this BigInteger.
4364      * @serialField magnitude int[]
4365      *              magnitude array of this BigInteger.
4366      * @serialField bitCount  int
4367      *              number of bits in this BigInteger
4368      * @serialField bitLength int
4369      *              the number of bits in the minimal two's-complement
4370      *              representation of this BigInteger
4371      * @serialField lowestSetBit int
4372      *              lowest set bit in the twos complement representation
4373      */

4374     private static final ObjectStreamField[] serialPersistentFields = {
4375         new ObjectStreamField("signum", Integer.TYPE),
4376         new ObjectStreamField("magnitude"byte[].class),
4377         new ObjectStreamField("bitCount", Integer.TYPE),
4378         new ObjectStreamField("bitLength", Integer.TYPE),
4379         new ObjectStreamField("firstNonzeroByteNum", Integer.TYPE),
4380         new ObjectStreamField("lowestSetBit", Integer.TYPE)
4381         };
4382
4383     /**
4384      * Reconstitute the {@code BigInteger} instance from a stream (that is,
4385      * deserialize it). The magnitude is read in as an array of bytes
4386      * for historical reasons, but it is converted to an array of ints
4387      * and the byte array is discarded.
4388      * Note:
4389      * The current convention is to initialize the cache fields, bitCount,
4390      * bitLength and lowestSetBit, to 0 rather than some other marker value.
4391      * Therefore, no explicit action to set these fields needs to be taken in
4392      * readObject because those fields already have a 0 value be default since
4393      * defaultReadObject is not being used.
4394      */

4395     private void readObject(java.io.ObjectInputStream s)
4396         throws java.io.IOException, ClassNotFoundException {
4397         /*
4398          * In order to maintain compatibility with previous serialized forms,
4399          * the magnitude of a BigInteger is serialized as an array of bytes.
4400          * The magnitude field is used as a temporary store for the byte array
4401          * that is deserialized. The cached computation fields should be
4402          * transient but are serialized for compatibility reasons.
4403          */

4404
4405         // prepare to read the alternate persistent fields
4406         ObjectInputStream.GetField fields = s.readFields();
4407
4408         // Read the alternate persistent fields that we care about
4409         int sign = fields.get("signum", -2);
4410         byte[] magnitude = (byte[])fields.get("magnitude"null);
4411
4412         // Validate signum
4413         if (sign < -1 || sign > 1) {
4414             String message = "BigInteger: Invalid signum value";
4415             if (fields.defaulted("signum"))
4416                 message = "BigInteger: Signum not present in stream";
4417             throw new java.io.StreamCorruptedException(message);
4418         }
4419         int[] mag = stripLeadingZeroBytes(magnitude);
4420         if ((mag.length == 0) != (sign == 0)) {
4421             String message = "BigInteger: signum-magnitude mismatch";
4422             if (fields.defaulted("magnitude"))
4423                 message = "BigInteger: Magnitude not present in stream";
4424             throw new java.io.StreamCorruptedException(message);
4425         }
4426
4427         // Commit final fields via Unsafe
4428         UnsafeHolder.putSign(this, sign);
4429
4430         // Calculate mag field from magnitude and discard magnitude
4431         UnsafeHolder.putMag(this, mag);
4432         if (mag.length >= MAX_MAG_LENGTH) {
4433             try {
4434                 checkRange();
4435             } catch (ArithmeticException e) {
4436                 throw new java.io.StreamCorruptedException("BigInteger: Out of the supported range");
4437             }
4438         }
4439     }
4440
4441     // Support for resetting final fields while deserializing
4442     private static class UnsafeHolder {
4443         private static final sun.misc.Unsafe unsafe;
4444         private static final long signumOffset;
4445         private static final long magOffset;
4446         static {
4447             try {
4448                 unsafe = sun.misc.Unsafe.getUnsafe();
4449                 signumOffset = unsafe.objectFieldOffset
4450                     (BigInteger.class.getDeclaredField("signum"));
4451                 magOffset = unsafe.objectFieldOffset
4452                     (BigInteger.class.getDeclaredField("mag"));
4453             } catch (Exception ex) {
4454                 throw new ExceptionInInitializerError(ex);
4455             }
4456         }
4457
4458         static void putSign(BigInteger bi, int sign) {
4459             unsafe.putIntVolatile(bi, signumOffset, sign);
4460         }
4461
4462         static void putMag(BigInteger bi, int[] magnitude) {
4463             unsafe.putObjectVolatile(bi, magOffset, magnitude);
4464         }
4465     }
4466
4467     /**
4468      * Save the {@code BigInteger} instance to a stream.
4469      * The magnitude of a BigInteger is serialized as a byte array for
4470      * historical reasons.
4471      *
4472      * @serialData two necessary fields are written as well as obsolete
4473      *             fields for compatibility with older versions.
4474      */

4475     private void writeObject(ObjectOutputStream s) throws IOException {
4476         // set the values of the Serializable fields
4477         ObjectOutputStream.PutField fields = s.putFields();
4478         fields.put("signum", signum);
4479         fields.put("magnitude", magSerializedForm());
4480         // The values written for cached fields are compatible with older
4481         // versions, but are ignored in readObject so don't otherwise matter.
4482         fields.put("bitCount", -1);
4483         fields.put("bitLength", -1);
4484         fields.put("lowestSetBit", -2);
4485         fields.put("firstNonzeroByteNum", -2);
4486
4487         // save them
4488         s.writeFields();
4489 }
4490
4491     /**
4492      * Returns the mag array as an array of bytes.
4493      */

4494     private byte[] magSerializedForm() {
4495         int len = mag.length;
4496
4497         int bitLen = (len == 0 ? 0 : ((len - 1) << 5) + bitLengthForInt(mag[0]));
4498         int byteLen = (bitLen + 7) >>> 3;
4499         byte[] result = new byte[byteLen];
4500
4501         for (int i = byteLen - 1, bytesCopied = 4, intIndex = len - 1, nextInt = 0;
4502              i >= 0; i--) {
4503             if (bytesCopied == 4) {
4504                 nextInt = mag[intIndex--];
4505                 bytesCopied = 1;
4506             } else {
4507                 nextInt >>>= 8;
4508                 bytesCopied++;
4509             }
4510             result[i] = (byte)nextInt;
4511         }
4512         return result;
4513     }
4514
4515     /**
4516      * Converts this {@code BigInteger} to a {@code long}, checking
4517      * for lost information.  If the value of this {@code BigInteger}
4518      * is out of the range of the {@code long} type, then an
4519      * {@code ArithmeticException} is thrown.
4520      *
4521      * @return this {@code BigInteger} converted to a {@code long}.
4522      * @throws ArithmeticException if the value of {@code this} will
4523      * not exactly fit in a {@code long}.
4524      * @see BigInteger#longValue
4525      * @since  1.8
4526      */

4527     public long longValueExact() {
4528         if (mag.length <= 2 && bitLength() <= 63)
4529             return longValue();
4530         else
4531             throw new ArithmeticException("BigInteger out of long range");
4532     }
4533
4534     /**
4535      * Converts this {@code BigInteger} to an {@code int}, checking
4536      * for lost information.  If the value of this {@code BigInteger}
4537      * is out of the range of the {@code int} type, then an
4538      * {@code ArithmeticException} is thrown.
4539      *
4540      * @return this {@code BigInteger} converted to an {@code int}.
4541      * @throws ArithmeticException if the value of {@code this} will
4542      * not exactly fit in a {@code int}.
4543      * @see BigInteger#intValue
4544      * @since  1.8
4545      */

4546     public int intValueExact() {
4547         if (mag.length <= 1 && bitLength() <= 31)
4548             return intValue();
4549         else
4550             throw new ArithmeticException("BigInteger out of int range");
4551     }
4552
4553     /**
4554      * Converts this {@code BigInteger} to a {@code short}, checking
4555      * for lost information.  If the value of this {@code BigInteger}
4556      * is out of the range of the {@code short} type, then an
4557      * {@code ArithmeticException} is thrown.
4558      *
4559      * @return this {@code BigInteger} converted to a {@code short}.
4560      * @throws ArithmeticException if the value of {@code this} will
4561      * not exactly fit in a {@code short}.
4562      * @see BigInteger#shortValue
4563      * @since  1.8
4564      */

4565     public short shortValueExact() {
4566         if (mag.length <= 1 && bitLength() <= 31) {
4567             int value = intValue();
4568             if (value >= Short.MIN_VALUE && value <= Short.MAX_VALUE)
4569                 return shortValue();
4570         }
4571         throw new ArithmeticException("BigInteger out of short range");
4572     }
4573
4574     /**
4575      * Converts this {@code BigInteger} to a {@code byte}, checking
4576      * for lost information.  If the value of this {@code BigInteger}
4577      * is out of the range of the {@code byte} type, then an
4578      * {@code ArithmeticException} is thrown.
4579      *
4580      * @return this {@code BigInteger} converted to a {@code byte}.
4581      * @throws ArithmeticException if the value of {@code this} will
4582      * not exactly fit in a {@code byte}.
4583      * @see BigInteger#byteValue
4584      * @since  1.8
4585      */

4586     public byte byteValueExact() {
4587         if (mag.length <= 1 && bitLength() <= 31) {
4588             int value = intValue();
4589             if (value >= Byte.MIN_VALUE && value <= Byte.MAX_VALUE)
4590                 return byteValue();
4591         }
4592         throw new ArithmeticException("BigInteger out of byte range");
4593     }
4594 }
4595
Powered by JavaMelody