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