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