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