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