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