src/share/classes/java/math/MutableBigInteger.java

Print this page
rev 7663 : 8014319: Faster division of large integers
Summary: Implement Burnickel-Ziegler division algorithm in BigInteger
Reviewed-by: bpb
Contributed-by: Tim Buktu <tbuktu@hotmail.com>
rev 7664 : 8020641: Clean up some code style in recent BigInteger contributions
Summary: Some minor cleanup to adhere better to Java coding conventions.
Reviewed-by: darcy
Contributed-by: Brian Burkhalter <brian.burkhalter@oracle.com>


 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;


 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 


 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;


 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 


 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)


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      /**


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) {


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     /**


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) {


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             }


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);


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;


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;


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;


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 


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 




 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;


 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     /**
 534      * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
 535      * in normal form.
 536      */
 537     void rightShift(int n) {
 538         if (intLen == 0)
 539             return;
 540         int nInts = n >>> 5;
 541         int nBits = n & 0x1F;
 542         this.intLen -= nInts;
 543         if (nBits == 0)
 544             return;
 545         int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
 546         if (nBits >= bitsInHighWord) {
 547             this.primitiveLeftShift(32 - nBits);
 548             this.intLen--;
 549         } else {
 550             primitiveRightShift(nBits);
 551         }
 552     }
 553 
 554     /**
 555      * Like {@link #leftShift(int)} but {@code n} can be zero.
 556      */
 557     void safeLeftShift(int n) {
 558         if (n > 0) {
 559             leftShift(n);
 560         }
 561     }
 562 
 563     /**
 564      * Left shift this MutableBigInteger n bits.
 565      */
 566     void leftShift(int n) {
 567         /*
 568          * If there is enough storage space in this MutableBigInteger already
 569          * the available space will be used. Space to the right of the used
 570          * ints in the value array is faster to utilize, so the extra space
 571          * will be taken from the right if possible.
 572          */
 573         if (intLen == 0)
 574            return;
 575         int nInts = n >>> 5;
 576         int nBits = n&0x1F;
 577         int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
 578 
 579         // If shift can be done without moving words, do so
 580         if (n <= (32-bitsInHighWord)) {
 581             primitiveLeftShift(nBits);
 582             return;
 583         }
 584 
 585         int newLen = intLen + nInts +1;
 586         if (nBits <= (32-bitsInHighWord))
 587             newLen--;
 588         if (value.length < newLen) {
 589             // The array must grow
 590             int[] result = new int[newLen];
 591             for (int i=0; i < intLen; i++)
 592                 result[i] = value[offset+i];
 593             setValue(result, newLen);
 594         } else if (value.length - offset >= newLen) {
 595             // Use space on right
 596             for(int i=0; i < newLen - intLen; i++)
 597                 value[offset+intLen+i] = 0;
 598         } else {
 599             // Must use space on left
 600             for (int i=0; i < intLen; i++)
 601                 value[i] = value[offset+i];
 602             for (int i=intLen; i < newLen; i++)
 603                 value[i] = 0;
 604             offset = 0;
 605         }
 606         intLen = newLen;
 607         if (nBits == 0)
 608             return;
 609         if (nBits <= (32-bitsInHighWord))
 610             primitiveLeftShift(nBits);
 611         else
 612             primitiveRightShift(32 -nBits);
 613     }
 614 
 615     /**
 616      * A primitive used for division. This method adds in one multiple of the
 617      * divisor a back to the dividend result at a specified offset. It is used
 618      * when qhat was estimated too large, and must be adjusted.
 619      */
 620     private int divadd(int[] a, int[] result, int offset) {
 621         long carry = 0;
 622 


 659         long carry = 0;
 660         offset += len;
 661         for (int j=len-1; j >= 0; j--) {
 662             long product = (a[j] & LONG_MASK) * xLong + carry;
 663             long difference = q[offset--] - product;
 664             carry = (product >>> 32)
 665                      + (((difference & LONG_MASK) >
 666                          (((~(int)product) & LONG_MASK))) ? 1:0);
 667         }
 668         return (int)carry;
 669     }
 670 
 671     /**
 672      * Right shift this MutableBigInteger n bits, where n is
 673      * less than 32.
 674      * Assumes that intLen > 0, n > 0 for speed
 675      */
 676     private final void primitiveRightShift(int n) {
 677         int[] val = value;
 678         int n2 = 32 - n;
 679         for (int i=offset+intLen-1, c=val[i]; i > offset; i--) {
 680             int b = c;
 681             c = val[i-1];
 682             val[i] = (c << n2) | (b >>> n);
 683         }
 684         val[offset] >>>= n;
 685     }
 686 
 687     /**
 688      * Left shift this MutableBigInteger n bits, where n is
 689      * less than 32.
 690      * Assumes that intLen > 0, n > 0 for speed
 691      */
 692     private final void primitiveLeftShift(int n) {
 693         int[] val = value;
 694         int n2 = 32 - n;
 695         for (int i=offset, c=val[i], m=i+intLen-1; i < m; i++) {
 696             int b = c;
 697             c = val[i+1];
 698             val[i] = (b << n) | (c >>> n2);
 699         }
 700         val[offset+intLen-1] <<= n;
 701     }
 702 
 703     /**
 704      * Returns a {@code BigInteger} equal to the {@code n}
 705      * low ints of this number.
 706      */
 707     private BigInteger getLower(int n) {
 708         if (isZero()) {
 709             return BigInteger.ZERO;
 710         } else if (intLen < n) {
 711             return toBigInteger(1);
 712         } else {
 713             // strip zeros
 714             int len = n;
 715             while (len > 0 && value[offset+intLen-len] == 0)
 716                 len--;
 717             int sign = len > 0 ? 1 : 0;
 718             return new BigInteger(Arrays.copyOfRange(value, offset+intLen-len, offset+intLen), sign);
 719         }
 720     }
 721 
 722     /**
 723      * Discards all ints whose index is greater than {@code n}.
 724      */
 725     private void keepLower(int n) {
 726         if (intLen >= n) {
 727             offset += intLen - n;
 728             intLen = n;
 729         }
 730     }
 731 
 732     /**
 733      * Adds the contents of two MutableBigInteger objects.The result
 734      * is placed within this MutableBigInteger.
 735      * The contents of the addend are not changed.
 736      */
 737     void add(MutableBigInteger addend) {
 738         int x = intLen;
 739         int y = addend.intLen;
 740         int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
 741         int[] result = (value.length < resultLen ? new int[resultLen] : value);
 742 
 743         int rstart = result.length-1;
 744         long sum;
 745         long carry = 0;
 746 
 747         // Add common parts of both numbers
 748         while(x > 0 && y > 0) {
 749             x--; y--;
 750             sum = (value[x+offset] & LONG_MASK) +
 751                 (addend.value[y+addend.offset] & LONG_MASK) + carry;
 752             result[rstart--] = (int)sum;
 753             carry = sum >>> 32;
 754         }
 755 
 756         // Add remainder of the longer number
 757         while(x > 0) {
 758             x--;
 759             if (carry == 0 && result == value && rstart == (x + offset))
 760                 return;
 761             sum = (value[x+offset] & LONG_MASK) + carry;
 762             result[rstart--] = (int)sum;
 763             carry = sum >>> 32;
 764         }
 765         while(y > 0) {
 766             y--;
 767             sum = (addend.value[y+addend.offset] & LONG_MASK) + carry;
 768             result[rstart--] = (int)sum;
 769             carry = sum >>> 32;
 770         }
 771 
 772         if (carry > 0) { // Result must grow in length
 773             resultLen++;
 774             if (result.length < resultLen) {
 775                 int temp[] = new int[resultLen];
 776                 // Result one word longer from carry-out; copy low-order
 777                 // bits into new result.
 778                 System.arraycopy(result, 0, temp, 1, result.length);
 779                 temp[0] = 1;
 780                 result = temp;
 781             } else {
 782                 result[rstart--] = 1;
 783             }
 784         }
 785 
 786         value = result;
 787         intLen = resultLen;
 788         offset = result.length - resultLen;
 789     }
 790 
 791     /**
 792      * Adds the value of {@code addend} shifted {@code n} ints to the left.
 793      * Has the same effect as {@code addend.leftShift(32*ints); add(addend);}
 794      * but doesn't change the value of {@code addend}.
 795      */
 796     void addShifted(MutableBigInteger addend, int n) {
 797         if (addend.isZero()) {
 798             return;
 799         }
 800 
 801         int x = intLen;
 802         int y = addend.intLen + n;
 803         int resultLen = (intLen > y ? intLen : y);
 804         int[] result = (value.length < resultLen ? new int[resultLen] : value);
 805 
 806         int rstart = result.length-1;
 807         long sum;
 808         long carry = 0;
 809 
 810         // Add common parts of both numbers
 811         while (x > 0 && y > 0) {
 812             x--; y--;
 813             int bval = y+addend.offset < addend.value.length ? addend.value[y+addend.offset] : 0;
 814             sum = (value[x+offset] & LONG_MASK) +
 815                 (bval & LONG_MASK) + carry;
 816             result[rstart--] = (int)sum;
 817             carry = sum >>> 32;
 818         }
 819 
 820         // Add remainder of the longer number
 821         while (x > 0) {
 822             x--;
 823             if (carry == 0 && result == value && rstart == (x + offset)) {
 824                 return;
 825             }
 826             sum = (value[x+offset] & LONG_MASK) + carry;
 827             result[rstart--] = (int)sum;
 828             carry = sum >>> 32;
 829         }
 830         while (y > 0) {
 831             y--;
 832             int bval = y+addend.offset < addend.value.length ? addend.value[y+addend.offset] : 0;
 833             sum = (bval & LONG_MASK) + carry;
 834             result[rstart--] = (int)sum;
 835             carry = sum >>> 32;
 836         }
 837 
 838         if (carry > 0) { // Result must grow in length
 839             resultLen++;
 840             if (result.length < resultLen) {
 841                 int temp[] = new int[resultLen];
 842                 // Result one word longer from carry-out; copy low-order
 843                 // bits into new result.
 844                 System.arraycopy(result, 0, temp, 1, result.length);
 845                 temp[0] = 1;
 846                 result = temp;
 847             } else {
 848                 result[rstart--] = 1;
 849             }
 850         }
 851 
 852         value = result;


 868         int resultLen = (intLen > y ? intLen : y);
 869         int[] result;
 870         if (value.length < resultLen)
 871             result = new int[resultLen];
 872         else {
 873             result = value;
 874             Arrays.fill(value, offset+intLen, value.length, 0);
 875         }
 876 
 877         int rstart = result.length-1;
 878 
 879         // copy from this if needed
 880         System.arraycopy(value, offset, result, rstart+1-x, x);
 881         y -= x;
 882         rstart -= x;
 883 
 884         int len = Math.min(y, addend.value.length-addend.offset);
 885         System.arraycopy(addend.value, addend.offset, result, rstart+1-y, len);
 886 
 887         // zero the gap
 888         for (int i=rstart+1-y+len; i < rstart+1; i++)
 889             result[i] = 0;
 890 
 891         value = result;
 892         intLen = resultLen;
 893         offset = result.length - resultLen;
 894     }
 895 
 896     /**
 897      * Adds the low {@code n} ints of {@code addend}.
 898      */
 899     void addLower(MutableBigInteger addend, int n) {
 900         MutableBigInteger a = new MutableBigInteger(addend);
 901         if (a.offset + a.intLen >= n) {
 902             a.offset = a.offset + a.intLen - n;
 903             a.intLen = n;
 904         }
 905         a.normalize();
 906         add(a);
 907     }
 908 


 919         if (sign == 0) {
 920             reset();
 921             return 0;
 922         }
 923         if (sign < 0) {
 924             MutableBigInteger tmp = a;
 925             a = b;
 926             b = tmp;
 927         }
 928 
 929         int resultLen = a.intLen;
 930         if (result.length < resultLen)
 931             result = new int[resultLen];
 932 
 933         long diff = 0;
 934         int x = a.intLen;
 935         int y = b.intLen;
 936         int rstart = result.length - 1;
 937 
 938         // Subtract common parts of both numbers
 939         while (y > 0) {
 940             x--; y--;
 941 
 942             diff = (a.value[x+a.offset] & LONG_MASK) -
 943                    (b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32));
 944             result[rstart--] = (int)diff;
 945         }
 946         // Subtract remainder of longer number
 947         while (x > 0) {
 948             x--;
 949             diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32));
 950             result[rstart--] = (int)diff;
 951         }
 952 
 953         value = result;
 954         intLen = resultLen;
 955         offset = value.length - resultLen;
 956         normalize();
 957         return sign;
 958     }
 959 
 960     /**
 961      * Subtracts the smaller of a and b from the larger and places the result
 962      * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
 963      * operation was performed.
 964      */
 965     private int difference(MutableBigInteger b) {
 966         MutableBigInteger a = this;
 967         int sign = a.compare(b);
 968         if (sign == 0)
 969             return 0;
 970         if (sign < 0) {
 971             MutableBigInteger tmp = a;
 972             a = b;
 973             b = tmp;
 974         }
 975 
 976         long diff = 0;
 977         int x = a.intLen;
 978         int y = b.intLen;
 979 
 980         // Subtract common parts of both numbers
 981         while (y > 0) {
 982             x--; y--;
 983             diff = (a.value[a.offset+ x] & LONG_MASK) -
 984                 (b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32));
 985             a.value[a.offset+x] = (int)diff;
 986         }
 987         // Subtract remainder of longer number
 988         while (x > 0) {
 989             x--;
 990             diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32));
 991             a.value[a.offset+x] = (int)diff;
 992         }
 993 
 994         a.normalize();
 995         return sign;
 996     }
 997 
 998     /**
 999      * Multiply the contents of two MutableBigInteger objects. The result is
1000      * placed into MutableBigInteger z. The contents of y are not changed.
1001      */
1002     void multiply(MutableBigInteger y, MutableBigInteger z) {
1003         int xLen = intLen;
1004         int yLen = y.intLen;
1005         int newLen = xLen + yLen;
1006 
1007         // Put z into an appropriate state to receive product
1008         if (z.value.length < newLen)


1037         z.normalize();
1038     }
1039 
1040     /**
1041      * Multiply the contents of this MutableBigInteger by the word y. The
1042      * result is placed into z.
1043      */
1044     void mul(int y, MutableBigInteger z) {
1045         if (y == 1) {
1046             z.copyValue(this);
1047             return;
1048         }
1049 
1050         if (y == 0) {
1051             z.clear();
1052             return;
1053         }
1054 
1055         // Perform the multiplication word by word
1056         long ylong = y & LONG_MASK;
1057         int[] zval = (z.value.length < intLen+1 ? new int[intLen + 1]
1058                                               : z.value);
1059         long carry = 0;
1060         for (int i = intLen-1; i >= 0; i--) {
1061             long product = ylong * (value[i+offset] & LONG_MASK) + carry;
1062             zval[i+1] = (int)product;
1063             carry = product >>> 32;
1064         }
1065 
1066         if (carry == 0) {
1067             z.offset = 1;
1068             z.intLen = intLen;
1069         } else {
1070             z.offset = 0;
1071             z.intLen = intLen + 1;
1072             zval[0] = (int)carry;
1073         }
1074         z.value = zval;
1075     }
1076 
1077      /**


1131         }
1132 
1133         quotient.normalize();
1134         // Unnormalize
1135         if (shift > 0)
1136             return rem % divisor;
1137         else
1138             return rem;
1139     }
1140 
1141     /**
1142      * Calculates the quotient of this div b and places the quotient in the
1143      * provided MutableBigInteger objects and the remainder object is returned.
1144      *
1145      */
1146     MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
1147         return divide(b,quotient,true);
1148     }
1149 
1150     MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
1151         if (intLen < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD ||
1152                 b.intLen < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD) {
1153             return divideKnuth(b, quotient, needRemainder);
1154         } else {
1155             return divideAndRemainderBurnikelZiegler(b, quotient);
1156         }
1157     }
1158 
1159     /**
1160      * @see #divideKnuth(MutableBigInteger, MutableBigInteger, boolean)
1161      */
1162     MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient) {
1163         return divideKnuth(b,quotient,true);
1164     }
1165 
1166     /**
1167      * Calculates the quotient of this div b and places the quotient in the
1168      * provided MutableBigInteger objects and the remainder object is returned.
1169      *
1170      * Uses Algorithm D in Knuth section 4.3.1.
1171      * Many optimizations to that algorithm have been adapted from the Colin
1172      * Plumb C library.
1173      * It special cases one word divisors for speed. The content of b is not
1174      * changed.
1175      *
1176      */
1177     MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {


1225         }
1226 
1227         return divideMagnitude(b, quotient, needRemainder);
1228     }
1229 
1230     /**
1231      * Computes {@code this/b} and {@code this%b} using the
1232      * <a href="http://cr.yp.to/bib/1998/burnikel.ps"> Burnikel-Ziegler algorithm</a>.
1233      * This method implements algorithm 3 from pg. 9 of the Burnikel-Ziegler paper.
1234      * The parameter beta was chosen to b 2<sup>32</sup> so almost all shifts are
1235      * multiples of 32 bits.<br/>
1236      * {@code this} and {@code b} must be nonnegative.
1237      * @param b the divisor
1238      * @param quotient output parameter for {@code this/b}
1239      * @return the remainder
1240      */
1241     MutableBigInteger divideAndRemainderBurnikelZiegler(MutableBigInteger b, MutableBigInteger quotient) {
1242         int r = intLen;
1243         int s = b.intLen;
1244 
1245         if (r < s) {
1246             return this;
1247         } else {
1248             // Unlike Knuth division, we don't check for common powers of two here because
1249             // BZ already runs faster if both numbers contain powers of two and cancelling them has no
1250             // additional benefit.
1251 
1252             // step 1: let m = min{2^k | (2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s}
1253             int m = 1 << (32-Integer.numberOfLeadingZeros(s/BigInteger.BURNIKEL_ZIEGLER_THRESHOLD));
1254 
1255             int j = (s+m-1) / m;      // step 2a: j = ceil(s/m)
1256             int n = j * m;            // step 2b: block length in 32-bit units
1257             int n32 = 32 * n;         // block length in bits
1258             int sigma = Math.max(0, n32 - b.bitLength());   // step 3: sigma = max{T | (2^T)*B < beta^n}
1259             MutableBigInteger bShifted = new MutableBigInteger(b);
1260             bShifted.safeLeftShift(sigma);   // step 4a: shift b so its length is a multiple of n
1261             safeLeftShift(sigma);     // step 4b: shift this by the same amount
1262 
1263             // step 5: t is the number of blocks needed to accommodate this plus one additional bit
1264             int t = (bitLength()+n32) / n32;
1265             if (t < 2) {
1266                 t = 2;
1267             }
1268 
1269             // step 6: conceptually split this into blocks a[t-1], ..., a[0]
1270             MutableBigInteger a1 = getBlock(t-1, t, n);   // the most significant block of this
1271 
1272             // step 7: z[t-2] = [a[t-1], a[t-2]]
1273             MutableBigInteger z = getBlock(t-2, t, n);    // the second to most significant block
1274             z.addDisjoint(a1, n);   // z[t-2]
1275 
1276             // do schoolbook division on blocks, dividing 2-block numbers by 1-block numbers
1277             MutableBigInteger qi = new MutableBigInteger();
1278             MutableBigInteger ri;
1279             quotient.offset = quotient.intLen = 0;
1280             for (int i=t-2; i > 0; i--) {
1281                 // step 8a: compute (qi,ri) such that z=b*qi+ri
1282                 ri = z.divide2n1n(bShifted, qi);
1283 
1284                 // step 8b: z = [ri, a[i-1]]
1285                 z = getBlock(i-1, t, n);   // a[i-1]
1286                 z.addDisjoint(ri, n);
1287                 quotient.addShifted(qi, i*n);   // update q (part of step 9)
1288             }
1289             // final iteration of step 8: do the loop one more time for i=0 but leave z unchanged
1290             ri = z.divide2n1n(bShifted, qi);
1291             quotient.add(qi);
1292 
1293             ri.rightShift(sigma);   // step 9: this and b were shifted, so shift back
1294             return ri;
1295         }
1296     }
1297 
1298     /**
1299      * This method implements algorithm 1 from pg. 4 of the Burnikel-Ziegler paper.
1300      * It divides a 2n-digit number by a n-digit number.<br/>
1301      * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.
1302      * <br/>
1303      * {@code this} must be a nonnegative number such that {@code this.bitLength() <= 2*b.bitLength()}
1304      * @param b a positive number such that {@code b.bitLength()} is even
1305      * @param quotient output parameter for {@code this/b}
1306      * @return {@code this%b}
1307      */
1308     private MutableBigInteger divide2n1n(MutableBigInteger b, MutableBigInteger quotient) {
1309         int n = b.intLen;
1310 
1311         // step 1: base case
1312         if (n%2 != 0 || n < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD) {
1313             return divideKnuth(b, quotient);
1314         }
1315 
1316         // step 2: view this as [a1,a2,a3,a4] where each ai is n/2 ints or less
1317         MutableBigInteger aUpper = new MutableBigInteger(this);
1318         aUpper.safeRightShift(32*(n/2));   // aUpper = [a1,a2,a3]
1319         keepLower(n/2);   // this = a4
1320 
1321         // step 3: q1=aUpper/b, r1=aUpper%b
1322         MutableBigInteger q1 = new MutableBigInteger();
1323         MutableBigInteger r1 = aUpper.divide3n2n(b, q1);
1324 
1325         // step 4: quotient=[r1,this]/b, r2=[r1,this]%b
1326         addDisjoint(r1, n/2);   // this = [r1,this]
1327         MutableBigInteger r2 = divide3n2n(b, quotient);
1328 
1329         // step 5: let quotient=[q1,quotient] and return r2
1330         quotient.addDisjoint(q1, n/2);
1331         return r2;
1332     }
1333 
1334     /**


1343     private MutableBigInteger divide3n2n(MutableBigInteger b, MutableBigInteger quotient) {
1344         int n = b.intLen / 2;   // half the length of b in ints
1345 
1346         // step 1: view this as [a1,a2,a3] where each ai is n ints or less; let a12=[a1,a2]
1347         MutableBigInteger a12 = new MutableBigInteger(this);
1348         a12.safeRightShift(32*n);
1349 
1350         // step 2: view b as [b1,b2] where each bi is n ints or less
1351         MutableBigInteger b1 = new MutableBigInteger(b);
1352         b1.safeRightShift(n * 32);
1353         BigInteger b2 = b.getLower(n);
1354 
1355         MutableBigInteger r;
1356         MutableBigInteger d;
1357         if (compareShifted(b, n) < 0) {
1358             // step 3a: if a1<b1, let quotient=a12/b1 and r=a12%b1
1359             r = a12.divide2n1n(b1, quotient);
1360 
1361             // step 4: d=quotient*b2
1362             d = new MutableBigInteger(quotient.toBigInteger().multiply(b2));
1363         } else {

1364             // step 3b: if a1>=b1, let quotient=beta^n-1 and r=a12-b1*2^n+b1
1365             quotient.ones(n);
1366             a12.add(b1);
1367             b1.leftShift(32*n);
1368             a12.subtract(b1);
1369             r = a12;
1370 
1371             // step 4: d=quotient*b2=(b2 << 32*n) - b2
1372             d = new MutableBigInteger(b2);
1373             d.leftShift(32 * n);
1374             d.subtract(new MutableBigInteger(b2));
1375         }
1376 
1377         // step 5: r = r*beta^n + a3 - d (paper says a4)
1378         // However, don't subtract d until after the while loop so r doesn't become negative
1379         r.leftShift(32 * n);
1380         r.addLower(this, n);
1381 
1382         // step 6: add b until r>=d
1383         while (r.compare(d) < 0) {
1384             r.add(b);
1385             quotient.subtract(MutableBigInteger.ONE);
1386         }
1387         r.subtract(d);
1388 
1389         return r;
1390     }
1391 
1392     /**
1393      * Returns a {@code MutableBigInteger} containing {@code blockLength} ints from
1394      * {@code this} number, starting at {@code index*blockLength}.<br/>
1395      * Used by Burnikel-Ziegler division.
1396      * @param index the block index
1397      * @param numBlocks the total number of blocks in {@code this} number
1398      * @param blockLength length of one block in units of 32 bits
1399      * @return
1400      */
1401     private MutableBigInteger getBlock(int index, int numBlocks, int blockLength) {
1402         int blockStart = index * blockLength;
1403         if (blockStart >= intLen) {
1404             return new MutableBigInteger();
1405         }
1406 
1407         int blockEnd;
1408         if (index == numBlocks-1) {
1409             blockEnd = intLen;
1410         } else {
1411             blockEnd = (index+1) * blockLength;
1412         }
1413         if (blockEnd > intLen) {
1414             return new MutableBigInteger();
1415         }
1416 
1417         int[] newVal = Arrays.copyOfRange(value, offset+intLen-blockEnd, offset+intLen-blockStart);
1418         return new MutableBigInteger(newVal);
1419     }
1420 
1421     /** @see BigInteger#bitLength() */
1422     int bitLength() {
1423         if (intLen == 0)
1424             return 0;
1425         return intLen*32 - Integer.numberOfLeadingZeros(value[offset]);
1426     }
1427 
1428     /**
1429      * Internally used  to calculate the quotient of this div v and places the
1430      * quotient in the provided MutableBigInteger object and the remainder is
1431      * returned.
1432      *
1433      * @return the remainder of the division will be returned.
1434      */
1435     long divide(long v, MutableBigInteger quotient) {


1466     }
1467 
1468     /**
1469      * Divide this MutableBigInteger by the divisor.
1470      * The quotient will be placed into the provided quotient object &
1471      * the remainder object is returned.
1472      */
1473     private MutableBigInteger divideMagnitude(MutableBigInteger div,
1474                                               MutableBigInteger quotient,
1475                                               boolean needRemainder ) {
1476         // assert div.intLen > 1
1477         // D1 normalize the divisor
1478         int shift = Integer.numberOfLeadingZeros(div.value[div.offset]);
1479         // Copy divisor value to protect divisor
1480         final int dlen = div.intLen;
1481         int[] divisor;
1482         MutableBigInteger rem; // Remainder starts as dividend with space for a leading zero
1483         if (shift > 0) {
1484             divisor = new int[dlen];
1485             copyAndShift(div.value,div.offset,dlen,divisor,0,shift);
1486             if (Integer.numberOfLeadingZeros(value[offset]) >= shift) {
1487                 int[] remarr = new int[intLen + 1];
1488                 rem = new MutableBigInteger(remarr);
1489                 rem.intLen = intLen;
1490                 rem.offset = 1;
1491                 copyAndShift(value,offset,intLen,remarr,1,shift);
1492             } else {
1493                 int[] remarr = new int[intLen + 2];
1494                 rem = new MutableBigInteger(remarr);
1495                 rem.intLen = intLen+1;
1496                 rem.offset = 1;
1497                 int rFrom = offset;
1498                 int c=0;
1499                 int n2 = 32 - shift;
1500                 for (int i=1; i < intLen+1; i++,rFrom++) {
1501                     int b = c;
1502                     c = value[rFrom];
1503                     remarr[i] = (b << shift) | (c >>> n2);
1504                 }
1505                 remarr[intLen+1] = c << shift;
1506             }


1519         if (quotient.value.length < limit) {
1520             quotient.value = new int[limit];
1521             quotient.offset = 0;
1522         }
1523         quotient.intLen = limit;
1524         int[] q = quotient.value;
1525 
1526 
1527         // Must insert leading 0 in rem if its length did not change
1528         if (rem.intLen == nlen) {
1529             rem.offset = 0;
1530             rem.value[0] = 0;
1531             rem.intLen++;
1532         }
1533 
1534         int dh = divisor[0];
1535         long dhLong = dh & LONG_MASK;
1536         int dl = divisor[1];
1537 
1538         // D2 Initialize j
1539         for (int j=0; j < limit-1; j++) {
1540             // D3 Calculate qhat
1541             // estimate qhat
1542             int qhat = 0;
1543             int qrem = 0;
1544             boolean skipCorrection = false;
1545             int nh = rem.value[j+rem.offset];
1546             int nh2 = nh + 0x80000000;
1547             int nm = rem.value[j+1+rem.offset];
1548 
1549             if (nh == dh) {
1550                 qhat = ~0;
1551                 qrem = nh + nm;
1552                 skipCorrection = qrem + 0x80000000 < nh2;
1553             } else {
1554                 long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
1555                 if (nChunk >= 0) {
1556                     qhat = (int) (nChunk / dhLong);
1557                     qrem = (int) (nChunk - (qhat * dhLong));
1558                 } else {
1559                     long tmp = divWord(nChunk, dh);


1643             int borrow;
1644             rem.value[limit - 1 + rem.offset] = 0;
1645             if(needRemainder)
1646                 borrow = mulsub(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
1647             else
1648                 borrow = mulsubBorrow(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
1649 
1650             // D5 Test remainder
1651             if (borrow + 0x80000000 > nh2) {
1652                 // D6 Add back
1653                 if(needRemainder)
1654                     divadd(divisor, rem.value, limit - 1 + 1 + rem.offset);
1655                 qhat--;
1656             }
1657 
1658             // Store the quotient digit
1659             q[(limit - 1)] = qhat;
1660         }
1661 
1662 
1663         if (needRemainder) {
1664             // D8 Unnormalize
1665             if (shift > 0)
1666                 rem.rightShift(shift);
1667             rem.normalize();
1668         }
1669         quotient.normalize();
1670         return needRemainder ? rem : null;
1671     }
1672 
1673     /**
1674      * Divide this MutableBigInteger by the divisor represented by positive long
1675      * value. The quotient will be placed into the provided quotient object &
1676      * the remainder object is returned.
1677      */
1678     private MutableBigInteger divideLongMagnitude(long ldivisor, MutableBigInteger quotient) {
1679         // Remainder starts as dividend with space for a leading zero
1680         MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
1681         System.arraycopy(value, offset, rem.value, 1, intLen);
1682         rem.intLen = intLen;
1683         rem.offset = 1;


1885 
1886     /**
1887      * Calculate GCD of this and v.
1888      * Assumes that this and v are not zero.
1889      */
1890     private MutableBigInteger binaryGCD(MutableBigInteger v) {
1891         // Algorithm B from Knuth section 4.5.2
1892         MutableBigInteger u = this;
1893         MutableBigInteger r = new MutableBigInteger();
1894 
1895         // step B1
1896         int s1 = u.getLowestSetBit();
1897         int s2 = v.getLowestSetBit();
1898         int k = (s1 < s2) ? s1 : s2;
1899         if (k != 0) {
1900             u.rightShift(k);
1901             v.rightShift(k);
1902         }
1903 
1904         // step B2
1905         boolean uOdd = (k == s1);
1906         MutableBigInteger t = uOdd ? v: u;
1907         int tsign = uOdd ? -1 : 1;
1908 
1909         int lb;
1910         while ((lb = t.getLowestSetBit()) >= 0) {
1911             // steps B3 and B4
1912             t.rightShift(lb);
1913             // step B5
1914             if (tsign > 0)
1915                 u = t;
1916             else
1917                 v = t;
1918 
1919             // Special case one word numbers
1920             if (u.intLen < 2 && v.intLen < 2) {
1921                 int x = u.value[u.offset];
1922                 int y = v.value[v.offset];
1923                 x  = binaryGcd(x, y);
1924                 r.value[0] = x;
1925                 r.intLen = 1;


1927                 if (k > 0)
1928                     r.leftShift(k);
1929                 return r;
1930             }
1931 
1932             // step B6
1933             if ((tsign = u.difference(v)) == 0)
1934                 break;
1935             t = (tsign >= 0) ? u : v;
1936         }
1937 
1938         if (k > 0)
1939             u.leftShift(k);
1940         return u;
1941     }
1942 
1943     /**
1944      * Calculate GCD of a and b interpreted as unsigned integers.
1945      */
1946     static int binaryGcd(int a, int b) {
1947         if (b == 0)
1948             return a;
1949         if (a == 0)
1950             return b;
1951 
1952         // Right shift a & b till their last bits equal to 1.
1953         int aZeros = Integer.numberOfTrailingZeros(a);
1954         int bZeros = Integer.numberOfTrailingZeros(b);
1955         a >>>= aZeros;
1956         b >>>= bZeros;
1957 
1958         int t = (aZeros < bZeros ? aZeros : bZeros);
1959 
1960         while (a != b) {
1961             if ((a+0x80000000) > (b+0x80000000)) {  // a > b as unsigned
1962                 a -= b;
1963                 a >>>= Integer.numberOfTrailingZeros(a);
1964             } else {
1965                 b -= a;
1966                 b >>>= Integer.numberOfTrailingZeros(b);
1967             }
1968         }
1969         return a<<t;


2080      */
2081     private MutableBigInteger modInverse(MutableBigInteger mod) {
2082         MutableBigInteger p = new MutableBigInteger(mod);
2083         MutableBigInteger f = new MutableBigInteger(this);
2084         MutableBigInteger g = new MutableBigInteger(p);
2085         SignedMutableBigInteger c = new SignedMutableBigInteger(1);
2086         SignedMutableBigInteger d = new SignedMutableBigInteger();
2087         MutableBigInteger temp = null;
2088         SignedMutableBigInteger sTemp = null;
2089 
2090         int k = 0;
2091         // Right shift f k times until odd, left shift d k times
2092         if (f.isEven()) {
2093             int trailingZeros = f.getLowestSetBit();
2094             f.rightShift(trailingZeros);
2095             d.leftShift(trailingZeros);
2096             k = trailingZeros;
2097         }
2098 
2099         // The Almost Inverse Algorithm
2100         while (!f.isOne()) {
2101             // If gcd(f, g) != 1, number is not invertible modulo mod
2102             if (f.isZero())
2103                 throw new ArithmeticException("BigInteger not invertible.");
2104 
2105             // If f < g exchange f, g and c, d
2106             if (f.compare(g) < 0) {
2107                 temp = f; f = g; g = temp;
2108                 sTemp = d; d = c; c = sTemp;
2109             }
2110 
2111             // If f == g (mod 4)
2112             if (((f.value[f.offset + f.intLen - 1] ^
2113                  g.value[g.offset + g.intLen - 1]) & 3) == 0) {
2114                 f.subtract(g);
2115                 c.signedSubtract(d);
2116             } else { // If f != g (mod 4)
2117                 f.add(g);
2118                 c.signedAdd(d);
2119             }
2120 


2125             k += trailingZeros;
2126         }
2127 
2128         while (c.sign < 0)
2129            c.signedAdd(p);
2130 
2131         return fixup(c, p, k);
2132     }
2133 
2134     /**
2135      * The Fixup Algorithm
2136      * Calculates X such that X = C * 2^(-k) (mod P)
2137      * Assumes C<P and P is odd.
2138      */
2139     static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,
2140                                                                       int k) {
2141         MutableBigInteger temp = new MutableBigInteger();
2142         // Set r to the multiplicative inverse of p mod 2^32
2143         int r = -inverseMod32(p.value[p.offset+p.intLen-1]);
2144 
2145         for (int i=0, numWords = k >> 5; i < numWords; i++) {
2146             // V = R * c (mod 2^j)
2147             int  v = r * c.value[c.offset + c.intLen-1];
2148             // c = c + (v * p)
2149             p.mul(v, temp);
2150             c.add(temp);
2151             // c = c / 2^j
2152             c.intLen--;
2153         }
2154         int numBits = k & 0x1f;
2155         if (numBits != 0) {
2156             // V = R * c (mod 2^j)
2157             int v = r * c.value[c.offset + c.intLen-1];
2158             v &= ((1<<numBits) - 1);
2159             // c = c + (v * p)
2160             p.mul(v, temp);
2161             c.add(temp);
2162             // c = c / 2^j
2163             c.rightShift(numBits);
2164         }
2165