1 /* 2 * Copyright (c) 1999, 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 package java.math; 27 28 /** 29 * A class used to represent multiprecision integers that makes efficient 30 * use of allocated space by allowing a number to occupy only part of 31 * an array so that the arrays do not have to be reallocated as often. 32 * When performing an operation with many iterations the array used to 33 * hold a number is only reallocated when necessary and does not have to 34 * be the same size as the number it represents. A mutable number allows 35 * calculations to occur on the same number without having to create 36 * a new number for every step of the calculation as occurs with 37 * BigIntegers. 38 * 39 * @see BigInteger 40 * @author Michael McCloskey 41 * @author Timothy Buktu 42 * @since 1.3 43 */ 44 45 import static java.math.BigDecimal.INFLATED; 46 import static java.math.BigInteger.LONG_MASK; 47 import java.util.Arrays; 48 49 class MutableBigInteger { 50 /** 51 * Holds the magnitude of this MutableBigInteger in big endian order. 52 * The magnitude may start at an offset into the value array, and it may 53 * end before the length of the value array. 54 */ 55 int[] value; 56 57 /** 58 * The number of ints of the value array that are currently used 59 * to hold the magnitude of this MutableBigInteger. The magnitude starts 60 * at an offset and offset + intLen may be less than value.length. 61 */ 62 int intLen; 63 64 /** 65 * The offset into the value array where the magnitude of this 66 * MutableBigInteger begins. 67 */ 68 int offset = 0; 69 70 // Constants 71 /** 72 * MutableBigInteger with one element value array with the value 1. Used by 73 * BigDecimal divideAndRound to increment the quotient. Use this constant 74 * only when the method is not going to modify this object. 75 */ 76 static final MutableBigInteger ONE = new MutableBigInteger(1); 77 78 /** 79 * The minimum {@code intLen} for cancelling powers of two before 80 * dividing. 81 * If the number of ints is less than this threshold, 82 * {@code divideKnuth} does not eliminate common powers of two from 83 * the dividend and divisor. 84 */ 85 static final int KNUTH_POW2_THRESH_LEN = 6; 86 87 /** 88 * The minimum number of trailing zero ints for cancelling powers of two 89 * before dividing. 90 * If the dividend and divisor don't share at least this many zero ints 91 * at the end, {@code divideKnuth} does not eliminate common powers 92 * of two from the dividend and divisor. 93 */ 94 static final int KNUTH_POW2_THRESH_ZEROS = 3; 95 96 // Constructors 97 98 /** 99 * The default constructor. An empty MutableBigInteger is created with 100 * a one word capacity. 101 */ 102 MutableBigInteger() { 103 value = new int[1]; 104 intLen = 0; 105 } 106 107 /** 108 * Construct a new MutableBigInteger with a magnitude specified by 109 * the int val. 110 */ 111 MutableBigInteger(int val) { 112 value = new int[1]; 113 intLen = 1; 114 value[0] = val; 115 } 116 117 /** 118 * Construct a new MutableBigInteger with the specified value array 119 * up to the length of the array supplied. 120 */ 121 MutableBigInteger(int[] val) { 122 value = val; 123 intLen = val.length; 124 } 125 126 /** 127 * Construct a new MutableBigInteger with a magnitude equal to the 128 * specified BigInteger. 129 */ 130 MutableBigInteger(BigInteger b) { 131 intLen = b.mag.length; 132 value = Arrays.copyOf(b.mag, intLen); 133 } 134 135 /** 136 * Construct a new MutableBigInteger with a magnitude equal to the 137 * specified MutableBigInteger. 138 */ 139 MutableBigInteger(MutableBigInteger val) { 140 intLen = val.intLen; 141 value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen); 142 } 143 144 /** 145 * Makes this number an {@code n}-int number all of whose bits are ones. 146 * Used by Burnikel-Ziegler division. 147 * @param n number of ints in the {@code value} array 148 * @return a number equal to {@code ((1<<(32*n)))-1} 149 */ 150 private void ones(int n) { 151 if (n > value.length) 152 value = new int[n]; 153 Arrays.fill(value, -1); 154 offset = 0; 155 intLen = n; 156 } 157 158 /** 159 * Internal helper method to return the magnitude array. The caller is not 160 * supposed to modify the returned array. 161 */ 162 private int[] getMagnitudeArray() { 163 if (offset > 0 || value.length != intLen) 164 return Arrays.copyOfRange(value, offset, offset + intLen); 165 return value; 166 } 167 168 /** 169 * Convert this MutableBigInteger to a long value. The caller has to make 170 * sure this MutableBigInteger can be fit into long. 171 */ 172 private long toLong() { 173 assert (intLen <= 2) : "this MutableBigInteger exceeds the range of long"; 174 if (intLen == 0) 175 return 0; 176 long d = value[offset] & LONG_MASK; 177 return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d; 178 } 179 180 /** 181 * Convert this MutableBigInteger to a BigInteger object. 182 */ 183 BigInteger toBigInteger(int sign) { 184 if (intLen == 0 || sign == 0) 185 return BigInteger.ZERO; 186 return new BigInteger(getMagnitudeArray(), sign); 187 } 188 189 /** 190 * Converts this number to a nonnegative {@code BigInteger}. 191 */ 192 BigInteger toBigInteger() { 193 normalize(); 194 return toBigInteger(isZero() ? 0 : 1); 195 } 196 197 /** 198 * Convert this MutableBigInteger to BigDecimal object with the specified sign 199 * and scale. 200 */ 201 BigDecimal toBigDecimal(int sign, int scale) { 202 if (intLen == 0 || sign == 0) 203 return BigDecimal.zeroValueOf(scale); 204 int[] mag = getMagnitudeArray(); 205 int len = mag.length; 206 int d = mag[0]; 207 // If this MutableBigInteger can't be fit into long, we need to 208 // make a BigInteger object for the resultant BigDecimal object. 209 if (len > 2 || (d < 0 && len == 2)) 210 return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0); 211 long v = (len == 2) ? 212 ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) : 213 d & LONG_MASK; 214 return BigDecimal.valueOf(sign == -1 ? -v : v, scale); 215 } 216 217 /** 218 * This is for internal use in converting from a MutableBigInteger 219 * object into a long value given a specified sign. 220 * returns INFLATED if value is not fit into long 221 */ 222 long toCompactValue(int sign) { 223 if (intLen == 0 || sign == 0) 224 return 0L; 225 int[] mag = getMagnitudeArray(); 226 int len = mag.length; 227 int d = mag[0]; 228 // If this MutableBigInteger can not be fitted into long, we need to 229 // make a BigInteger object for the resultant BigDecimal object. 230 if (len > 2 || (d < 0 && len == 2)) 231 return INFLATED; 232 long v = (len == 2) ? 233 ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) : 234 d & LONG_MASK; 235 return sign == -1 ? -v : v; 236 } 237 238 /** 239 * Clear out a MutableBigInteger for reuse. 240 */ 241 void clear() { 242 offset = intLen = 0; 243 for (int index=0, n=value.length; index < n; index++) 244 value[index] = 0; 245 } 246 247 /** 248 * Set a MutableBigInteger to zero, removing its offset. 249 */ 250 void reset() { 251 offset = intLen = 0; 252 } 253 254 /** 255 * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1 256 * as this MutableBigInteger is numerically less than, equal to, or 257 * greater than <tt>b</tt>. 258 */ 259 final int compare(MutableBigInteger b) { 260 int blen = b.intLen; 261 if (intLen < blen) 262 return -1; 263 if (intLen > blen) 264 return 1; 265 266 // Add Integer.MIN_VALUE to make the comparison act as unsigned integer 267 // comparison. 268 int[] bval = b.value; 269 for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) { 270 int b1 = value[i] + 0x80000000; 271 int b2 = bval[j] + 0x80000000; 272 if (b1 < b2) 273 return -1; 274 if (b1 > b2) 275 return 1; 276 } 277 return 0; 278 } 279 280 /** 281 * Returns a value equal to what {@code b.leftShift(32*ints); return compare(b);} 282 * would return, but doesn't change the value of {@code b}. 283 */ 284 private int compareShifted(MutableBigInteger b, int ints) { 285 int blen = b.intLen; 286 int alen = intLen - ints; 287 if (alen < blen) 288 return -1; 289 if (alen > blen) 290 return 1; 291 292 // Add Integer.MIN_VALUE to make the comparison act as unsigned integer 293 // comparison. 294 int[] bval = b.value; 295 for (int i = offset, j = b.offset; i < alen + offset; i++, j++) { 296 int b1 = value[i] + 0x80000000; 297 int b2 = bval[j] + 0x80000000; 298 if (b1 < b2) 299 return -1; 300 if (b1 > b2) 301 return 1; 302 } 303 return 0; 304 } 305 306 /** 307 * Compare this against half of a MutableBigInteger object (Needed for 308 * remainder tests). 309 * Assumes no leading unnecessary zeros, which holds for results 310 * from divide(). 311 */ 312 final int compareHalf(MutableBigInteger b) { 313 int blen = b.intLen; 314 int len = intLen; 315 if (len <= 0) 316 return blen <=0 ? 0 : -1; 317 if (len > blen) 318 return 1; 319 if (len < blen - 1) 320 return -1; 321 int[] bval = b.value; 322 int bstart = 0; 323 int carry = 0; 324 // Only 2 cases left:len == blen or len == blen - 1 325 if (len != blen) { // len == blen - 1 326 if (bval[bstart] == 1) { 327 ++bstart; 328 carry = 0x80000000; 329 } else 330 return -1; 331 } 332 // compare values with right-shifted values of b, 333 // carrying shifted-out bits across words 334 int[] val = value; 335 for (int i = offset, j = bstart; i < len + offset;) { 336 int bv = bval[j++]; 337 long hb = ((bv >>> 1) + carry) & LONG_MASK; 338 long v = val[i++] & LONG_MASK; 339 if (v != hb) 340 return v < hb ? -1 : 1; 341 carry = (bv & 1) << 31; // carray will be either 0x80000000 or 0 342 } 343 return carry == 0? 0 : -1; 344 } 345 346 /** 347 * Return the index of the lowest set bit in this MutableBigInteger. If the 348 * magnitude of this MutableBigInteger is zero, -1 is returned. 349 */ 350 private final int getLowestSetBit() { 351 if (intLen == 0) 352 return -1; 353 int j, b; 354 for (j=intLen-1; (j>0) && (value[j+offset]==0); j--) 355 ; 356 b = value[j+offset]; 357 if (b==0) 358 return -1; 359 return ((intLen-1-j)<<5) + Integer.numberOfTrailingZeros(b); 360 } 361 362 /** 363 * Return the int in use in this MutableBigInteger at the specified 364 * index. This method is not used because it is not inlined on all 365 * platforms. 366 */ 367 private final int getInt(int index) { 368 return value[offset+index]; 369 } 370 371 /** 372 * Return a long which is equal to the unsigned value of the int in 373 * use in this MutableBigInteger at the specified index. This method is 374 * not used because it is not inlined on all platforms. 375 */ 376 private final long getLong(int index) { 377 return value[offset+index] & LONG_MASK; 378 } 379 380 /** 381 * Ensure that the MutableBigInteger is in normal form, specifically 382 * making sure that there are no leading zeros, and that if the 383 * magnitude is zero, then intLen is zero. 384 */ 385 final void normalize() { 386 if (intLen == 0) { 387 offset = 0; 388 return; 389 } 390 391 int index = offset; 392 if (value[index] != 0) 393 return; 394 395 int indexBound = index+intLen; 396 do { 397 index++; 398 } while(index < indexBound && value[index]==0); 399 400 int numZeros = index - offset; 401 intLen -= numZeros; 402 offset = (intLen==0 ? 0 : offset+numZeros); 403 } 404 405 /** 406 * If this MutableBigInteger cannot hold len words, increase the size 407 * of the value array to len words. 408 */ 409 private final void ensureCapacity(int len) { 410 if (value.length < len) { 411 value = new int[len]; 412 offset = 0; 413 intLen = len; 414 } 415 } 416 417 /** 418 * Convert this MutableBigInteger into an int array with no leading 419 * zeros, of a length that is equal to this MutableBigInteger's intLen. 420 */ 421 int[] toIntArray() { 422 int[] result = new int[intLen]; 423 for(int i=0; i<intLen; i++) 424 result[i] = value[offset+i]; 425 return result; 426 } 427 428 /** 429 * Sets the int at index+offset in this MutableBigInteger to val. 430 * This does not get inlined on all platforms so it is not used 431 * as often as originally intended. 432 */ 433 void setInt(int index, int val) { 434 value[offset + index] = val; 435 } 436 437 /** 438 * Sets this MutableBigInteger's value array to the specified array. 439 * The intLen is set to the specified length. 440 */ 441 void setValue(int[] val, int length) { 442 value = val; 443 intLen = length; 444 offset = 0; 445 } 446 447 /** 448 * Sets this MutableBigInteger's value array to a copy of the specified 449 * array. The intLen is set to the length of the new array. 450 */ 451 void copyValue(MutableBigInteger src) { 452 int len = src.intLen; 453 if (value.length < len) 454 value = new int[len]; 455 System.arraycopy(src.value, src.offset, value, 0, len); 456 intLen = len; 457 offset = 0; 458 } 459 460 /** 461 * Sets this MutableBigInteger's value array to a copy of the specified 462 * array. The intLen is set to the length of the specified array. 463 */ 464 void copyValue(int[] val) { 465 int len = val.length; 466 if (value.length < len) 467 value = new int[len]; 468 System.arraycopy(val, 0, value, 0, len); 469 intLen = len; 470 offset = 0; 471 } 472 473 /** 474 * Returns true iff this MutableBigInteger has a value of one. 475 */ 476 boolean isOne() { 477 return (intLen == 1) && (value[offset] == 1); 478 } 479 480 /** 481 * Returns true iff this MutableBigInteger has a value of zero. 482 */ 483 boolean isZero() { 484 return (intLen == 0); 485 } 486 487 /** 488 * Returns true iff this MutableBigInteger is even. 489 */ 490 boolean isEven() { 491 return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0); 492 } 493 494 /** 495 * Returns true iff this MutableBigInteger is odd. 496 */ 497 boolean isOdd() { 498 return isZero() ? false : ((value[offset + intLen - 1] & 1) == 1); 499 } 500 501 /** 502 * Returns true iff this MutableBigInteger is in normal form. A 503 * MutableBigInteger is in normal form if it has no leading zeros 504 * after the offset, and intLen + offset <= value.length. 505 */ 506 boolean isNormal() { 507 if (intLen + offset > value.length) 508 return false; 509 if (intLen ==0) 510 return true; 511 return (value[offset] != 0); 512 } 513 514 /** 515 * Returns a String representation of this MutableBigInteger in radix 10. 516 */ 517 public String toString() { 518 BigInteger b = toBigInteger(1); 519 return b.toString(); 520 } 521 522 /** 523 * Like {@link #rightShift(int)} but {@code n} can be greater than the length of the number. 524 */ 525 void safeRightShift(int n) { 526 if (n/32 >= intLen) 527 reset(); 528 else 529 rightShift(n); 530 } 531 532 /** 533 * Right shift this MutableBigInteger n bits. The MutableBigInteger is left 534 * in normal form. 535 */ 536 void rightShift(int n) { 537 if (intLen == 0) 538 return; 539 int nInts = n >>> 5; 540 int nBits = n & 0x1F; 541 this.intLen -= nInts; 542 if (nBits == 0) 543 return; 544 int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]); 545 if (nBits >= bitsInHighWord) { 546 this.primitiveLeftShift(32 - nBits); 547 this.intLen--; 548 } else { 549 primitiveRightShift(nBits); 550 } 551 } 552 553 /** 554 * Like {@link #leftShift(int)} but {@code n} can be zero. 555 */ 556 void safeLeftShift(int n) { 557 if (n > 0) 558 leftShift(n); 559 } 560 561 /** 562 * Left shift this MutableBigInteger n bits. 563 */ 564 void leftShift(int n) { 565 /* 566 * If there is enough storage space in this MutableBigInteger already 567 * the available space will be used. Space to the right of the used 568 * ints in the value array is faster to utilize, so the extra space 569 * will be taken from the right if possible. 570 */ 571 if (intLen == 0) 572 return; 573 int nInts = n >>> 5; 574 int nBits = n&0x1F; 575 int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]); 576 577 // If shift can be done without moving words, do so 578 if (n <= (32-bitsInHighWord)) { 579 primitiveLeftShift(nBits); 580 return; 581 } 582 583 int newLen = intLen + nInts +1; 584 if (nBits <= (32-bitsInHighWord)) 585 newLen--; 586 if (value.length < newLen) { 587 // The array must grow 588 int[] result = new int[newLen]; 589 for (int i=0; i<intLen; i++) 590 result[i] = value[offset+i]; 591 setValue(result, newLen); 592 } else if (value.length - offset >= newLen) { 593 // Use space on right 594 for(int i=0; i<newLen - intLen; i++) 595 value[offset+intLen+i] = 0; 596 } else { 597 // Must use space on left 598 for (int i=0; i<intLen; i++) 599 value[i] = value[offset+i]; 600 for (int i=intLen; i<newLen; i++) 601 value[i] = 0; 602 offset = 0; 603 } 604 intLen = newLen; 605 if (nBits == 0) 606 return; 607 if (nBits <= (32-bitsInHighWord)) 608 primitiveLeftShift(nBits); 609 else 610 primitiveRightShift(32 -nBits); 611 } 612 613 /** 614 * A primitive used for division. This method adds in one multiple of the 615 * divisor a back to the dividend result at a specified offset. It is used 616 * when qhat was estimated too large, and must be adjusted. 617 */ 618 private int divadd(int[] a, int[] result, int offset) { 619 long carry = 0; 620 621 for (int j=a.length-1; j >= 0; j--) { 622 long sum = (a[j] & LONG_MASK) + 623 (result[j+offset] & LONG_MASK) + carry; 624 result[j+offset] = (int)sum; 625 carry = sum >>> 32; 626 } 627 return (int)carry; 628 } 629 630 /** 631 * This method is used for division. It multiplies an n word input a by one 632 * word input x, and subtracts the n word product from q. This is needed 633 * when subtracting qhat*divisor from dividend. 634 */ 635 private int mulsub(int[] q, int[] a, int x, int len, int offset) { 636 long xLong = x & LONG_MASK; 637 long carry = 0; 638 offset += len; 639 640 for (int j=len-1; j >= 0; j--) { 641 long product = (a[j] & LONG_MASK) * xLong + carry; 642 long difference = q[offset] - product; 643 q[offset--] = (int)difference; 644 carry = (product >>> 32) 645 + (((difference & LONG_MASK) > 646 (((~(int)product) & LONG_MASK))) ? 1:0); 647 } 648 return (int)carry; 649 } 650 651 /** 652 * The method is the same as mulsun, except the fact that q array is not 653 * updated, the only result of the method is borrow flag. 654 */ 655 private int mulsubBorrow(int[] q, int[] a, int x, int len, int offset) { 656 long xLong = x & LONG_MASK; 657 long carry = 0; 658 offset += len; 659 for (int j=len-1; j >= 0; j--) { 660 long product = (a[j] & LONG_MASK) * xLong + carry; 661 long difference = q[offset--] - product; 662 carry = (product >>> 32) 663 + (((difference & LONG_MASK) > 664 (((~(int)product) & LONG_MASK))) ? 1:0); 665 } 666 return (int)carry; 667 } 668 669 /** 670 * Right shift this MutableBigInteger n bits, where n is 671 * less than 32. 672 * Assumes that intLen > 0, n > 0 for speed 673 */ 674 private final void primitiveRightShift(int n) { 675 int[] val = value; 676 int n2 = 32 - n; 677 for (int i=offset+intLen-1, c=val[i]; i>offset; i--) { 678 int b = c; 679 c = val[i-1]; 680 val[i] = (c << n2) | (b >>> n); 681 } 682 val[offset] >>>= n; 683 } 684 685 /** 686 * Left shift this MutableBigInteger n bits, where n is 687 * less than 32. 688 * Assumes that intLen > 0, n > 0 for speed 689 */ 690 private final void primitiveLeftShift(int n) { 691 int[] val = value; 692 int n2 = 32 - n; 693 for (int i=offset, c=val[i], m=i+intLen-1; i<m; i++) { 694 int b = c; 695 c = val[i+1]; 696 val[i] = (b << n) | (c >>> n2); 697 } 698 val[offset+intLen-1] <<= n; 699 } 700 701 /** 702 * Returns a {@code BigInteger} equal to the {@code n} 703 * low ints of this number. 704 */ 705 private BigInteger getLower(int n) { 706 if (isZero()) 707 return BigInteger.ZERO; 708 else if (intLen < n) 709 return toBigInteger(1); 710 else { 711 // strip zeros 712 int len = n; 713 while (len>0 && value[offset+intLen-len]==0) 714 len--; 715 int sign = len>0 ? 1 : 0; 716 return new BigInteger(Arrays.copyOfRange(value, offset+intLen-len, offset+intLen), sign); 717 } 718 } 719 720 /** 721 * Discards all ints whose index is greater than {@code n}. 722 */ 723 private void keepLower(int n) { 724 if (intLen >= n) { 725 offset += intLen - n; 726 intLen = n; 727 } 728 } 729 730 /** 731 * Adds the contents of two MutableBigInteger objects.The result 732 * is placed within this MutableBigInteger. 733 * The contents of the addend are not changed. 734 */ 735 void add(MutableBigInteger addend) { 736 int x = intLen; 737 int y = addend.intLen; 738 int resultLen = (intLen > addend.intLen ? intLen : addend.intLen); 739 int[] result = (value.length < resultLen ? new int[resultLen] : value); 740 741 int rstart = result.length-1; 742 long sum; 743 long carry = 0; 744 745 // Add common parts of both numbers 746 while(x>0 && y>0) { 747 x--; y--; 748 sum = (value[x+offset] & LONG_MASK) + 749 (addend.value[y+addend.offset] & LONG_MASK) + carry; 750 result[rstart--] = (int)sum; 751 carry = sum >>> 32; 752 } 753 754 // Add remainder of the longer number 755 while(x>0) { 756 x--; 757 if (carry == 0 && result == value && rstart == (x + offset)) 758 return; 759 sum = (value[x+offset] & LONG_MASK) + carry; 760 result[rstart--] = (int)sum; 761 carry = sum >>> 32; 762 } 763 while(y>0) { 764 y--; 765 sum = (addend.value[y+addend.offset] & LONG_MASK) + carry; 766 result[rstart--] = (int)sum; 767 carry = sum >>> 32; 768 } 769 770 if (carry > 0) { // Result must grow in length 771 resultLen++; 772 if (result.length < resultLen) { 773 int temp[] = new int[resultLen]; 774 // Result one word longer from carry-out; copy low-order 775 // bits into new result. 776 System.arraycopy(result, 0, temp, 1, result.length); 777 temp[0] = 1; 778 result = temp; 779 } else { 780 result[rstart--] = 1; 781 } 782 } 783 784 value = result; 785 intLen = resultLen; 786 offset = result.length - resultLen; 787 } 788 789 /** 790 * Adds the value of {@code addend} shifted {@code n} ints to the left. 791 * Has the same effect as {@code addend.leftShift(32*ints); add(b);} 792 * but doesn't change the value of {@code b}. 793 */ 794 void addShifted(MutableBigInteger addend, int n) { 795 if (addend.isZero()) 796 return; 797 798 int x = intLen; 799 int y = addend.intLen + n; 800 int resultLen = (intLen > y ? intLen : y); 801 int[] result = (value.length < resultLen ? new int[resultLen] : value); 802 803 int rstart = result.length-1; 804 long sum; 805 long carry = 0; 806 807 // Add common parts of both numbers 808 while(x>0 && y>0) { 809 x--; y--; 810 int bval = y+addend.offset<addend.value.length ? addend.value[y+addend.offset] : 0; 811 sum = (value[x+offset] & LONG_MASK) + 812 (bval & LONG_MASK) + carry; 813 result[rstart--] = (int)sum; 814 carry = sum >>> 32; 815 } 816 817 // Add remainder of the longer number 818 while(x>0) { 819 x--; 820 if (carry == 0 && result == value && rstart == (x + offset)) 821 return; 822 sum = (value[x+offset] & LONG_MASK) + carry; 823 result[rstart--] = (int)sum; 824 carry = sum >>> 32; 825 } 826 while(y>0) { 827 y--; 828 int bval = y+addend.offset<addend.value.length ? addend.value[y+addend.offset] : 0; 829 sum = (bval & LONG_MASK) + carry; 830 result[rstart--] = (int)sum; 831 carry = sum >>> 32; 832 } 833 834 if (carry > 0) { // Result must grow in length 835 resultLen++; 836 if (result.length < resultLen) { 837 int temp[] = new int[resultLen]; 838 // Result one word longer from carry-out; copy low-order 839 // bits into new result. 840 System.arraycopy(result, 0, temp, 1, result.length); 841 temp[0] = 1; 842 result = temp; 843 } else { 844 result[rstart--] = 1; 845 } 846 } 847 848 value = result; 849 intLen = resultLen; 850 offset = result.length - resultLen; 851 } 852 853 /** 854 * Like {@link #addShifted(MutableBigInteger, int)} but {@code this.intLen} must 855 * not be greater than {@code n}. In other words, concatenates {@code this} 856 * and {@code addend}. 857 */ 858 void addDisjoint(MutableBigInteger addend, int n) { 859 if (addend.isZero()) 860 return; 861 862 int x = intLen; 863 int y = addend.intLen + n; 864 int resultLen = (intLen > y ? intLen : y); 865 int[] result; 866 if (value.length < resultLen) 867 result = new int[resultLen]; 868 else { 869 result = value; 870 Arrays.fill(value, offset+intLen, value.length, 0); 871 } 872 873 int rstart = result.length-1; 874 875 // copy from this if needed 876 System.arraycopy(value, offset, result, rstart+1-x, x); 877 y -= x; 878 rstart -= x; 879 880 int len = Math.min(y, addend.value.length-addend.offset); 881 System.arraycopy(addend.value, addend.offset, result, rstart+1-y, len); 882 883 // zero the gap 884 for (int i=rstart+1-y+len; i<rstart+1; i++) 885 result[i] = 0; 886 887 value = result; 888 intLen = resultLen; 889 offset = result.length - resultLen; 890 } 891 892 /** 893 * Adds the low {@code n} ints of {@code addend}. 894 */ 895 void addLower(MutableBigInteger addend, int n) { 896 MutableBigInteger a = new MutableBigInteger(addend); 897 if (a.offset + a.intLen >= n) { 898 a.offset = a.offset + a.intLen - n; 899 a.intLen = n; 900 } 901 a.normalize(); 902 add(a); 903 } 904 905 /** 906 * Subtracts the smaller of this and b from the larger and places the 907 * result into this MutableBigInteger. 908 */ 909 int subtract(MutableBigInteger b) { 910 MutableBigInteger a = this; 911 912 int[] result = value; 913 int sign = a.compare(b); 914 915 if (sign == 0) { 916 reset(); 917 return 0; 918 } 919 if (sign < 0) { 920 MutableBigInteger tmp = a; 921 a = b; 922 b = tmp; 923 } 924 925 int resultLen = a.intLen; 926 if (result.length < resultLen) 927 result = new int[resultLen]; 928 929 long diff = 0; 930 int x = a.intLen; 931 int y = b.intLen; 932 int rstart = result.length - 1; 933 934 // Subtract common parts of both numbers 935 while (y>0) { 936 x--; y--; 937 938 diff = (a.value[x+a.offset] & LONG_MASK) - 939 (b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32)); 940 result[rstart--] = (int)diff; 941 } 942 // Subtract remainder of longer number 943 while (x>0) { 944 x--; 945 diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32)); 946 result[rstart--] = (int)diff; 947 } 948 949 value = result; 950 intLen = resultLen; 951 offset = value.length - resultLen; 952 normalize(); 953 return sign; 954 } 955 956 /** 957 * Subtracts the smaller of a and b from the larger and places the result 958 * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no 959 * operation was performed. 960 */ 961 private int difference(MutableBigInteger b) { 962 MutableBigInteger a = this; 963 int sign = a.compare(b); 964 if (sign ==0) 965 return 0; 966 if (sign < 0) { 967 MutableBigInteger tmp = a; 968 a = b; 969 b = tmp; 970 } 971 972 long diff = 0; 973 int x = a.intLen; 974 int y = b.intLen; 975 976 // Subtract common parts of both numbers 977 while (y>0) { 978 x--; y--; 979 diff = (a.value[a.offset+ x] & LONG_MASK) - 980 (b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32)); 981 a.value[a.offset+x] = (int)diff; 982 } 983 // Subtract remainder of longer number 984 while (x>0) { 985 x--; 986 diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32)); 987 a.value[a.offset+x] = (int)diff; 988 } 989 990 a.normalize(); 991 return sign; 992 } 993 994 /** 995 * Multiply the contents of two MutableBigInteger objects. The result is 996 * placed into MutableBigInteger z. The contents of y are not changed. 997 */ 998 void multiply(MutableBigInteger y, MutableBigInteger z) { 999 int xLen = intLen; 1000 int yLen = y.intLen; 1001 int newLen = xLen + yLen; 1002 1003 // Put z into an appropriate state to receive product 1004 if (z.value.length < newLen) 1005 z.value = new int[newLen]; 1006 z.offset = 0; 1007 z.intLen = newLen; 1008 1009 // The first iteration is hoisted out of the loop to avoid extra add 1010 long carry = 0; 1011 for (int j=yLen-1, k=yLen+xLen-1; j >= 0; j--, k--) { 1012 long product = (y.value[j+y.offset] & LONG_MASK) * 1013 (value[xLen-1+offset] & LONG_MASK) + carry; 1014 z.value[k] = (int)product; 1015 carry = product >>> 32; 1016 } 1017 z.value[xLen-1] = (int)carry; 1018 1019 // Perform the multiplication word by word 1020 for (int i = xLen-2; i >= 0; i--) { 1021 carry = 0; 1022 for (int j=yLen-1, k=yLen+i; j >= 0; j--, k--) { 1023 long product = (y.value[j+y.offset] & LONG_MASK) * 1024 (value[i+offset] & LONG_MASK) + 1025 (z.value[k] & LONG_MASK) + carry; 1026 z.value[k] = (int)product; 1027 carry = product >>> 32; 1028 } 1029 z.value[i] = (int)carry; 1030 } 1031 1032 // Remove leading zeros from product 1033 z.normalize(); 1034 } 1035 1036 /** 1037 * Multiply the contents of this MutableBigInteger by the word y. The 1038 * result is placed into z. 1039 */ 1040 void mul(int y, MutableBigInteger z) { 1041 if (y == 1) { 1042 z.copyValue(this); 1043 return; 1044 } 1045 1046 if (y == 0) { 1047 z.clear(); 1048 return; 1049 } 1050 1051 // Perform the multiplication word by word 1052 long ylong = y & LONG_MASK; 1053 int[] zval = (z.value.length<intLen+1 ? new int[intLen + 1] 1054 : z.value); 1055 long carry = 0; 1056 for (int i = intLen-1; i >= 0; i--) { 1057 long product = ylong * (value[i+offset] & LONG_MASK) + carry; 1058 zval[i+1] = (int)product; 1059 carry = product >>> 32; 1060 } 1061 1062 if (carry == 0) { 1063 z.offset = 1; 1064 z.intLen = intLen; 1065 } else { 1066 z.offset = 0; 1067 z.intLen = intLen + 1; 1068 zval[0] = (int)carry; 1069 } 1070 z.value = zval; 1071 } 1072 1073 /** 1074 * This method is used for division of an n word dividend by a one word 1075 * divisor. The quotient is placed into quotient. The one word divisor is 1076 * specified by divisor. 1077 * 1078 * @return the remainder of the division is returned. 1079 * 1080 */ 1081 int divideOneWord(int divisor, MutableBigInteger quotient) { 1082 long divisorLong = divisor & LONG_MASK; 1083 1084 // Special case of one word dividend 1085 if (intLen == 1) { 1086 long dividendValue = value[offset] & LONG_MASK; 1087 int q = (int) (dividendValue / divisorLong); 1088 int r = (int) (dividendValue - q * divisorLong); 1089 quotient.value[0] = q; 1090 quotient.intLen = (q == 0) ? 0 : 1; 1091 quotient.offset = 0; 1092 return r; 1093 } 1094 1095 if (quotient.value.length < intLen) 1096 quotient.value = new int[intLen]; 1097 quotient.offset = 0; 1098 quotient.intLen = intLen; 1099 1100 // Normalize the divisor 1101 int shift = Integer.numberOfLeadingZeros(divisor); 1102 1103 int rem = value[offset]; 1104 long remLong = rem & LONG_MASK; 1105 if (remLong < divisorLong) { 1106 quotient.value[0] = 0; 1107 } else { 1108 quotient.value[0] = (int)(remLong / divisorLong); 1109 rem = (int) (remLong - (quotient.value[0] * divisorLong)); 1110 remLong = rem & LONG_MASK; 1111 } 1112 int xlen = intLen; 1113 while (--xlen > 0) { 1114 long dividendEstimate = (remLong << 32) | 1115 (value[offset + intLen - xlen] & LONG_MASK); 1116 int q; 1117 if (dividendEstimate >= 0) { 1118 q = (int) (dividendEstimate / divisorLong); 1119 rem = (int) (dividendEstimate - q * divisorLong); 1120 } else { 1121 long tmp = divWord(dividendEstimate, divisor); 1122 q = (int) (tmp & LONG_MASK); 1123 rem = (int) (tmp >>> 32); 1124 } 1125 quotient.value[intLen - xlen] = q; 1126 remLong = rem & LONG_MASK; 1127 } 1128 1129 quotient.normalize(); 1130 // Unnormalize 1131 if (shift > 0) 1132 return rem % divisor; 1133 else 1134 return rem; 1135 } 1136 1137 /** 1138 * Calculates the quotient of this div b and places the quotient in the 1139 * provided MutableBigInteger objects and the remainder object is returned. 1140 * 1141 */ 1142 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) { 1143 return divide(b,quotient,true); 1144 } 1145 1146 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) { 1147 if (intLen<BigInteger.BURNIKEL_ZIEGLER_THRESHOLD || b.intLen<BigInteger.BURNIKEL_ZIEGLER_THRESHOLD) 1148 return divideKnuth(b, quotient, needRemainder); 1149 else 1150 return divideAndRemainderBurnikelZiegler(b, quotient); 1151 } 1152 1153 /** 1154 * @see #divideKnuth(MutableBigInteger, MutableBigInteger, boolean) 1155 */ 1156 MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient) { 1157 return divideKnuth(b,quotient,true); 1158 } 1159 1160 /** 1161 * Calculates the quotient of this div b and places the quotient in the 1162 * provided MutableBigInteger objects and the remainder object is returned. 1163 * 1164 * Uses Algorithm D in Knuth section 4.3.1. 1165 * Many optimizations to that algorithm have been adapted from the Colin 1166 * Plumb C library. 1167 * It special cases one word divisors for speed. The content of b is not 1168 * changed. 1169 * 1170 */ 1171 MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) { 1172 if (b.intLen == 0) 1173 throw new ArithmeticException("BigInteger divide by zero"); 1174 1175 // Dividend is zero 1176 if (intLen == 0) { 1177 quotient.intLen = quotient.offset = 0; 1178 return needRemainder ? new MutableBigInteger() : null; 1179 } 1180 1181 int cmp = compare(b); 1182 // Dividend less than divisor 1183 if (cmp < 0) { 1184 quotient.intLen = quotient.offset = 0; 1185 return needRemainder ? new MutableBigInteger(this) : null; 1186 } 1187 // Dividend equal to divisor 1188 if (cmp == 0) { 1189 quotient.value[0] = quotient.intLen = 1; 1190 quotient.offset = 0; 1191 return needRemainder ? new MutableBigInteger() : null; 1192 } 1193 1194 quotient.clear(); 1195 // Special case one word divisor 1196 if (b.intLen == 1) { 1197 int r = divideOneWord(b.value[b.offset], quotient); 1198 if(needRemainder) { 1199 if (r == 0) 1200 return new MutableBigInteger(); 1201 return new MutableBigInteger(r); 1202 } else { 1203 return null; 1204 } 1205 } 1206 1207 // Cancel common powers of two if we're above the KNUTH_POW2_* thresholds 1208 if (intLen >= KNUTH_POW2_THRESH_LEN) { 1209 int trailingZeroBits = Math.min(getLowestSetBit(), b.getLowestSetBit()); 1210 if (trailingZeroBits >= KNUTH_POW2_THRESH_ZEROS*32) { 1211 MutableBigInteger a = new MutableBigInteger(this); 1212 b = new MutableBigInteger(b); 1213 a.rightShift(trailingZeroBits); 1214 b.rightShift(trailingZeroBits); 1215 MutableBigInteger r = a.divideKnuth(b, quotient); 1216 r.leftShift(trailingZeroBits); 1217 return r; 1218 } 1219 } 1220 1221 return divideMagnitude(b, quotient, needRemainder); 1222 } 1223 1224 /** 1225 * Computes {@code this/b} and {@code this%b} using the 1226 * <a href="http://cr.yp.to/bib/1998/burnikel.ps"> Burnikel-Ziegler algorithm</a>. 1227 * This method implements algorithm 3 from pg. 9 of the Burnikel-Ziegler paper. 1228 * The parameter beta was chosen to b 2<sup>32</sup> so almost all shifts are 1229 * multiples of 32 bits.<br/> 1230 * {@code this} and {@code b} must be nonnegative. 1231 * @param b the divisor 1232 * @param quotient output parameter for {@code this/b} 1233 * @return the remainder 1234 */ 1235 MutableBigInteger divideAndRemainderBurnikelZiegler(MutableBigInteger b, MutableBigInteger quotient) { 1236 int r = intLen; 1237 int s = b.intLen; 1238 1239 if (r < s) 1240 return this; 1241 else { 1242 // Unlike Knuth division, we don't check for common powers of two here because 1243 // BZ already runs faster if both numbers contain powers of two and cancelling them has no 1244 // additional benefit. 1245 1246 // step 1: let m = min{2^k | (2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s} 1247 int m = 1 << (32-Integer.numberOfLeadingZeros(s/BigInteger.BURNIKEL_ZIEGLER_THRESHOLD)); 1248 1249 int j = (s+m-1) / m; // step 2a: j = ceil(s/m) 1250 int n = j * m; // step 2b: block length in 32-bit units 1251 int n32 = 32 * n; // block length in bits 1252 int sigma = Math.max(0, n32 - b.bitLength()); // step 3: sigma = max{T | (2^T)*B < beta^n} 1253 MutableBigInteger bShifted = new MutableBigInteger(b); 1254 bShifted.safeLeftShift(sigma); // step 4a: shift b so its length is a multiple of n 1255 safeLeftShift(sigma); // step 4b: shift this by the same amount 1256 1257 // step 5: t is the number of blocks needed to accommodate this plus one additional bit 1258 int t = (bitLength()+n32) / n32; 1259 if (t < 2) 1260 t = 2; 1261 1262 // step 6: conceptually split this into blocks a[t-1], ..., a[0] 1263 MutableBigInteger a1 = getBlock(t-1, t, n); // the most significant block of this 1264 1265 // step 7: z[t-2] = [a[t-1], a[t-2]] 1266 MutableBigInteger z = getBlock(t-2, t, n); // the second to most significant block 1267 z.addDisjoint(a1, n); // z[t-2] 1268 1269 // do schoolbook division on blocks, dividing 2-block numbers by 1-block numbers 1270 MutableBigInteger qi = new MutableBigInteger(); 1271 MutableBigInteger ri; 1272 quotient.offset = quotient.intLen = 0; 1273 for (int i=t-2; i>0; i--) { 1274 // step 8a: compute (qi,ri) such that z=b*qi+ri 1275 ri = z.divide2n1n(bShifted, qi); 1276 1277 // step 8b: z = [ri, a[i-1]] 1278 z = getBlock(i-1, t, n); // a[i-1] 1279 z.addDisjoint(ri, n); 1280 quotient.addShifted(qi, i*n); // update q (part of step 9) 1281 } 1282 // final iteration of step 8: do the loop one more time for i=0 but leave z unchanged 1283 ri = z.divide2n1n(bShifted, qi); 1284 quotient.add(qi); 1285 1286 ri.rightShift(sigma); // step 9: this and b were shifted, so shift back 1287 return ri; 1288 } 1289 } 1290 1291 /** 1292 * This method implements algorithm 1 from pg. 4 of the Burnikel-Ziegler paper. 1293 * It divides a 2n-digit number by a n-digit number.<br/> 1294 * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits. 1295 * <br/> 1296 * {@code this} must be a nonnegative number such that {@code this.bitLength() <= 2*b.bitLength()} 1297 * @param b a positive number such that {@code b.bitLength()} is even 1298 * @param quotient output parameter for {@code this/b} 1299 * @return {@code this%b} 1300 */ 1301 private MutableBigInteger divide2n1n(MutableBigInteger b, MutableBigInteger quotient) { 1302 int n = b.intLen; 1303 1304 // step 1: base case 1305 if (n%2!=0 || n<BigInteger.BURNIKEL_ZIEGLER_THRESHOLD) 1306 return divideKnuth(b, quotient); 1307 1308 // step 2: view this as [a1,a2,a3,a4] where each ai is n/2 ints or less 1309 MutableBigInteger aUpper = new MutableBigInteger(this); 1310 aUpper.safeRightShift(32*(n/2)); // aUpper = [a1,a2,a3] 1311 keepLower(n/2); // this = a4 1312 1313 // step 3: q1=aUpper/b, r1=aUpper%b 1314 MutableBigInteger q1 = new MutableBigInteger(); 1315 MutableBigInteger r1 = aUpper.divide3n2n(b, q1); 1316 1317 // step 4: quotient=[r1,this]/b, r2=[r1,this]%b 1318 addDisjoint(r1, n/2); // this = [r1,this] 1319 MutableBigInteger r2 = divide3n2n(b, quotient); 1320 1321 // step 5: let quotient=[q1,quotient] and return r2 1322 quotient.addDisjoint(q1, n/2); 1323 return r2; 1324 } 1325 1326 /** 1327 * This method implements algorithm 2 from pg. 5 of the Burnikel-Ziegler paper. 1328 * It divides a 3n-digit number by a 2n-digit number.<br/> 1329 * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.<br/> 1330 * <br/> 1331 * {@code this} must be a nonnegative number such that {@code 2*this.bitLength() <= 3*b.bitLength()} 1332 * @param quotient output parameter for {@code this/b} 1333 * @return {@code this%b} 1334 */ 1335 private MutableBigInteger divide3n2n(MutableBigInteger b, MutableBigInteger quotient) { 1336 int n = b.intLen / 2; // half the length of b in ints 1337 1338 // step 1: view this as [a1,a2,a3] where each ai is n ints or less; let a12=[a1,a2] 1339 MutableBigInteger a12 = new MutableBigInteger(this); 1340 a12.safeRightShift(32*n); 1341 1342 // step 2: view b as [b1,b2] where each bi is n ints or less 1343 MutableBigInteger b1 = new MutableBigInteger(b); 1344 b1.safeRightShift(n * 32); 1345 BigInteger b2 = b.getLower(n); 1346 1347 MutableBigInteger r; 1348 MutableBigInteger d; 1349 if (compareShifted(b, n) < 0) { 1350 // step 3a: if a1<b1, let quotient=a12/b1 and r=a12%b1 1351 r = a12.divide2n1n(b1, quotient); 1352 1353 // step 4: d=quotient*b2 1354 d = new MutableBigInteger(quotient.toBigInteger().multiply(b2)); 1355 } 1356 else { 1357 // step 3b: if a1>=b1, let quotient=beta^n-1 and r=a12-b1*2^n+b1 1358 quotient.ones(n); 1359 a12.add(b1); 1360 b1.leftShift(32*n); 1361 a12.subtract(b1); 1362 r = a12; 1363 1364 // step 4: d=quotient*b2=(b2 << 32*n) - b2 1365 d = new MutableBigInteger(b2); 1366 d.leftShift(32 * n); 1367 d.subtract(new MutableBigInteger(b2)); 1368 } 1369 1370 // step 5: r = r*beta^n + a3 - d (paper says a4) 1371 // However, don't subtract d until after the while loop so r doesn't become negative 1372 r.leftShift(32 * n); 1373 r.addLower(this, n); 1374 1375 // step 6: add b until r>=d 1376 while (r.compare(d) < 0) { 1377 r.add(b); 1378 quotient.subtract(MutableBigInteger.ONE); 1379 } 1380 r.subtract(d); 1381 1382 return r; 1383 } 1384 1385 /** 1386 * Returns a {@code MutableBigInteger} containing {@code blockLength} ints from 1387 * {@code this} number, starting at {@code index*blockLength}.<br/> 1388 * Used by Burnikel-Ziegler division. 1389 * @param index the block index 1390 * @param numBlocks the total number of blocks in {@code this} number 1391 * @param blockLength length of one block in units of 32 bits 1392 * @return 1393 */ 1394 private MutableBigInteger getBlock(int index, int numBlocks, int blockLength) { 1395 int blockStart = index * blockLength; 1396 if (blockStart >= intLen) 1397 return new MutableBigInteger(); 1398 1399 int blockEnd; 1400 if (index == numBlocks-1) 1401 blockEnd = intLen; 1402 else 1403 blockEnd = (index+1) * blockLength; 1404 if (blockEnd > intLen) 1405 return new MutableBigInteger(); 1406 1407 int[] newVal = Arrays.copyOfRange(value, offset+intLen-blockEnd, offset+intLen-blockStart); 1408 return new MutableBigInteger(newVal); 1409 } 1410 1411 /** @see BigInteger#bitLength() */ 1412 int bitLength() { 1413 if (intLen == 0) 1414 return 0; 1415 return intLen*32 - Integer.numberOfLeadingZeros(value[offset]); 1416 } 1417 1418 /** 1419 * Internally used to calculate the quotient of this div v and places the 1420 * quotient in the provided MutableBigInteger object and the remainder is 1421 * returned. 1422 * 1423 * @return the remainder of the division will be returned. 1424 */ 1425 long divide(long v, MutableBigInteger quotient) { 1426 if (v == 0) 1427 throw new ArithmeticException("BigInteger divide by zero"); 1428 1429 // Dividend is zero 1430 if (intLen == 0) { 1431 quotient.intLen = quotient.offset = 0; 1432 return 0; 1433 } 1434 if (v < 0) 1435 v = -v; 1436 1437 int d = (int)(v >>> 32); 1438 quotient.clear(); 1439 // Special case on word divisor 1440 if (d == 0) 1441 return divideOneWord((int)v, quotient) & LONG_MASK; 1442 else { 1443 return divideLongMagnitude(v, quotient).toLong(); 1444 } 1445 } 1446 1447 private static void copyAndShift(int[] src, int srcFrom, int srcLen, int[] dst, int dstFrom, int shift) { 1448 int n2 = 32 - shift; 1449 int c=src[srcFrom]; 1450 for (int i=0; i < srcLen-1; i++) { 1451 int b = c; 1452 c = src[++srcFrom]; 1453 dst[dstFrom+i] = (b << shift) | (c >>> n2); 1454 } 1455 dst[dstFrom+srcLen-1] = c << shift; 1456 } 1457 1458 /** 1459 * Divide this MutableBigInteger by the divisor. 1460 * The quotient will be placed into the provided quotient object & 1461 * the remainder object is returned. 1462 */ 1463 private MutableBigInteger divideMagnitude(MutableBigInteger div, 1464 MutableBigInteger quotient, 1465 boolean needRemainder ) { 1466 // assert div.intLen > 1 1467 // D1 normalize the divisor 1468 int shift = Integer.numberOfLeadingZeros(div.value[div.offset]); 1469 // Copy divisor value to protect divisor 1470 final int dlen = div.intLen; 1471 int[] divisor; 1472 MutableBigInteger rem; // Remainder starts as dividend with space for a leading zero 1473 if (shift > 0) { 1474 divisor = new int[dlen]; 1475 copyAndShift(div.value,div.offset,dlen,divisor,0,shift); 1476 if(Integer.numberOfLeadingZeros(value[offset])>=shift) { 1477 int[] remarr = new int[intLen + 1]; 1478 rem = new MutableBigInteger(remarr); 1479 rem.intLen = intLen; 1480 rem.offset = 1; 1481 copyAndShift(value,offset,intLen,remarr,1,shift); 1482 } else { 1483 int[] remarr = new int[intLen + 2]; 1484 rem = new MutableBigInteger(remarr); 1485 rem.intLen = intLen+1; 1486 rem.offset = 1; 1487 int rFrom = offset; 1488 int c=0; 1489 int n2 = 32 - shift; 1490 for (int i=1; i < intLen+1; i++,rFrom++) { 1491 int b = c; 1492 c = value[rFrom]; 1493 remarr[i] = (b << shift) | (c >>> n2); 1494 } 1495 remarr[intLen+1] = c << shift; 1496 } 1497 } else { 1498 divisor = Arrays.copyOfRange(div.value, div.offset, div.offset + div.intLen); 1499 rem = new MutableBigInteger(new int[intLen + 1]); 1500 System.arraycopy(value, offset, rem.value, 1, intLen); 1501 rem.intLen = intLen; 1502 rem.offset = 1; 1503 } 1504 1505 int nlen = rem.intLen; 1506 1507 // Set the quotient size 1508 final int limit = nlen - dlen + 1; 1509 if (quotient.value.length < limit) { 1510 quotient.value = new int[limit]; 1511 quotient.offset = 0; 1512 } 1513 quotient.intLen = limit; 1514 int[] q = quotient.value; 1515 1516 1517 // Must insert leading 0 in rem if its length did not change 1518 if (rem.intLen == nlen) { 1519 rem.offset = 0; 1520 rem.value[0] = 0; 1521 rem.intLen++; 1522 } 1523 1524 int dh = divisor[0]; 1525 long dhLong = dh & LONG_MASK; 1526 int dl = divisor[1]; 1527 1528 // D2 Initialize j 1529 for(int j=0; j<limit-1; j++) { 1530 // D3 Calculate qhat 1531 // estimate qhat 1532 int qhat = 0; 1533 int qrem = 0; 1534 boolean skipCorrection = false; 1535 int nh = rem.value[j+rem.offset]; 1536 int nh2 = nh + 0x80000000; 1537 int nm = rem.value[j+1+rem.offset]; 1538 1539 if (nh == dh) { 1540 qhat = ~0; 1541 qrem = nh + nm; 1542 skipCorrection = qrem + 0x80000000 < nh2; 1543 } else { 1544 long nChunk = (((long)nh) << 32) | (nm & LONG_MASK); 1545 if (nChunk >= 0) { 1546 qhat = (int) (nChunk / dhLong); 1547 qrem = (int) (nChunk - (qhat * dhLong)); 1548 } else { 1549 long tmp = divWord(nChunk, dh); 1550 qhat = (int) (tmp & LONG_MASK); 1551 qrem = (int) (tmp >>> 32); 1552 } 1553 } 1554 1555 if (qhat == 0) 1556 continue; 1557 1558 if (!skipCorrection) { // Correct qhat 1559 long nl = rem.value[j+2+rem.offset] & LONG_MASK; 1560 long rs = ((qrem & LONG_MASK) << 32) | nl; 1561 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK); 1562 1563 if (unsignedLongCompare(estProduct, rs)) { 1564 qhat--; 1565 qrem = (int)((qrem & LONG_MASK) + dhLong); 1566 if ((qrem & LONG_MASK) >= dhLong) { 1567 estProduct -= (dl & LONG_MASK); 1568 rs = ((qrem & LONG_MASK) << 32) | nl; 1569 if (unsignedLongCompare(estProduct, rs)) 1570 qhat--; 1571 } 1572 } 1573 } 1574 1575 // D4 Multiply and subtract 1576 rem.value[j+rem.offset] = 0; 1577 int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset); 1578 1579 // D5 Test remainder 1580 if (borrow + 0x80000000 > nh2) { 1581 // D6 Add back 1582 divadd(divisor, rem.value, j+1+rem.offset); 1583 qhat--; 1584 } 1585 1586 // Store the quotient digit 1587 q[j] = qhat; 1588 } // D7 loop on j 1589 // D3 Calculate qhat 1590 // estimate qhat 1591 int qhat = 0; 1592 int qrem = 0; 1593 boolean skipCorrection = false; 1594 int nh = rem.value[limit - 1 + rem.offset]; 1595 int nh2 = nh + 0x80000000; 1596 int nm = rem.value[limit + rem.offset]; 1597 1598 if (nh == dh) { 1599 qhat = ~0; 1600 qrem = nh + nm; 1601 skipCorrection = qrem + 0x80000000 < nh2; 1602 } else { 1603 long nChunk = (((long) nh) << 32) | (nm & LONG_MASK); 1604 if (nChunk >= 0) { 1605 qhat = (int) (nChunk / dhLong); 1606 qrem = (int) (nChunk - (qhat * dhLong)); 1607 } else { 1608 long tmp = divWord(nChunk, dh); 1609 qhat = (int) (tmp & LONG_MASK); 1610 qrem = (int) (tmp >>> 32); 1611 } 1612 } 1613 if (qhat != 0) { 1614 if (!skipCorrection) { // Correct qhat 1615 long nl = rem.value[limit + 1 + rem.offset] & LONG_MASK; 1616 long rs = ((qrem & LONG_MASK) << 32) | nl; 1617 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK); 1618 1619 if (unsignedLongCompare(estProduct, rs)) { 1620 qhat--; 1621 qrem = (int) ((qrem & LONG_MASK) + dhLong); 1622 if ((qrem & LONG_MASK) >= dhLong) { 1623 estProduct -= (dl & LONG_MASK); 1624 rs = ((qrem & LONG_MASK) << 32) | nl; 1625 if (unsignedLongCompare(estProduct, rs)) 1626 qhat--; 1627 } 1628 } 1629 } 1630 1631 1632 // D4 Multiply and subtract 1633 int borrow; 1634 rem.value[limit - 1 + rem.offset] = 0; 1635 if(needRemainder) 1636 borrow = mulsub(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset); 1637 else 1638 borrow = mulsubBorrow(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset); 1639 1640 // D5 Test remainder 1641 if (borrow + 0x80000000 > nh2) { 1642 // D6 Add back 1643 if(needRemainder) 1644 divadd(divisor, rem.value, limit - 1 + 1 + rem.offset); 1645 qhat--; 1646 } 1647 1648 // Store the quotient digit 1649 q[(limit - 1)] = qhat; 1650 } 1651 1652 1653 if(needRemainder) { 1654 // D8 Unnormalize 1655 if (shift > 0) 1656 rem.rightShift(shift); 1657 rem.normalize(); 1658 } 1659 quotient.normalize(); 1660 return needRemainder ? rem : null; 1661 } 1662 1663 /** 1664 * Divide this MutableBigInteger by the divisor represented by positive long 1665 * value. The quotient will be placed into the provided quotient object & 1666 * the remainder object is returned. 1667 */ 1668 private MutableBigInteger divideLongMagnitude(long ldivisor, MutableBigInteger quotient) { 1669 // Remainder starts as dividend with space for a leading zero 1670 MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]); 1671 System.arraycopy(value, offset, rem.value, 1, intLen); 1672 rem.intLen = intLen; 1673 rem.offset = 1; 1674 1675 int nlen = rem.intLen; 1676 1677 int limit = nlen - 2 + 1; 1678 if (quotient.value.length < limit) { 1679 quotient.value = new int[limit]; 1680 quotient.offset = 0; 1681 } 1682 quotient.intLen = limit; 1683 int[] q = quotient.value; 1684 1685 // D1 normalize the divisor 1686 int shift = Long.numberOfLeadingZeros(ldivisor); 1687 if (shift > 0) { 1688 ldivisor<<=shift; 1689 rem.leftShift(shift); 1690 } 1691 1692 // Must insert leading 0 in rem if its length did not change 1693 if (rem.intLen == nlen) { 1694 rem.offset = 0; 1695 rem.value[0] = 0; 1696 rem.intLen++; 1697 } 1698 1699 int dh = (int)(ldivisor >>> 32); 1700 long dhLong = dh & LONG_MASK; 1701 int dl = (int)(ldivisor & LONG_MASK); 1702 1703 // D2 Initialize j 1704 for (int j = 0; j < limit; j++) { 1705 // D3 Calculate qhat 1706 // estimate qhat 1707 int qhat = 0; 1708 int qrem = 0; 1709 boolean skipCorrection = false; 1710 int nh = rem.value[j + rem.offset]; 1711 int nh2 = nh + 0x80000000; 1712 int nm = rem.value[j + 1 + rem.offset]; 1713 1714 if (nh == dh) { 1715 qhat = ~0; 1716 qrem = nh + nm; 1717 skipCorrection = qrem + 0x80000000 < nh2; 1718 } else { 1719 long nChunk = (((long) nh) << 32) | (nm & LONG_MASK); 1720 if (nChunk >= 0) { 1721 qhat = (int) (nChunk / dhLong); 1722 qrem = (int) (nChunk - (qhat * dhLong)); 1723 } else { 1724 long tmp = divWord(nChunk, dh); 1725 qhat =(int)(tmp & LONG_MASK); 1726 qrem = (int)(tmp>>>32); 1727 } 1728 } 1729 1730 if (qhat == 0) 1731 continue; 1732 1733 if (!skipCorrection) { // Correct qhat 1734 long nl = rem.value[j + 2 + rem.offset] & LONG_MASK; 1735 long rs = ((qrem & LONG_MASK) << 32) | nl; 1736 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK); 1737 1738 if (unsignedLongCompare(estProduct, rs)) { 1739 qhat--; 1740 qrem = (int) ((qrem & LONG_MASK) + dhLong); 1741 if ((qrem & LONG_MASK) >= dhLong) { 1742 estProduct -= (dl & LONG_MASK); 1743 rs = ((qrem & LONG_MASK) << 32) | nl; 1744 if (unsignedLongCompare(estProduct, rs)) 1745 qhat--; 1746 } 1747 } 1748 } 1749 1750 // D4 Multiply and subtract 1751 rem.value[j + rem.offset] = 0; 1752 int borrow = mulsubLong(rem.value, dh, dl, qhat, j + rem.offset); 1753 1754 // D5 Test remainder 1755 if (borrow + 0x80000000 > nh2) { 1756 // D6 Add back 1757 divaddLong(dh,dl, rem.value, j + 1 + rem.offset); 1758 qhat--; 1759 } 1760 1761 // Store the quotient digit 1762 q[j] = qhat; 1763 } // D7 loop on j 1764 1765 // D8 Unnormalize 1766 if (shift > 0) 1767 rem.rightShift(shift); 1768 1769 quotient.normalize(); 1770 rem.normalize(); 1771 return rem; 1772 } 1773 1774 /** 1775 * A primitive used for division by long. 1776 * Specialized version of the method divadd. 1777 * dh is a high part of the divisor, dl is a low part 1778 */ 1779 private int divaddLong(int dh, int dl, int[] result, int offset) { 1780 long carry = 0; 1781 1782 long sum = (dl & LONG_MASK) + (result[1+offset] & LONG_MASK); 1783 result[1+offset] = (int)sum; 1784 1785 sum = (dh & LONG_MASK) + (result[offset] & LONG_MASK) + carry; 1786 result[offset] = (int)sum; 1787 carry = sum >>> 32; 1788 return (int)carry; 1789 } 1790 1791 /** 1792 * This method is used for division by long. 1793 * Specialized version of the method sulsub. 1794 * dh is a high part of the divisor, dl is a low part 1795 */ 1796 private int mulsubLong(int[] q, int dh, int dl, int x, int offset) { 1797 long xLong = x & LONG_MASK; 1798 offset += 2; 1799 long product = (dl & LONG_MASK) * xLong; 1800 long difference = q[offset] - product; 1801 q[offset--] = (int)difference; 1802 long carry = (product >>> 32) 1803 + (((difference & LONG_MASK) > 1804 (((~(int)product) & LONG_MASK))) ? 1:0); 1805 product = (dh & LONG_MASK) * xLong + carry; 1806 difference = q[offset] - product; 1807 q[offset--] = (int)difference; 1808 carry = (product >>> 32) 1809 + (((difference & LONG_MASK) > 1810 (((~(int)product) & LONG_MASK))) ? 1:0); 1811 return (int)carry; 1812 } 1813 1814 /** 1815 * Compare two longs as if they were unsigned. 1816 * Returns true iff one is bigger than two. 1817 */ 1818 private boolean unsignedLongCompare(long one, long two) { 1819 return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE); 1820 } 1821 1822 /** 1823 * This method divides a long quantity by an int to estimate 1824 * qhat for two multi precision numbers. It is used when 1825 * the signed value of n is less than zero. 1826 * Returns long value where high 32 bits contain remainder value and 1827 * low 32 bits contain quotient value. 1828 */ 1829 static long divWord(long n, int d) { 1830 long dLong = d & LONG_MASK; 1831 long r; 1832 long q; 1833 if (dLong == 1) { 1834 q = (int)n; 1835 r = 0; 1836 return (r << 32) | (q & LONG_MASK); 1837 } 1838 1839 // Approximate the quotient and remainder 1840 q = (n >>> 1) / (dLong >>> 1); 1841 r = n - q*dLong; 1842 1843 // Correct the approximation 1844 while (r < 0) { 1845 r += dLong; 1846 q--; 1847 } 1848 while (r >= dLong) { 1849 r -= dLong; 1850 q++; 1851 } 1852 // n - q*dlong == r && 0 <= r <dLong, hence we're done. 1853 return (r << 32) | (q & LONG_MASK); 1854 } 1855 1856 /** 1857 * Calculate GCD of this and b. This and b are changed by the computation. 1858 */ 1859 MutableBigInteger hybridGCD(MutableBigInteger b) { 1860 // Use Euclid's algorithm until the numbers are approximately the 1861 // same length, then use the binary GCD algorithm to find the GCD. 1862 MutableBigInteger a = this; 1863 MutableBigInteger q = new MutableBigInteger(); 1864 1865 while (b.intLen != 0) { 1866 if (Math.abs(a.intLen - b.intLen) < 2) 1867 return a.binaryGCD(b); 1868 1869 MutableBigInteger r = a.divide(b, q); 1870 a = b; 1871 b = r; 1872 } 1873 return a; 1874 } 1875 1876 /** 1877 * Calculate GCD of this and v. 1878 * Assumes that this and v are not zero. 1879 */ 1880 private MutableBigInteger binaryGCD(MutableBigInteger v) { 1881 // Algorithm B from Knuth section 4.5.2 1882 MutableBigInteger u = this; 1883 MutableBigInteger r = new MutableBigInteger(); 1884 1885 // step B1 1886 int s1 = u.getLowestSetBit(); 1887 int s2 = v.getLowestSetBit(); 1888 int k = (s1 < s2) ? s1 : s2; 1889 if (k != 0) { 1890 u.rightShift(k); 1891 v.rightShift(k); 1892 } 1893 1894 // step B2 1895 boolean uOdd = (k==s1); 1896 MutableBigInteger t = uOdd ? v: u; 1897 int tsign = uOdd ? -1 : 1; 1898 1899 int lb; 1900 while ((lb = t.getLowestSetBit()) >= 0) { 1901 // steps B3 and B4 1902 t.rightShift(lb); 1903 // step B5 1904 if (tsign > 0) 1905 u = t; 1906 else 1907 v = t; 1908 1909 // Special case one word numbers 1910 if (u.intLen < 2 && v.intLen < 2) { 1911 int x = u.value[u.offset]; 1912 int y = v.value[v.offset]; 1913 x = binaryGcd(x, y); 1914 r.value[0] = x; 1915 r.intLen = 1; 1916 r.offset = 0; 1917 if (k > 0) 1918 r.leftShift(k); 1919 return r; 1920 } 1921 1922 // step B6 1923 if ((tsign = u.difference(v)) == 0) 1924 break; 1925 t = (tsign >= 0) ? u : v; 1926 } 1927 1928 if (k > 0) 1929 u.leftShift(k); 1930 return u; 1931 } 1932 1933 /** 1934 * Calculate GCD of a and b interpreted as unsigned integers. 1935 */ 1936 static int binaryGcd(int a, int b) { 1937 if (b==0) 1938 return a; 1939 if (a==0) 1940 return b; 1941 1942 // Right shift a & b till their last bits equal to 1. 1943 int aZeros = Integer.numberOfTrailingZeros(a); 1944 int bZeros = Integer.numberOfTrailingZeros(b); 1945 a >>>= aZeros; 1946 b >>>= bZeros; 1947 1948 int t = (aZeros < bZeros ? aZeros : bZeros); 1949 1950 while (a != b) { 1951 if ((a+0x80000000) > (b+0x80000000)) { // a > b as unsigned 1952 a -= b; 1953 a >>>= Integer.numberOfTrailingZeros(a); 1954 } else { 1955 b -= a; 1956 b >>>= Integer.numberOfTrailingZeros(b); 1957 } 1958 } 1959 return a<<t; 1960 } 1961 1962 /** 1963 * Returns the modInverse of this mod p. This and p are not affected by 1964 * the operation. 1965 */ 1966 MutableBigInteger mutableModInverse(MutableBigInteger p) { 1967 // Modulus is odd, use Schroeppel's algorithm 1968 if (p.isOdd()) 1969 return modInverse(p); 1970 1971 // Base and modulus are even, throw exception 1972 if (isEven()) 1973 throw new ArithmeticException("BigInteger not invertible."); 1974 1975 // Get even part of modulus expressed as a power of 2 1976 int powersOf2 = p.getLowestSetBit(); 1977 1978 // Construct odd part of modulus 1979 MutableBigInteger oddMod = new MutableBigInteger(p); 1980 oddMod.rightShift(powersOf2); 1981 1982 if (oddMod.isOne()) 1983 return modInverseMP2(powersOf2); 1984 1985 // Calculate 1/a mod oddMod 1986 MutableBigInteger oddPart = modInverse(oddMod); 1987 1988 // Calculate 1/a mod evenMod 1989 MutableBigInteger evenPart = modInverseMP2(powersOf2); 1990 1991 // Combine the results using Chinese Remainder Theorem 1992 MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2); 1993 MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2); 1994 1995 MutableBigInteger temp1 = new MutableBigInteger(); 1996 MutableBigInteger temp2 = new MutableBigInteger(); 1997 MutableBigInteger result = new MutableBigInteger(); 1998 1999 oddPart.leftShift(powersOf2); 2000 oddPart.multiply(y1, result); 2001 2002 evenPart.multiply(oddMod, temp1); 2003 temp1.multiply(y2, temp2); 2004 2005 result.add(temp2); 2006 return result.divide(p, temp1); 2007 } 2008 2009 /* 2010 * Calculate the multiplicative inverse of this mod 2^k. 2011 */ 2012 MutableBigInteger modInverseMP2(int k) { 2013 if (isEven()) 2014 throw new ArithmeticException("Non-invertible. (GCD != 1)"); 2015 2016 if (k > 64) 2017 return euclidModInverse(k); 2018 2019 int t = inverseMod32(value[offset+intLen-1]); 2020 2021 if (k < 33) { 2022 t = (k == 32 ? t : t & ((1 << k) - 1)); 2023 return new MutableBigInteger(t); 2024 } 2025 2026 long pLong = (value[offset+intLen-1] & LONG_MASK); 2027 if (intLen > 1) 2028 pLong |= ((long)value[offset+intLen-2] << 32); 2029 long tLong = t & LONG_MASK; 2030 tLong = tLong * (2 - pLong * tLong); // 1 more Newton iter step 2031 tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1)); 2032 2033 MutableBigInteger result = new MutableBigInteger(new int[2]); 2034 result.value[0] = (int)(tLong >>> 32); 2035 result.value[1] = (int)tLong; 2036 result.intLen = 2; 2037 result.normalize(); 2038 return result; 2039 } 2040 2041 /** 2042 * Returns the multiplicative inverse of val mod 2^32. Assumes val is odd. 2043 */ 2044 static int inverseMod32(int val) { 2045 // Newton's iteration! 2046 int t = val; 2047 t *= 2 - val*t; 2048 t *= 2 - val*t; 2049 t *= 2 - val*t; 2050 t *= 2 - val*t; 2051 return t; 2052 } 2053 2054 /** 2055 * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd. 2056 */ 2057 static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) { 2058 // Copy the mod to protect original 2059 return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k); 2060 } 2061 2062 /** 2063 * Calculate the multiplicative inverse of this mod mod, where mod is odd. 2064 * This and mod are not changed by the calculation. 2065 * 2066 * This method implements an algorithm due to Richard Schroeppel, that uses 2067 * the same intermediate representation as Montgomery Reduction 2068 * ("Montgomery Form"). The algorithm is described in an unpublished 2069 * manuscript entitled "Fast Modular Reciprocals." 2070 */ 2071 private MutableBigInteger modInverse(MutableBigInteger mod) { 2072 MutableBigInteger p = new MutableBigInteger(mod); 2073 MutableBigInteger f = new MutableBigInteger(this); 2074 MutableBigInteger g = new MutableBigInteger(p); 2075 SignedMutableBigInteger c = new SignedMutableBigInteger(1); 2076 SignedMutableBigInteger d = new SignedMutableBigInteger(); 2077 MutableBigInteger temp = null; 2078 SignedMutableBigInteger sTemp = null; 2079 2080 int k = 0; 2081 // Right shift f k times until odd, left shift d k times 2082 if (f.isEven()) { 2083 int trailingZeros = f.getLowestSetBit(); 2084 f.rightShift(trailingZeros); 2085 d.leftShift(trailingZeros); 2086 k = trailingZeros; 2087 } 2088 2089 // The Almost Inverse Algorithm 2090 while(!f.isOne()) { 2091 // If gcd(f, g) != 1, number is not invertible modulo mod 2092 if (f.isZero()) 2093 throw new ArithmeticException("BigInteger not invertible."); 2094 2095 // If f < g exchange f, g and c, d 2096 if (f.compare(g) < 0) { 2097 temp = f; f = g; g = temp; 2098 sTemp = d; d = c; c = sTemp; 2099 } 2100 2101 // If f == g (mod 4) 2102 if (((f.value[f.offset + f.intLen - 1] ^ 2103 g.value[g.offset + g.intLen - 1]) & 3) == 0) { 2104 f.subtract(g); 2105 c.signedSubtract(d); 2106 } else { // If f != g (mod 4) 2107 f.add(g); 2108 c.signedAdd(d); 2109 } 2110 2111 // Right shift f k times until odd, left shift d k times 2112 int trailingZeros = f.getLowestSetBit(); 2113 f.rightShift(trailingZeros); 2114 d.leftShift(trailingZeros); 2115 k += trailingZeros; 2116 } 2117 2118 while (c.sign < 0) 2119 c.signedAdd(p); 2120 2121 return fixup(c, p, k); 2122 } 2123 2124 /** 2125 * The Fixup Algorithm 2126 * Calculates X such that X = C * 2^(-k) (mod P) 2127 * Assumes C<P and P is odd. 2128 */ 2129 static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p, 2130 int k) { 2131 MutableBigInteger temp = new MutableBigInteger(); 2132 // Set r to the multiplicative inverse of p mod 2^32 2133 int r = -inverseMod32(p.value[p.offset+p.intLen-1]); 2134 2135 for(int i=0, numWords = k >> 5; i<numWords; i++) { 2136 // V = R * c (mod 2^j) 2137 int v = r * c.value[c.offset + c.intLen-1]; 2138 // c = c + (v * p) 2139 p.mul(v, temp); 2140 c.add(temp); 2141 // c = c / 2^j 2142 c.intLen--; 2143 } 2144 int numBits = k & 0x1f; 2145 if (numBits != 0) { 2146 // V = R * c (mod 2^j) 2147 int v = r * c.value[c.offset + c.intLen-1]; 2148 v &= ((1<<numBits) - 1); 2149 // c = c + (v * p) 2150 p.mul(v, temp); 2151 c.add(temp); 2152 // c = c / 2^j 2153 c.rightShift(numBits); 2154 } 2155 2156 // In theory, c may be greater than p at this point (Very rare!) 2157 while (c.compare(p) >= 0) 2158 c.subtract(p); 2159 2160 return c; 2161 } 2162 2163 /** 2164 * Uses the extended Euclidean algorithm to compute the modInverse of base 2165 * mod a modulus that is a power of 2. The modulus is 2^k. 2166 */ 2167 MutableBigInteger euclidModInverse(int k) { 2168 MutableBigInteger b = new MutableBigInteger(1); 2169 b.leftShift(k); 2170 MutableBigInteger mod = new MutableBigInteger(b); 2171 2172 MutableBigInteger a = new MutableBigInteger(this); 2173 MutableBigInteger q = new MutableBigInteger(); 2174 MutableBigInteger r = b.divide(a, q); 2175 2176 MutableBigInteger swapper = b; 2177 // swap b & r 2178 b = r; 2179 r = swapper; 2180 2181 MutableBigInteger t1 = new MutableBigInteger(q); 2182 MutableBigInteger t0 = new MutableBigInteger(1); 2183 MutableBigInteger temp = new MutableBigInteger(); 2184 2185 while (!b.isOne()) { 2186 r = a.divide(b, q); 2187 2188 if (r.intLen == 0) 2189 throw new ArithmeticException("BigInteger not invertible."); 2190 2191 swapper = r; 2192 a = swapper; 2193 2194 if (q.intLen == 1) 2195 t1.mul(q.value[q.offset], temp); 2196 else 2197 q.multiply(t1, temp); 2198 swapper = q; 2199 q = temp; 2200 temp = swapper; 2201 t0.add(q); 2202 2203 if (a.isOne()) 2204 return t0; 2205 2206 r = b.divide(a, q); 2207 2208 if (r.intLen == 0) 2209 throw new ArithmeticException("BigInteger not invertible."); 2210 2211 swapper = b; 2212 b = r; 2213 2214 if (q.intLen == 1) 2215 t0.mul(q.value[q.offset], temp); 2216 else 2217 q.multiply(t0, temp); 2218 swapper = q; q = temp; temp = swapper; 2219 2220 t1.add(q); 2221 } 2222 mod.subtract(t1); 2223 return mod; 2224 } 2225 }