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