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
|