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