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