1 /* 2 * Copyright (c) 1999, 2007, 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 * @since 1.3 42 */ 43 44 import java.util.Arrays; 45 46 import static java.math.BigInteger.LONG_MASK; 47 import static java.math.BigDecimal.INFLATED; 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 // Constructors 79 80 /** 81 * The default constructor. An empty MutableBigInteger is created with 82 * a one word capacity. 83 */ 84 MutableBigInteger() { 85 value = new int[1]; 86 intLen = 0; 87 } 88 89 /** 90 * Construct a new MutableBigInteger with a magnitude specified by 91 * the int val. 92 */ 93 MutableBigInteger(int val) { 94 value = new int[1]; 95 intLen = 1; 96 value[0] = val; 97 } 98 99 /** 100 * Construct a new MutableBigInteger with the specified value array 101 * up to the length of the array supplied. 102 */ 103 MutableBigInteger(int[] val) { 104 value = val; 105 intLen = val.length; 106 } 107 108 /** 109 * Construct a new MutableBigInteger with a magnitude equal to the 110 * specified BigInteger. 111 */ 112 MutableBigInteger(BigInteger b) { 113 intLen = b.mag.length; 114 value = Arrays.copyOf(b.mag, intLen); 115 } 116 117 /** 118 * Construct a new MutableBigInteger with a magnitude equal to the 119 * specified MutableBigInteger. 120 */ 121 MutableBigInteger(MutableBigInteger val) { 122 intLen = val.intLen; 123 value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen); 124 } 125 126 /** 127 * Internal helper method to return the magnitude array. The caller is not 128 * supposed to modify the returned array. 129 */ 130 private int[] getMagnitudeArray() { 131 if (offset > 0 || value.length != intLen) 132 return Arrays.copyOfRange(value, offset, offset + intLen); 133 return value; 134 } 135 136 /** 137 * Convert this MutableBigInteger to a long value. The caller has to make 138 * sure this MutableBigInteger can be fit into long. 139 */ 140 private long toLong() { 141 assert (intLen <= 2) : "this MutableBigInteger exceeds the range of long"; 142 if (intLen == 0) 143 return 0; 144 long d = value[offset] & LONG_MASK; 145 return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d; 146 } 147 148 /** 149 * Convert this MutableBigInteger to a BigInteger object. 150 */ 151 BigInteger toBigInteger(int sign) { 152 if (intLen == 0 || sign == 0) 153 return BigInteger.ZERO; 154 return new BigInteger(getMagnitudeArray(), sign); 155 } 156 157 /** 158 * Convert this MutableBigInteger to BigDecimal object with the specified sign 159 * and scale. 160 */ 161 BigDecimal toBigDecimal(int sign, int scale) { 162 if (intLen == 0 || sign == 0) 163 return BigDecimal.valueOf(0, scale); 164 int[] mag = getMagnitudeArray(); 165 int len = mag.length; 166 int d = mag[0]; 167 // If this MutableBigInteger can't be fit into long, we need to 168 // make a BigInteger object for the resultant BigDecimal object. 169 if (len > 2 || (d < 0 && len == 2)) 170 return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0); 171 long v = (len == 2) ? 172 ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) : 173 d & LONG_MASK; 174 return new BigDecimal(null, sign == -1 ? -v : v, scale, 0); 175 } 176 177 /** 178 * Clear out a MutableBigInteger for reuse. 179 */ 180 void clear() { 181 offset = intLen = 0; 182 for (int index=0, n=value.length; index < n; index++) 183 value[index] = 0; 184 } 185 186 /** 187 * Set a MutableBigInteger to zero, removing its offset. 188 */ 189 void reset() { 190 offset = intLen = 0; 191 } 192 193 /** 194 * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1 195 * as this MutableBigInteger is numerically less than, equal to, or 196 * greater than <tt>b</tt>. 197 */ 198 final int compare(MutableBigInteger b) { 199 int blen = b.intLen; 200 if (intLen < blen) 201 return -1; 202 if (intLen > blen) 203 return 1; 204 205 // Add Integer.MIN_VALUE to make the comparison act as unsigned integer 206 // comparison. 207 int[] bval = b.value; 208 for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) { 209 int b1 = value[i] + 0x80000000; 210 int b2 = bval[j] + 0x80000000; 211 if (b1 < b2) 212 return -1; 213 if (b1 > b2) 214 return 1; 215 } 216 return 0; 217 } 218 219 /** 220 * Compare this against half of a MutableBigInteger object (Needed for 221 * remainder tests). 222 * Assumes no leading unnecessary zeros, which holds for results 223 * from divide(). 224 */ 225 final int compareHalf(MutableBigInteger b) { 226 int blen = b.intLen; 227 int len = intLen; 228 if (len <= 0) 229 return blen <=0 ? 0 : -1; 230 if (len > blen) 231 return 1; 232 if (len < blen - 1) 233 return -1; 234 int[] bval = b.value; 235 int bstart = 0; 236 int carry = 0; 237 // Only 2 cases left:len == blen or len == blen - 1 238 if (len != blen) { // len == blen - 1 239 if (bval[bstart] == 1) { 240 ++bstart; 241 carry = 0x80000000; 242 } else 243 return -1; 244 } 245 // compare values with right-shifted values of b, 246 // carrying shifted-out bits across words 247 int[] val = value; 248 for (int i = offset, j = bstart; i < len + offset;) { 249 int bv = bval[j++]; 250 long hb = ((bv >>> 1) + carry) & LONG_MASK; 251 long v = val[i++] & LONG_MASK; 252 if (v != hb) 253 return v < hb ? -1 : 1; 254 carry = (bv & 1) << 31; // carray will be either 0x80000000 or 0 255 } 256 return carry == 0? 0 : -1; 257 } 258 259 /** 260 * Return the index of the lowest set bit in this MutableBigInteger. If the 261 * magnitude of this MutableBigInteger is zero, -1 is returned. 262 */ 263 private final int getLowestSetBit() { 264 if (intLen == 0) 265 return -1; 266 int j, b; 267 for (j=intLen-1; (j>0) && (value[j+offset]==0); j--) 268 ; 269 b = value[j+offset]; 270 if (b==0) 271 return -1; 272 return ((intLen-1-j)<<5) + Integer.numberOfTrailingZeros(b); 273 } 274 275 /** 276 * Return the int in use in this MutableBigInteger at the specified 277 * index. This method is not used because it is not inlined on all 278 * platforms. 279 */ 280 private final int getInt(int index) { 281 return value[offset+index]; 282 } 283 284 /** 285 * Return a long which is equal to the unsigned value of the int in 286 * use in this MutableBigInteger at the specified index. This method is 287 * not used because it is not inlined on all platforms. 288 */ 289 private final long getLong(int index) { 290 return value[offset+index] & LONG_MASK; 291 } 292 293 /** 294 * Ensure that the MutableBigInteger is in normal form, specifically 295 * making sure that there are no leading zeros, and that if the 296 * magnitude is zero, then intLen is zero. 297 */ 298 final void normalize() { 299 if (intLen == 0) { 300 offset = 0; 301 return; 302 } 303 304 int index = offset; 305 if (value[index] != 0) 306 return; 307 308 int indexBound = index+intLen; 309 do { 310 index++; 311 } while(index < indexBound && value[index]==0); 312 313 int numZeros = index - offset; 314 intLen -= numZeros; 315 offset = (intLen==0 ? 0 : offset+numZeros); 316 } 317 318 /** 319 * If this MutableBigInteger cannot hold len words, increase the size 320 * of the value array to len words. 321 */ 322 private final void ensureCapacity(int len) { 323 if (value.length < len) { 324 value = new int[len]; 325 offset = 0; 326 intLen = len; 327 } 328 } 329 330 /** 331 * Convert this MutableBigInteger into an int array with no leading 332 * zeros, of a length that is equal to this MutableBigInteger's intLen. 333 */ 334 int[] toIntArray() { 335 int[] result = new int[intLen]; 336 for(int i=0; i<intLen; i++) 337 result[i] = value[offset+i]; 338 return result; 339 } 340 341 /** 342 * Sets the int at index+offset in this MutableBigInteger to val. 343 * This does not get inlined on all platforms so it is not used 344 * as often as originally intended. 345 */ 346 void setInt(int index, int val) { 347 value[offset + index] = val; 348 } 349 350 /** 351 * Sets this MutableBigInteger's value array to the specified array. 352 * The intLen is set to the specified length. 353 */ 354 void setValue(int[] val, int length) { 355 value = val; 356 intLen = length; 357 offset = 0; 358 } 359 360 /** 361 * Sets this MutableBigInteger's value array to a copy of the specified 362 * array. The intLen is set to the length of the new array. 363 */ 364 void copyValue(MutableBigInteger src) { 365 int len = src.intLen; 366 if (value.length < len) 367 value = new int[len]; 368 System.arraycopy(src.value, src.offset, value, 0, len); 369 intLen = len; 370 offset = 0; 371 } 372 373 /** 374 * Sets this MutableBigInteger's value array to a copy of the specified 375 * array. The intLen is set to the length of the specified array. 376 */ 377 void copyValue(int[] val) { 378 int len = val.length; 379 if (value.length < len) 380 value = new int[len]; 381 System.arraycopy(val, 0, value, 0, len); 382 intLen = len; 383 offset = 0; 384 } 385 386 /** 387 * Returns true iff this MutableBigInteger has a value of one. 388 */ 389 boolean isOne() { 390 return (intLen == 1) && (value[offset] == 1); 391 } 392 393 /** 394 * Returns true iff this MutableBigInteger has a value of zero. 395 */ 396 boolean isZero() { 397 return (intLen == 0); 398 } 399 400 /** 401 * Returns true iff this MutableBigInteger is even. 402 */ 403 boolean isEven() { 404 return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0); 405 } 406 407 /** 408 * Returns true iff this MutableBigInteger is odd. 409 */ 410 boolean isOdd() { 411 return isZero() ? false : ((value[offset + intLen - 1] & 1) == 1); 412 } 413 414 /** 415 * Returns true iff this MutableBigInteger is in normal form. A 416 * MutableBigInteger is in normal form if it has no leading zeros 417 * after the offset, and intLen + offset <= value.length. 418 */ 419 boolean isNormal() { 420 if (intLen + offset > value.length) 421 return false; 422 if (intLen ==0) 423 return true; 424 return (value[offset] != 0); 425 } 426 427 /** 428 * Returns a String representation of this MutableBigInteger in radix 10. 429 */ 430 public String toString() { 431 BigInteger b = toBigInteger(1); 432 return b.toString(); 433 } 434 435 /** 436 * Right shift this MutableBigInteger n bits. The MutableBigInteger is left 437 * in normal form. 438 */ 439 void rightShift(int n) { 440 if (intLen == 0) 441 return; 442 int nInts = n >>> 5; 443 int nBits = n & 0x1F; 444 this.intLen -= nInts; 445 if (nBits == 0) 446 return; 447 int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]); 448 if (nBits >= bitsInHighWord) { 449 this.primitiveLeftShift(32 - nBits); 450 this.intLen--; 451 } else { 452 primitiveRightShift(nBits); 453 } 454 } 455 456 /** 457 * Left shift this MutableBigInteger n bits. 458 */ 459 void leftShift(int n) { 460 /* 461 * If there is enough storage space in this MutableBigInteger already 462 * the available space will be used. Space to the right of the used 463 * ints in the value array is faster to utilize, so the extra space 464 * will be taken from the right if possible. 465 */ 466 if (intLen == 0) 467 return; 468 int nInts = n >>> 5; 469 int nBits = n&0x1F; 470 int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]); 471 472 // If shift can be done without moving words, do so 473 if (n <= (32-bitsInHighWord)) { 474 primitiveLeftShift(nBits); 475 return; 476 } 477 478 int newLen = intLen + nInts +1; 479 if (nBits <= (32-bitsInHighWord)) 480 newLen--; 481 if (value.length < newLen) { 482 // The array must grow 483 int[] result = new int[newLen]; 484 for (int i=0; i<intLen; i++) 485 result[i] = value[offset+i]; 486 setValue(result, newLen); 487 } else if (value.length - offset >= newLen) { 488 // Use space on right 489 for(int i=0; i<newLen - intLen; i++) 490 value[offset+intLen+i] = 0; 491 } else { 492 // Must use space on left 493 for (int i=0; i<intLen; i++) 494 value[i] = value[offset+i]; 495 for (int i=intLen; i<newLen; i++) 496 value[i] = 0; 497 offset = 0; 498 } 499 intLen = newLen; 500 if (nBits == 0) 501 return; 502 if (nBits <= (32-bitsInHighWord)) 503 primitiveLeftShift(nBits); 504 else 505 primitiveRightShift(32 -nBits); 506 } 507 508 /** 509 * A primitive used for division. This method adds in one multiple of the 510 * divisor a back to the dividend result at a specified offset. It is used 511 * when qhat was estimated too large, and must be adjusted. 512 */ 513 private int divadd(int[] a, int[] result, int offset) { 514 long carry = 0; 515 516 for (int j=a.length-1; j >= 0; j--) { 517 long sum = (a[j] & LONG_MASK) + 518 (result[j+offset] & LONG_MASK) + carry; 519 result[j+offset] = (int)sum; 520 carry = sum >>> 32; 521 } 522 return (int)carry; 523 } 524 525 /** 526 * This method is used for division. It multiplies an n word input a by one 527 * word input x, and subtracts the n word product from q. This is needed 528 * when subtracting qhat*divisor from dividend. 529 */ 530 private int mulsub(int[] q, int[] a, int x, int len, int offset) { 531 long xLong = x & LONG_MASK; 532 long carry = 0; 533 offset += len; 534 535 for (int j=len-1; j >= 0; j--) { 536 long product = (a[j] & LONG_MASK) * xLong + carry; 537 long difference = q[offset] - product; 538 q[offset--] = (int)difference; 539 carry = (product >>> 32) 540 + (((difference & LONG_MASK) > 541 (((~(int)product) & LONG_MASK))) ? 1:0); 542 } 543 return (int)carry; 544 } 545 546 /** 547 * Right shift this MutableBigInteger n bits, where n is 548 * less than 32. 549 * Assumes that intLen > 0, n > 0 for speed 550 */ 551 private final void primitiveRightShift(int n) { 552 int[] val = value; 553 int n2 = 32 - n; 554 for (int i=offset+intLen-1, c=val[i]; i>offset; i--) { 555 int b = c; 556 c = val[i-1]; 557 val[i] = (c << n2) | (b >>> n); 558 } 559 val[offset] >>>= n; 560 } 561 562 /** 563 * Left shift this MutableBigInteger n bits, where n is 564 * less than 32. 565 * Assumes that intLen > 0, n > 0 for speed 566 */ 567 private final void primitiveLeftShift(int n) { 568 int[] val = value; 569 int n2 = 32 - n; 570 for (int i=offset, c=val[i], m=i+intLen-1; i<m; i++) { 571 int b = c; 572 c = val[i+1]; 573 val[i] = (b << n) | (c >>> n2); 574 } 575 val[offset+intLen-1] <<= n; 576 } 577 578 /** 579 * Adds the contents of two MutableBigInteger objects.The result 580 * is placed within this MutableBigInteger. 581 * The contents of the addend are not changed. 582 */ 583 void add(MutableBigInteger addend) { 584 int x = intLen; 585 int y = addend.intLen; 586 int resultLen = (intLen > addend.intLen ? intLen : addend.intLen); 587 int[] result = (value.length < resultLen ? new int[resultLen] : value); 588 589 int rstart = result.length-1; 590 long sum; 591 long carry = 0; 592 593 // Add common parts of both numbers 594 while(x>0 && y>0) { 595 x--; y--; 596 sum = (value[x+offset] & LONG_MASK) + 597 (addend.value[y+addend.offset] & LONG_MASK) + carry; 598 result[rstart--] = (int)sum; 599 carry = sum >>> 32; 600 } 601 602 // Add remainder of the longer number 603 while(x>0) { 604 x--; 605 if (carry == 0 && result == value && rstart == (x + offset)) 606 return; 607 sum = (value[x+offset] & LONG_MASK) + carry; 608 result[rstart--] = (int)sum; 609 carry = sum >>> 32; 610 } 611 while(y>0) { 612 y--; 613 sum = (addend.value[y+addend.offset] & LONG_MASK) + carry; 614 result[rstart--] = (int)sum; 615 carry = sum >>> 32; 616 } 617 618 if (carry > 0) { // Result must grow in length 619 resultLen++; 620 if (result.length < resultLen) { 621 int temp[] = new int[resultLen]; 622 // Result one word longer from carry-out; copy low-order 623 // bits into new result. 624 System.arraycopy(result, 0, temp, 1, result.length); 625 temp[0] = 1; 626 result = temp; 627 } else { 628 result[rstart--] = 1; 629 } 630 } 631 632 value = result; 633 intLen = resultLen; 634 offset = result.length - resultLen; 635 } 636 637 638 /** 639 * Subtracts the smaller of this and b from the larger and places the 640 * result into this MutableBigInteger. 641 */ 642 int subtract(MutableBigInteger b) { 643 MutableBigInteger a = this; 644 645 int[] result = value; 646 int sign = a.compare(b); 647 648 if (sign == 0) { 649 reset(); 650 return 0; 651 } 652 if (sign < 0) { 653 MutableBigInteger tmp = a; 654 a = b; 655 b = tmp; 656 } 657 658 int resultLen = a.intLen; 659 if (result.length < resultLen) 660 result = new int[resultLen]; 661 662 long diff = 0; 663 int x = a.intLen; 664 int y = b.intLen; 665 int rstart = result.length - 1; 666 667 // Subtract common parts of both numbers 668 while (y>0) { 669 x--; y--; 670 671 diff = (a.value[x+a.offset] & LONG_MASK) - 672 (b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32)); 673 result[rstart--] = (int)diff; 674 } 675 // Subtract remainder of longer number 676 while (x>0) { 677 x--; 678 diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32)); 679 result[rstart--] = (int)diff; 680 } 681 682 value = result; 683 intLen = resultLen; 684 offset = value.length - resultLen; 685 normalize(); 686 return sign; 687 } 688 689 /** 690 * Subtracts the smaller of a and b from the larger and places the result 691 * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no 692 * operation was performed. 693 */ 694 private int difference(MutableBigInteger b) { 695 MutableBigInteger a = this; 696 int sign = a.compare(b); 697 if (sign ==0) 698 return 0; 699 if (sign < 0) { 700 MutableBigInteger tmp = a; 701 a = b; 702 b = tmp; 703 } 704 705 long diff = 0; 706 int x = a.intLen; 707 int y = b.intLen; 708 709 // Subtract common parts of both numbers 710 while (y>0) { 711 x--; y--; 712 diff = (a.value[a.offset+ x] & LONG_MASK) - 713 (b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32)); 714 a.value[a.offset+x] = (int)diff; 715 } 716 // Subtract remainder of longer number 717 while (x>0) { 718 x--; 719 diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32)); 720 a.value[a.offset+x] = (int)diff; 721 } 722 723 a.normalize(); 724 return sign; 725 } 726 727 /** 728 * Multiply the contents of two MutableBigInteger objects. The result is 729 * placed into MutableBigInteger z. The contents of y are not changed. 730 */ 731 void multiply(MutableBigInteger y, MutableBigInteger z) { 732 int xLen = intLen; 733 int yLen = y.intLen; 734 int newLen = xLen + yLen; 735 736 // Put z into an appropriate state to receive product 737 if (z.value.length < newLen) 738 z.value = new int[newLen]; 739 z.offset = 0; 740 z.intLen = newLen; 741 742 // The first iteration is hoisted out of the loop to avoid extra add 743 long carry = 0; 744 for (int j=yLen-1, k=yLen+xLen-1; j >= 0; j--, k--) { 745 long product = (y.value[j+y.offset] & LONG_MASK) * 746 (value[xLen-1+offset] & LONG_MASK) + carry; 747 z.value[k] = (int)product; 748 carry = product >>> 32; 749 } 750 z.value[xLen-1] = (int)carry; 751 752 // Perform the multiplication word by word 753 for (int i = xLen-2; i >= 0; i--) { 754 carry = 0; 755 for (int j=yLen-1, k=yLen+i; j >= 0; j--, k--) { 756 long product = (y.value[j+y.offset] & LONG_MASK) * 757 (value[i+offset] & LONG_MASK) + 758 (z.value[k] & LONG_MASK) + carry; 759 z.value[k] = (int)product; 760 carry = product >>> 32; 761 } 762 z.value[i] = (int)carry; 763 } 764 765 // Remove leading zeros from product 766 z.normalize(); 767 } 768 769 /** 770 * Multiply the contents of this MutableBigInteger by the word y. The 771 * result is placed into z. 772 */ 773 void mul(int y, MutableBigInteger z) { 774 if (y == 1) { 775 z.copyValue(this); 776 return; 777 } 778 779 if (y == 0) { 780 z.clear(); 781 return; 782 } 783 784 // Perform the multiplication word by word 785 long ylong = y & LONG_MASK; 786 int[] zval = (z.value.length<intLen+1 ? new int[intLen + 1] 787 : z.value); 788 long carry = 0; 789 for (int i = intLen-1; i >= 0; i--) { 790 long product = ylong * (value[i+offset] & LONG_MASK) + carry; 791 zval[i+1] = (int)product; 792 carry = product >>> 32; 793 } 794 795 if (carry == 0) { 796 z.offset = 1; 797 z.intLen = intLen; 798 } else { 799 z.offset = 0; 800 z.intLen = intLen + 1; 801 zval[0] = (int)carry; 802 } 803 z.value = zval; 804 } 805 806 /** 807 * This method is used for division of an n word dividend by a one word 808 * divisor. The quotient is placed into quotient. The one word divisor is 809 * specified by divisor. 810 * 811 * @return the remainder of the division is returned. 812 * 813 */ 814 int divideOneWord(int divisor, MutableBigInteger quotient) { 815 long divisorLong = divisor & LONG_MASK; 816 817 // Special case of one word dividend 818 if (intLen == 1) { 819 long dividendValue = value[offset] & LONG_MASK; 820 int q = (int) (dividendValue / divisorLong); 821 int r = (int) (dividendValue - q * divisorLong); 822 quotient.value[0] = q; 823 quotient.intLen = (q == 0) ? 0 : 1; 824 quotient.offset = 0; 825 return r; 826 } 827 828 if (quotient.value.length < intLen) 829 quotient.value = new int[intLen]; 830 quotient.offset = 0; 831 quotient.intLen = intLen; 832 833 // Normalize the divisor 834 int shift = Integer.numberOfLeadingZeros(divisor); 835 836 int rem = value[offset]; 837 long remLong = rem & LONG_MASK; 838 if (remLong < divisorLong) { 839 quotient.value[0] = 0; 840 } else { 841 quotient.value[0] = (int)(remLong / divisorLong); 842 rem = (int) (remLong - (quotient.value[0] * divisorLong)); 843 remLong = rem & LONG_MASK; 844 } 845 846 int xlen = intLen; 847 int[] qWord = new int[2]; 848 while (--xlen > 0) { 849 long dividendEstimate = (remLong<<32) | 850 (value[offset + intLen - xlen] & LONG_MASK); 851 if (dividendEstimate >= 0) { 852 qWord[0] = (int) (dividendEstimate / divisorLong); 853 qWord[1] = (int) (dividendEstimate - qWord[0] * divisorLong); 854 } else { 855 divWord(qWord, dividendEstimate, divisor); 856 } 857 quotient.value[intLen - xlen] = qWord[0]; 858 rem = qWord[1]; 859 remLong = rem & LONG_MASK; 860 } 861 862 quotient.normalize(); 863 // Unnormalize 864 if (shift > 0) 865 return rem % divisor; 866 else 867 return rem; 868 } 869 870 /** 871 * Calculates the quotient of this div b and places the quotient in the 872 * provided MutableBigInteger objects and the remainder object is returned. 873 * 874 * Uses Algorithm D in Knuth section 4.3.1. 875 * Many optimizations to that algorithm have been adapted from the Colin 876 * Plumb C library. 877 * It special cases one word divisors for speed. The content of b is not 878 * changed. 879 * 880 */ 881 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) { 882 if (b.intLen == 0) 883 throw new ArithmeticException("BigInteger divide by zero"); 884 885 // Dividend is zero 886 if (intLen == 0) { 887 quotient.intLen = quotient.offset; 888 return new MutableBigInteger(); 889 } 890 891 int cmp = compare(b); 892 // Dividend less than divisor 893 if (cmp < 0) { 894 quotient.intLen = quotient.offset = 0; 895 return new MutableBigInteger(this); 896 } 897 // Dividend equal to divisor 898 if (cmp == 0) { 899 quotient.value[0] = quotient.intLen = 1; 900 quotient.offset = 0; 901 return new MutableBigInteger(); 902 } 903 904 quotient.clear(); 905 // Special case one word divisor 906 if (b.intLen == 1) { 907 int r = divideOneWord(b.value[b.offset], quotient); 908 if (r == 0) 909 return new MutableBigInteger(); 910 return new MutableBigInteger(r); 911 } 912 913 // Copy divisor value to protect divisor 914 int[] div = Arrays.copyOfRange(b.value, b.offset, b.offset + b.intLen); 915 return divideMagnitude(div, quotient); 916 } 917 918 /** 919 * Internally used to calculate the quotient of this div v and places the 920 * quotient in the provided MutableBigInteger object and the remainder is 921 * returned. 922 * 923 * @return the remainder of the division will be returned. 924 */ 925 long divide(long v, MutableBigInteger quotient) { 926 if (v == 0) 927 throw new ArithmeticException("BigInteger divide by zero"); 928 929 // Dividend is zero 930 if (intLen == 0) { 931 quotient.intLen = quotient.offset = 0; 932 return 0; 933 } 934 if (v < 0) 935 v = -v; 936 937 int d = (int)(v >>> 32); 938 quotient.clear(); 939 // Special case on word divisor 940 if (d == 0) 941 return divideOneWord((int)v, quotient) & LONG_MASK; 942 else { 943 int[] div = new int[]{ d, (int)(v & LONG_MASK) }; 944 return divideMagnitude(div, quotient).toLong(); 945 } 946 } 947 948 /** 949 * Divide this MutableBigInteger by the divisor represented by its magnitude 950 * array. The quotient will be placed into the provided quotient object & 951 * the remainder object is returned. 952 */ 953 private MutableBigInteger divideMagnitude(int[] divisor, 954 MutableBigInteger quotient) { 955 956 // Remainder starts as dividend with space for a leading zero 957 MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]); 958 System.arraycopy(value, offset, rem.value, 1, intLen); 959 rem.intLen = intLen; 960 rem.offset = 1; 961 962 int nlen = rem.intLen; 963 964 // Set the quotient size 965 int dlen = divisor.length; 966 int limit = nlen - dlen + 1; 967 if (quotient.value.length < limit) { 968 quotient.value = new int[limit]; 969 quotient.offset = 0; 970 } 971 quotient.intLen = limit; 972 int[] q = quotient.value; 973 974 // D1 normalize the divisor 975 int shift = Integer.numberOfLeadingZeros(divisor[0]); 976 if (shift > 0) { 977 // First shift will not grow array 978 BigInteger.primitiveLeftShift(divisor, dlen, shift); 979 // But this one might 980 rem.leftShift(shift); 981 } 982 983 // Must insert leading 0 in rem if its length did not change 984 if (rem.intLen == nlen) { 985 rem.offset = 0; 986 rem.value[0] = 0; 987 rem.intLen++; 988 } 989 990 int dh = divisor[0]; 991 long dhLong = dh & LONG_MASK; 992 int dl = divisor[1]; 993 int[] qWord = new int[2]; 994 995 // D2 Initialize j 996 for(int j=0; j<limit; j++) { 997 // D3 Calculate qhat 998 // estimate qhat 999 int qhat = 0; 1000 int qrem = 0; 1001 boolean skipCorrection = false; 1002 int nh = rem.value[j+rem.offset]; 1003 int nh2 = nh + 0x80000000; 1004 int nm = rem.value[j+1+rem.offset]; 1005 1006 if (nh == dh) { 1007 qhat = ~0; 1008 qrem = nh + nm; 1009 skipCorrection = qrem + 0x80000000 < nh2; 1010 } else { 1011 long nChunk = (((long)nh) << 32) | (nm & LONG_MASK); 1012 if (nChunk >= 0) { 1013 qhat = (int) (nChunk / dhLong); 1014 qrem = (int) (nChunk - (qhat * dhLong)); 1015 } else { 1016 divWord(qWord, nChunk, dh); 1017 qhat = qWord[0]; 1018 qrem = qWord[1]; 1019 } 1020 } 1021 1022 if (qhat == 0) 1023 continue; 1024 1025 if (!skipCorrection) { // Correct qhat 1026 long nl = rem.value[j+2+rem.offset] & LONG_MASK; 1027 long rs = ((qrem & LONG_MASK) << 32) | nl; 1028 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK); 1029 1030 if (unsignedLongCompare(estProduct, rs)) { 1031 qhat--; 1032 qrem = (int)((qrem & LONG_MASK) + dhLong); 1033 if ((qrem & LONG_MASK) >= dhLong) { 1034 estProduct -= (dl & LONG_MASK); 1035 rs = ((qrem & LONG_MASK) << 32) | nl; 1036 if (unsignedLongCompare(estProduct, rs)) 1037 qhat--; 1038 } 1039 } 1040 } 1041 1042 // D4 Multiply and subtract 1043 rem.value[j+rem.offset] = 0; 1044 int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset); 1045 1046 // D5 Test remainder 1047 if (borrow + 0x80000000 > nh2) { 1048 // D6 Add back 1049 divadd(divisor, rem.value, j+1+rem.offset); 1050 qhat--; 1051 } 1052 1053 // Store the quotient digit 1054 q[j] = qhat; 1055 } // D7 loop on j 1056 1057 // D8 Unnormalize 1058 if (shift > 0) 1059 rem.rightShift(shift); 1060 1061 quotient.normalize(); 1062 rem.normalize(); 1063 return rem; 1064 } 1065 1066 /** 1067 * Compare two longs as if they were unsigned. 1068 * Returns true iff one is bigger than two. 1069 */ 1070 private boolean unsignedLongCompare(long one, long two) { 1071 return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE); 1072 } 1073 1074 /** 1075 * This method divides a long quantity by an int to estimate 1076 * qhat for two multi precision numbers. It is used when 1077 * the signed value of n is less than zero. 1078 */ 1079 private void divWord(int[] result, long n, int d) { 1080 long dLong = d & LONG_MASK; 1081 1082 if (dLong == 1) { 1083 result[0] = (int)n; 1084 result[1] = 0; 1085 return; 1086 } 1087 1088 // Approximate the quotient and remainder 1089 long q = (n >>> 1) / (dLong >>> 1); 1090 long r = n - q*dLong; 1091 1092 // Correct the approximation 1093 while (r < 0) { 1094 r += dLong; 1095 q--; 1096 } 1097 while (r >= dLong) { 1098 r -= dLong; 1099 q++; 1100 } 1101 1102 // n - q*dlong == r && 0 <= r <dLong, hence we're done. 1103 result[0] = (int)q; 1104 result[1] = (int)r; 1105 } 1106 1107 /** 1108 * Calculate GCD of this and b. This and b are changed by the computation. 1109 */ 1110 MutableBigInteger hybridGCD(MutableBigInteger b) { 1111 // Use Euclid's algorithm until the numbers are approximately the 1112 // same length, then use the binary GCD algorithm to find the GCD. 1113 MutableBigInteger a = this; 1114 MutableBigInteger q = new MutableBigInteger(); 1115 1116 while (b.intLen != 0) { 1117 if (Math.abs(a.intLen - b.intLen) < 2) 1118 return a.binaryGCD(b); 1119 1120 MutableBigInteger r = a.divide(b, q); 1121 a = b; 1122 b = r; 1123 } 1124 return a; 1125 } 1126 1127 /** 1128 * Calculate GCD of this and v. 1129 * Assumes that this and v are not zero. 1130 */ 1131 private MutableBigInteger binaryGCD(MutableBigInteger v) { 1132 // Algorithm B from Knuth section 4.5.2 1133 MutableBigInteger u = this; 1134 MutableBigInteger r = new MutableBigInteger(); 1135 1136 // step B1 1137 int s1 = u.getLowestSetBit(); 1138 int s2 = v.getLowestSetBit(); 1139 int k = (s1 < s2) ? s1 : s2; 1140 if (k != 0) { 1141 u.rightShift(k); 1142 v.rightShift(k); 1143 } 1144 1145 // step B2 1146 boolean uOdd = (k==s1); 1147 MutableBigInteger t = uOdd ? v: u; 1148 int tsign = uOdd ? -1 : 1; 1149 1150 int lb; 1151 while ((lb = t.getLowestSetBit()) >= 0) { 1152 // steps B3 and B4 1153 t.rightShift(lb); 1154 // step B5 1155 if (tsign > 0) 1156 u = t; 1157 else 1158 v = t; 1159 1160 // Special case one word numbers 1161 if (u.intLen < 2 && v.intLen < 2) { 1162 int x = u.value[u.offset]; 1163 int y = v.value[v.offset]; 1164 x = binaryGcd(x, y); 1165 r.value[0] = x; 1166 r.intLen = 1; 1167 r.offset = 0; 1168 if (k > 0) 1169 r.leftShift(k); 1170 return r; 1171 } 1172 1173 // step B6 1174 if ((tsign = u.difference(v)) == 0) 1175 break; 1176 t = (tsign >= 0) ? u : v; 1177 } 1178 1179 if (k > 0) 1180 u.leftShift(k); 1181 return u; 1182 } 1183 1184 /** 1185 * Calculate GCD of a and b interpreted as unsigned integers. 1186 */ 1187 static int binaryGcd(int a, int b) { 1188 if (b==0) 1189 return a; 1190 if (a==0) 1191 return b; 1192 1193 // Right shift a & b till their last bits equal to 1. 1194 int aZeros = Integer.numberOfTrailingZeros(a); 1195 int bZeros = Integer.numberOfTrailingZeros(b); 1196 a >>>= aZeros; 1197 b >>>= bZeros; 1198 1199 int t = (aZeros < bZeros ? aZeros : bZeros); 1200 1201 while (a != b) { 1202 if ((a+0x80000000) > (b+0x80000000)) { // a > b as unsigned 1203 a -= b; 1204 a >>>= Integer.numberOfTrailingZeros(a); 1205 } else { 1206 b -= a; 1207 b >>>= Integer.numberOfTrailingZeros(b); 1208 } 1209 } 1210 return a<<t; 1211 } 1212 1213 /** 1214 * Returns the modInverse of this mod p. This and p are not affected by 1215 * the operation. 1216 */ 1217 MutableBigInteger mutableModInverse(MutableBigInteger p) { 1218 // Modulus is odd, use Schroeppel's algorithm 1219 if (p.isOdd()) 1220 return modInverse(p); 1221 1222 // Base and modulus are even, throw exception 1223 if (isEven()) 1224 throw new ArithmeticException("BigInteger not invertible."); 1225 1226 // Get even part of modulus expressed as a power of 2 1227 int powersOf2 = p.getLowestSetBit(); 1228 1229 // Construct odd part of modulus 1230 MutableBigInteger oddMod = new MutableBigInteger(p); 1231 oddMod.rightShift(powersOf2); 1232 1233 if (oddMod.isOne()) 1234 return modInverseMP2(powersOf2); 1235 1236 // Calculate 1/a mod oddMod 1237 MutableBigInteger oddPart = modInverse(oddMod); 1238 1239 // Calculate 1/a mod evenMod 1240 MutableBigInteger evenPart = modInverseMP2(powersOf2); 1241 1242 // Combine the results using Chinese Remainder Theorem 1243 MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2); 1244 MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2); 1245 1246 MutableBigInteger temp1 = new MutableBigInteger(); 1247 MutableBigInteger temp2 = new MutableBigInteger(); 1248 MutableBigInteger result = new MutableBigInteger(); 1249 1250 oddPart.leftShift(powersOf2); 1251 oddPart.multiply(y1, result); 1252 1253 evenPart.multiply(oddMod, temp1); 1254 temp1.multiply(y2, temp2); 1255 1256 result.add(temp2); 1257 return result.divide(p, temp1); 1258 } 1259 1260 /* 1261 * Calculate the multiplicative inverse of this mod 2^k. 1262 */ 1263 MutableBigInteger modInverseMP2(int k) { 1264 if (isEven()) 1265 throw new ArithmeticException("Non-invertible. (GCD != 1)"); 1266 1267 if (k > 64) 1268 return euclidModInverse(k); 1269 1270 int t = inverseMod32(value[offset+intLen-1]); 1271 1272 if (k < 33) { 1273 t = (k == 32 ? t : t & ((1 << k) - 1)); 1274 return new MutableBigInteger(t); 1275 } 1276 1277 long pLong = (value[offset+intLen-1] & LONG_MASK); 1278 if (intLen > 1) 1279 pLong |= ((long)value[offset+intLen-2] << 32); 1280 long tLong = t & LONG_MASK; 1281 tLong = tLong * (2 - pLong * tLong); // 1 more Newton iter step 1282 tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1)); 1283 1284 MutableBigInteger result = new MutableBigInteger(new int[2]); 1285 result.value[0] = (int)(tLong >>> 32); 1286 result.value[1] = (int)tLong; 1287 result.intLen = 2; 1288 result.normalize(); 1289 return result; 1290 } 1291 1292 /* 1293 * Returns the multiplicative inverse of val mod 2^32. Assumes val is odd. 1294 */ 1295 static int inverseMod32(int val) { 1296 // Newton's iteration! 1297 int t = val; 1298 t *= 2 - val*t; 1299 t *= 2 - val*t; 1300 t *= 2 - val*t; 1301 t *= 2 - val*t; 1302 return t; 1303 } 1304 1305 /* 1306 * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd. 1307 */ 1308 static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) { 1309 // Copy the mod to protect original 1310 return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k); 1311 } 1312 1313 /** 1314 * Calculate the multiplicative inverse of this mod mod, where mod is odd. 1315 * This and mod are not changed by the calculation. 1316 * 1317 * This method implements an algorithm due to Richard Schroeppel, that uses 1318 * the same intermediate representation as Montgomery Reduction 1319 * ("Montgomery Form"). The algorithm is described in an unpublished 1320 * manuscript entitled "Fast Modular Reciprocals." 1321 */ 1322 private MutableBigInteger modInverse(MutableBigInteger mod) { 1323 MutableBigInteger p = new MutableBigInteger(mod); 1324 MutableBigInteger f = new MutableBigInteger(this); 1325 MutableBigInteger g = new MutableBigInteger(p); 1326 SignedMutableBigInteger c = new SignedMutableBigInteger(1); 1327 SignedMutableBigInteger d = new SignedMutableBigInteger(); 1328 MutableBigInteger temp = null; 1329 SignedMutableBigInteger sTemp = null; 1330 1331 int k = 0; 1332 // Right shift f k times until odd, left shift d k times 1333 if (f.isEven()) { 1334 int trailingZeros = f.getLowestSetBit(); 1335 f.rightShift(trailingZeros); 1336 d.leftShift(trailingZeros); 1337 k = trailingZeros; 1338 } 1339 1340 // The Almost Inverse Algorithm 1341 while(!f.isOne()) { 1342 // If gcd(f, g) != 1, number is not invertible modulo mod 1343 if (f.isZero()) 1344 throw new ArithmeticException("BigInteger not invertible."); 1345 1346 // If f < g exchange f, g and c, d 1347 if (f.compare(g) < 0) { 1348 temp = f; f = g; g = temp; 1349 sTemp = d; d = c; c = sTemp; 1350 } 1351 1352 // If f == g (mod 4) 1353 if (((f.value[f.offset + f.intLen - 1] ^ 1354 g.value[g.offset + g.intLen - 1]) & 3) == 0) { 1355 f.subtract(g); 1356 c.signedSubtract(d); 1357 } else { // If f != g (mod 4) 1358 f.add(g); 1359 c.signedAdd(d); 1360 } 1361 1362 // Right shift f k times until odd, left shift d k times 1363 int trailingZeros = f.getLowestSetBit(); 1364 f.rightShift(trailingZeros); 1365 d.leftShift(trailingZeros); 1366 k += trailingZeros; 1367 } 1368 1369 while (c.sign < 0) 1370 c.signedAdd(p); 1371 1372 return fixup(c, p, k); 1373 } 1374 1375 /* 1376 * The Fixup Algorithm 1377 * Calculates X such that X = C * 2^(-k) (mod P) 1378 * Assumes C<P and P is odd. 1379 */ 1380 static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p, 1381 int k) { 1382 MutableBigInteger temp = new MutableBigInteger(); 1383 // Set r to the multiplicative inverse of p mod 2^32 1384 int r = -inverseMod32(p.value[p.offset+p.intLen-1]); 1385 1386 for(int i=0, numWords = k >> 5; i<numWords; i++) { 1387 // V = R * c (mod 2^j) 1388 int v = r * c.value[c.offset + c.intLen-1]; 1389 // c = c + (v * p) 1390 p.mul(v, temp); 1391 c.add(temp); 1392 // c = c / 2^j 1393 c.intLen--; 1394 } 1395 int numBits = k & 0x1f; 1396 if (numBits != 0) { 1397 // V = R * c (mod 2^j) 1398 int v = r * c.value[c.offset + c.intLen-1]; 1399 v &= ((1<<numBits) - 1); 1400 // c = c + (v * p) 1401 p.mul(v, temp); 1402 c.add(temp); 1403 // c = c / 2^j 1404 c.rightShift(numBits); 1405 } 1406 1407 // In theory, c may be greater than p at this point (Very rare!) 1408 while (c.compare(p) >= 0) 1409 c.subtract(p); 1410 1411 return c; 1412 } 1413 1414 /** 1415 * Uses the extended Euclidean algorithm to compute the modInverse of base 1416 * mod a modulus that is a power of 2. The modulus is 2^k. 1417 */ 1418 MutableBigInteger euclidModInverse(int k) { 1419 MutableBigInteger b = new MutableBigInteger(1); 1420 b.leftShift(k); 1421 MutableBigInteger mod = new MutableBigInteger(b); 1422 1423 MutableBigInteger a = new MutableBigInteger(this); 1424 MutableBigInteger q = new MutableBigInteger(); 1425 MutableBigInteger r = b.divide(a, q); 1426 1427 MutableBigInteger swapper = b; 1428 // swap b & r 1429 b = r; 1430 r = swapper; 1431 1432 MutableBigInteger t1 = new MutableBigInteger(q); 1433 MutableBigInteger t0 = new MutableBigInteger(1); 1434 MutableBigInteger temp = new MutableBigInteger(); 1435 1436 while (!b.isOne()) { 1437 r = a.divide(b, q); 1438 1439 if (r.intLen == 0) 1440 throw new ArithmeticException("BigInteger not invertible."); 1441 1442 swapper = r; 1443 a = swapper; 1444 1445 if (q.intLen == 1) 1446 t1.mul(q.value[q.offset], temp); 1447 else 1448 q.multiply(t1, temp); 1449 swapper = q; 1450 q = temp; 1451 temp = swapper; 1452 t0.add(q); 1453 1454 if (a.isOne()) 1455 return t0; 1456 1457 r = b.divide(a, q); 1458 1459 if (r.intLen == 0) 1460 throw new ArithmeticException("BigInteger not invertible."); 1461 1462 swapper = b; 1463 b = r; 1464 1465 if (q.intLen == 1) 1466 t0.mul(q.value[q.offset], temp); 1467 else 1468 q.multiply(t0, temp); 1469 swapper = q; q = temp; temp = swapper; 1470 1471 t1.add(q); 1472 } 1473 mod.subtract(t1); 1474 return mod; 1475 } 1476 1477 }