2 * Copyright (c) 1999, 2007, Oracle and/or its affiliates. All rights reserved.
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
5 * This code is free software; you can redistribute it and/or modify it
6 * under the terms of the GNU General Public License version 2 only, as
7 * published by the Free Software Foundation. Oracle designates this
8 * particular file as subject to the "Classpath" exception as provided
9 * by Oracle in the LICENSE file that accompanied this code.
11 * This code is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 * version 2 for more details (a copy is included in the LICENSE file that
15 * accompanied this code).
17 * You should have received a copy of the GNU General Public License version
18 * 2 along with this work; if not, write to the Free Software Foundation,
19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22 * or visit www.oracle.com if you need additional information or have any
29 * A class used to represent multiprecision integers that makes efficient
30 * use of allocated space by allowing a number to occupy only part of
31 * an array so that the arrays do not have to be reallocated as often.
32 * When performing an operation with many iterations the array used to
33 * hold a number is only reallocated when necessary and does not have to
34 * be the same size as the number it represents. A mutable number allows
35 * calculations to occur on the same number without having to create
36 * a new number for every step of the calculation as occurs with
40 * @author Michael McCloskey
44 import java.util.Arrays;
46 import static java.math.BigInteger.LONG_MASK;
47 import static java.math.BigDecimal.INFLATED;
49 class MutableBigInteger {
51 * Holds the magnitude of this MutableBigInteger in big endian order.
52 * The magnitude may start at an offset into the value array, and it may
53 * end before the length of the value array.
58 * The number of ints of the value array that are currently used
59 * to hold the magnitude of this MutableBigInteger. The magnitude starts
60 * at an offset and offset + intLen may be less than value.length.
65 * The offset into the value array where the magnitude of this
66 * MutableBigInteger begins.
72 * MutableBigInteger with one element value array with the value 1. Used by
73 * BigDecimal divideAndRound to increment the quotient. Use this constant
74 * only when the method is not going to modify this object.
76 static final MutableBigInteger ONE = new MutableBigInteger(1);
81 * The default constructor. An empty MutableBigInteger is created with
82 * a one word capacity.
90 * Construct a new MutableBigInteger with a magnitude specified by
93 MutableBigInteger(int val) {
100 * Construct a new MutableBigInteger with the specified value array
101 * up to the length of the array supplied.
103 MutableBigInteger(int[] val) {
109 * Construct a new MutableBigInteger with a magnitude equal to the
110 * specified BigInteger.
112 MutableBigInteger(BigInteger b) {
113 intLen = b.mag.length;
114 value = Arrays.copyOf(b.mag, intLen);
118 * Construct a new MutableBigInteger with a magnitude equal to the
119 * specified MutableBigInteger.
121 MutableBigInteger(MutableBigInteger val) {
123 value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen);
127 * Internal helper method to return the magnitude array. The caller is not
128 * supposed to modify the returned array.
130 private int[] getMagnitudeArray() {
131 if (offset > 0 || value.length != intLen)
132 return Arrays.copyOfRange(value, offset, offset + intLen);
137 * Convert this MutableBigInteger to a long value. The caller has to make
138 * sure this MutableBigInteger can be fit into long.
140 private long toLong() {
141 assert (intLen <= 2) : "this MutableBigInteger exceeds the range of long";
144 long d = value[offset] & LONG_MASK;
145 return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d;
149 * Convert this MutableBigInteger to a BigInteger object.
151 BigInteger toBigInteger(int sign) {
152 if (intLen == 0 || sign == 0)
153 return BigInteger.ZERO;
154 return new BigInteger(getMagnitudeArray(), sign);
158 * Convert this MutableBigInteger to BigDecimal object with the specified sign
161 BigDecimal toBigDecimal(int sign, int scale) {
162 if (intLen == 0 || sign == 0)
163 return BigDecimal.valueOf(0, scale);
164 int[] mag = getMagnitudeArray();
165 int len = mag.length;
167 // If this MutableBigInteger can't be fit into long, we need to
168 // make a BigInteger object for the resultant BigDecimal object.
169 if (len > 2 || (d < 0 && len == 2))
170 return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0);
171 long v = (len == 2) ?
172 ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
174 return new BigDecimal(null, sign == -1 ? -v : v, scale, 0);
178 * Clear out a MutableBigInteger for reuse.
182 for (int index=0, n=value.length; index < n; index++)
187 * Set a MutableBigInteger to zero, removing its offset.
194 * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
195 * as this MutableBigInteger is numerically less than, equal to, or
196 * greater than <tt>b</tt>.
198 final int compare(MutableBigInteger b) {
205 // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
207 int[] bval = b.value;
208 for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) {
209 int b1 = value[i] + 0x80000000;
210 int b2 = bval[j] + 0x80000000;
220 * Compare this against half of a MutableBigInteger object (Needed for
222 * Assumes no leading unnecessary zeros, which holds for results
225 final int compareHalf(MutableBigInteger b) {
229 return blen <=0 ? 0 : -1;
234 int[] bval = b.value;
237 // Only 2 cases left:len == blen or len == blen - 1
238 if (len != blen) { // len == blen - 1
239 if (bval[bstart] == 1) {
245 // compare values with right-shifted values of b,
246 // carrying shifted-out bits across words
248 for (int i = offset, j = bstart; i < len + offset;) {
250 long hb = ((bv >>> 1) + carry) & LONG_MASK;
251 long v = val[i++] & LONG_MASK;
253 return v < hb ? -1 : 1;
254 carry = (bv & 1) << 31; // carray will be either 0x80000000 or 0
256 return carry == 0? 0 : -1;
260 * Return the index of the lowest set bit in this MutableBigInteger. If the
261 * magnitude of this MutableBigInteger is zero, -1 is returned.
263 private final int getLowestSetBit() {
267 for (j=intLen-1; (j>0) && (value[j+offset]==0); j--)
272 return ((intLen-1-j)<<5) + Integer.numberOfTrailingZeros(b);
276 * Return the int in use in this MutableBigInteger at the specified
277 * index. This method is not used because it is not inlined on all
280 private final int getInt(int index) {
281 return value[offset+index];
285 * Return a long which is equal to the unsigned value of the int in
286 * use in this MutableBigInteger at the specified index. This method is
287 * not used because it is not inlined on all platforms.
289 private final long getLong(int index) {
290 return value[offset+index] & LONG_MASK;
294 * Ensure that the MutableBigInteger is in normal form, specifically
295 * making sure that there are no leading zeros, and that if the
296 * magnitude is zero, then intLen is zero.
298 final void normalize() {
305 if (value[index] != 0)
308 int indexBound = index+intLen;
311 } while(index < indexBound && value[index]==0);
313 int numZeros = index - offset;
315 offset = (intLen==0 ? 0 : offset+numZeros);
319 * If this MutableBigInteger cannot hold len words, increase the size
320 * of the value array to len words.
322 private final void ensureCapacity(int len) {
323 if (value.length < len) {
324 value = new int[len];
331 * Convert this MutableBigInteger into an int array with no leading
332 * zeros, of a length that is equal to this MutableBigInteger's intLen.
335 int[] result = new int[intLen];
336 for(int i=0; i<intLen; i++)
337 result[i] = value[offset+i];
342 * Sets the int at index+offset in this MutableBigInteger to val.
343 * This does not get inlined on all platforms so it is not used
344 * as often as originally intended.
346 void setInt(int index, int val) {
347 value[offset + index] = val;
351 * Sets this MutableBigInteger's value array to the specified array.
352 * The intLen is set to the specified length.
354 void setValue(int[] val, int length) {
361 * Sets this MutableBigInteger's value array to a copy of the specified
362 * array. The intLen is set to the length of the new array.
364 void copyValue(MutableBigInteger src) {
365 int len = src.intLen;
366 if (value.length < len)
367 value = new int[len];
368 System.arraycopy(src.value, src.offset, value, 0, len);
374 * Sets this MutableBigInteger's value array to a copy of the specified
375 * array. The intLen is set to the length of the specified array.
377 void copyValue(int[] val) {
378 int len = val.length;
379 if (value.length < len)
380 value = new int[len];
381 System.arraycopy(val, 0, value, 0, len);
387 * Returns true iff this MutableBigInteger has a value of one.
390 return (intLen == 1) && (value[offset] == 1);
394 * Returns true iff this MutableBigInteger has a value of zero.
397 return (intLen == 0);
401 * Returns true iff this MutableBigInteger is even.
404 return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0);
408 * Returns true iff this MutableBigInteger is odd.
411 return isZero() ? false : ((value[offset + intLen - 1] & 1) == 1);
415 * Returns true iff this MutableBigInteger is in normal form. A
416 * MutableBigInteger is in normal form if it has no leading zeros
417 * after the offset, and intLen + offset <= value.length.
420 if (intLen + offset > value.length)
424 return (value[offset] != 0);
428 * Returns a String representation of this MutableBigInteger in radix 10.
430 public String toString() {
431 BigInteger b = toBigInteger(1);
436 * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
439 void rightShift(int n) {
443 int nBits = n & 0x1F;
444 this.intLen -= nInts;
447 int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
448 if (nBits >= bitsInHighWord) {
449 this.primitiveLeftShift(32 - nBits);
452 primitiveRightShift(nBits);
457 * Left shift this MutableBigInteger n bits.
459 void leftShift(int n) {
461 * If there is enough storage space in this MutableBigInteger already
462 * the available space will be used. Space to the right of the used
463 * ints in the value array is faster to utilize, so the extra space
464 * will be taken from the right if possible.
470 int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
472 // If shift can be done without moving words, do so
473 if (n <= (32-bitsInHighWord)) {
474 primitiveLeftShift(nBits);
478 int newLen = intLen + nInts +1;
479 if (nBits <= (32-bitsInHighWord))
481 if (value.length < newLen) {
482 // The array must grow
483 int[] result = new int[newLen];
484 for (int i=0; i<intLen; i++)
485 result[i] = value[offset+i];
486 setValue(result, newLen);
487 } else if (value.length - offset >= newLen) {
488 // Use space on right
489 for(int i=0; i<newLen - intLen; i++)
490 value[offset+intLen+i] = 0;
492 // Must use space on left
493 for (int i=0; i<intLen; i++)
494 value[i] = value[offset+i];
495 for (int i=intLen; i<newLen; i++)
502 if (nBits <= (32-bitsInHighWord))
503 primitiveLeftShift(nBits);
505 primitiveRightShift(32 -nBits);
509 * A primitive used for division. This method adds in one multiple of the
510 * divisor a back to the dividend result at a specified offset. It is used
511 * when qhat was estimated too large, and must be adjusted.
513 private int divadd(int[] a, int[] result, int offset) {
516 for (int j=a.length-1; j >= 0; j--) {
517 long sum = (a[j] & LONG_MASK) +
518 (result[j+offset] & LONG_MASK) + carry;
519 result[j+offset] = (int)sum;
526 * This method is used for division. It multiplies an n word input a by one
527 * word input x, and subtracts the n word product from q. This is needed
528 * when subtracting qhat*divisor from dividend.
530 private int mulsub(int[] q, int[] a, int x, int len, int offset) {
531 long xLong = x & LONG_MASK;
535 for (int j=len-1; j >= 0; j--) {
536 long product = (a[j] & LONG_MASK) * xLong + carry;
537 long difference = q[offset] - product;
538 q[offset--] = (int)difference;
539 carry = (product >>> 32)
540 + (((difference & LONG_MASK) >
541 (((~(int)product) & LONG_MASK))) ? 1:0);
547 * Right shift this MutableBigInteger n bits, where n is
549 * Assumes that intLen > 0, n > 0 for speed
551 private final void primitiveRightShift(int n) {
554 for (int i=offset+intLen-1, c=val[i]; i>offset; i--) {
557 val[i] = (c << n2) | (b >>> n);
563 * Left shift this MutableBigInteger n bits, where n is
565 * Assumes that intLen > 0, n > 0 for speed
567 private final void primitiveLeftShift(int n) {
570 for (int i=offset, c=val[i], m=i+intLen-1; i<m; i++) {
573 val[i] = (b << n) | (c >>> n2);
575 val[offset+intLen-1] <<= n;
579 * Adds the contents of two MutableBigInteger objects.The result
580 * is placed within this MutableBigInteger.
581 * The contents of the addend are not changed.
583 void add(MutableBigInteger addend) {
585 int y = addend.intLen;
586 int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
587 int[] result = (value.length < resultLen ? new int[resultLen] : value);
589 int rstart = result.length-1;
593 // Add common parts of both numbers
596 sum = (value[x+offset] & LONG_MASK) +
597 (addend.value[y+addend.offset] & LONG_MASK) + carry;
598 result[rstart--] = (int)sum;
602 // Add remainder of the longer number
605 if (carry == 0 && result == value && rstart == (x + offset))
607 sum = (value[x+offset] & LONG_MASK) + carry;
608 result[rstart--] = (int)sum;
613 sum = (addend.value[y+addend.offset] & LONG_MASK) + carry;
614 result[rstart--] = (int)sum;
618 if (carry > 0) { // Result must grow in length
620 if (result.length < resultLen) {
621 int temp[] = new int[resultLen];
622 // Result one word longer from carry-out; copy low-order
623 // bits into new result.
624 System.arraycopy(result, 0, temp, 1, result.length);
628 result[rstart--] = 1;
634 offset = result.length - resultLen;
639 * Subtracts the smaller of this and b from the larger and places the
640 * result into this MutableBigInteger.
642 int subtract(MutableBigInteger b) {
643 MutableBigInteger a = this;
645 int[] result = value;
646 int sign = a.compare(b);
653 MutableBigInteger tmp = a;
658 int resultLen = a.intLen;
659 if (result.length < resultLen)
660 result = new int[resultLen];
665 int rstart = result.length - 1;
667 // Subtract common parts of both numbers
671 diff = (a.value[x+a.offset] & LONG_MASK) -
672 (b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32));
673 result[rstart--] = (int)diff;
675 // Subtract remainder of longer number
678 diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32));
679 result[rstart--] = (int)diff;
684 offset = value.length - resultLen;
690 * Subtracts the smaller of a and b from the larger and places the result
691 * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
692 * operation was performed.
694 private int difference(MutableBigInteger b) {
695 MutableBigInteger a = this;
696 int sign = a.compare(b);
700 MutableBigInteger tmp = a;
709 // Subtract common parts of both numbers
712 diff = (a.value[a.offset+ x] & LONG_MASK) -
713 (b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32));
714 a.value[a.offset+x] = (int)diff;
716 // Subtract remainder of longer number
719 diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32));
720 a.value[a.offset+x] = (int)diff;
728 * Multiply the contents of two MutableBigInteger objects. The result is
729 * placed into MutableBigInteger z. The contents of y are not changed.
731 void multiply(MutableBigInteger y, MutableBigInteger z) {
734 int newLen = xLen + yLen;
736 // Put z into an appropriate state to receive product
737 if (z.value.length < newLen)
738 z.value = new int[newLen];
742 // The first iteration is hoisted out of the loop to avoid extra add
744 for (int j=yLen-1, k=yLen+xLen-1; j >= 0; j--, k--) {
745 long product = (y.value[j+y.offset] & LONG_MASK) *
746 (value[xLen-1+offset] & LONG_MASK) + carry;
747 z.value[k] = (int)product;
748 carry = product >>> 32;
750 z.value[xLen-1] = (int)carry;
752 // Perform the multiplication word by word
753 for (int i = xLen-2; i >= 0; i--) {
755 for (int j=yLen-1, k=yLen+i; j >= 0; j--, k--) {
756 long product = (y.value[j+y.offset] & LONG_MASK) *
757 (value[i+offset] & LONG_MASK) +
758 (z.value[k] & LONG_MASK) + carry;
759 z.value[k] = (int)product;
760 carry = product >>> 32;
762 z.value[i] = (int)carry;
765 // Remove leading zeros from product
770 * Multiply the contents of this MutableBigInteger by the word y. The
771 * result is placed into z.
773 void mul(int y, MutableBigInteger z) {
784 // Perform the multiplication word by word
785 long ylong = y & LONG_MASK;
786 int[] zval = (z.value.length<intLen+1 ? new int[intLen + 1]
789 for (int i = intLen-1; i >= 0; i--) {
790 long product = ylong * (value[i+offset] & LONG_MASK) + carry;
791 zval[i+1] = (int)product;
792 carry = product >>> 32;
800 z.intLen = intLen + 1;
801 zval[0] = (int)carry;
807 * This method is used for division of an n word dividend by a one word
808 * divisor. The quotient is placed into quotient. The one word divisor is
809 * specified by divisor.
811 * @return the remainder of the division is returned.
814 int divideOneWord(int divisor, MutableBigInteger quotient) {
815 long divisorLong = divisor & LONG_MASK;
817 // Special case of one word dividend
819 long dividendValue = value[offset] & LONG_MASK;
820 int q = (int) (dividendValue / divisorLong);
821 int r = (int) (dividendValue - q * divisorLong);
822 quotient.value[0] = q;
823 quotient.intLen = (q == 0) ? 0 : 1;
828 if (quotient.value.length < intLen)
829 quotient.value = new int[intLen];
831 quotient.intLen = intLen;
833 // Normalize the divisor
834 int shift = Integer.numberOfLeadingZeros(divisor);
836 int rem = value[offset];
837 long remLong = rem & LONG_MASK;
838 if (remLong < divisorLong) {
839 quotient.value[0] = 0;
841 quotient.value[0] = (int)(remLong / divisorLong);
842 rem = (int) (remLong - (quotient.value[0] * divisorLong));
843 remLong = rem & LONG_MASK;
847 int[] qWord = new int[2];
849 long dividendEstimate = (remLong<<32) |
850 (value[offset + intLen - xlen] & LONG_MASK);
851 if (dividendEstimate >= 0) {
852 qWord[0] = (int) (dividendEstimate / divisorLong);
853 qWord[1] = (int) (dividendEstimate - qWord[0] * divisorLong);
855 divWord(qWord, dividendEstimate, divisor);
857 quotient.value[intLen - xlen] = qWord[0];
859 remLong = rem & LONG_MASK;
862 quotient.normalize();
865 return rem % divisor;
871 * Calculates the quotient of this div b and places the quotient in the
872 * provided MutableBigInteger objects and the remainder object is returned.
874 * Uses Algorithm D in Knuth section 4.3.1.
875 * Many optimizations to that algorithm have been adapted from the Colin
877 * It special cases one word divisors for speed. The content of b is not
881 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
883 throw new ArithmeticException("BigInteger divide by zero");
887 quotient.intLen = quotient.offset;
888 return new MutableBigInteger();
891 int cmp = compare(b);
892 // Dividend less than divisor
894 quotient.intLen = quotient.offset = 0;
895 return new MutableBigInteger(this);
897 // Dividend equal to divisor
899 quotient.value[0] = quotient.intLen = 1;
901 return new MutableBigInteger();
905 // Special case one word divisor
907 int r = divideOneWord(b.value[b.offset], quotient);
909 return new MutableBigInteger();
910 return new MutableBigInteger(r);
913 // Copy divisor value to protect divisor
914 int[] div = Arrays.copyOfRange(b.value, b.offset, b.offset + b.intLen);
915 return divideMagnitude(div, quotient);
919 * Internally used to calculate the quotient of this div v and places the
920 * quotient in the provided MutableBigInteger object and the remainder is
923 * @return the remainder of the division will be returned.
925 long divide(long v, MutableBigInteger quotient) {
927 throw new ArithmeticException("BigInteger divide by zero");
931 quotient.intLen = quotient.offset = 0;
937 int d = (int)(v >>> 32);
939 // Special case on word divisor
941 return divideOneWord((int)v, quotient) & LONG_MASK;
943 int[] div = new int[]{ d, (int)(v & LONG_MASK) };
944 return divideMagnitude(div, quotient).toLong();
949 * Divide this MutableBigInteger by the divisor represented by its magnitude
950 * array. The quotient will be placed into the provided quotient object &
951 * the remainder object is returned.
953 private MutableBigInteger divideMagnitude(int[] divisor,
954 MutableBigInteger quotient) {
956 // Remainder starts as dividend with space for a leading zero
957 MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
958 System.arraycopy(value, offset, rem.value, 1, intLen);
962 int nlen = rem.intLen;
964 // Set the quotient size
965 int dlen = divisor.length;
966 int limit = nlen - dlen + 1;
967 if (quotient.value.length < limit) {
968 quotient.value = new int[limit];
971 quotient.intLen = limit;
972 int[] q = quotient.value;
974 // D1 normalize the divisor
975 int shift = Integer.numberOfLeadingZeros(divisor[0]);
977 // First shift will not grow array
978 BigInteger.primitiveLeftShift(divisor, dlen, shift);
979 // But this one might
980 rem.leftShift(shift);
983 // Must insert leading 0 in rem if its length did not change
984 if (rem.intLen == nlen) {
991 long dhLong = dh & LONG_MASK;
993 int[] qWord = new int[2];
996 for(int j=0; j<limit; j++) {
1001 boolean skipCorrection = false;
1002 int nh = rem.value[j+rem.offset];
1003 int nh2 = nh + 0x80000000;
1004 int nm = rem.value[j+1+rem.offset];
1009 skipCorrection = qrem + 0x80000000 < nh2;
1011 long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
1013 qhat = (int) (nChunk / dhLong);
1014 qrem = (int) (nChunk - (qhat * dhLong));
1016 divWord(qWord, nChunk, dh);
1025 if (!skipCorrection) { // Correct qhat
1026 long nl = rem.value[j+2+rem.offset] & LONG_MASK;
1027 long rs = ((qrem & LONG_MASK) << 32) | nl;
1028 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1030 if (unsignedLongCompare(estProduct, rs)) {
1032 qrem = (int)((qrem & LONG_MASK) + dhLong);
1033 if ((qrem & LONG_MASK) >= dhLong) {
1034 estProduct -= (dl & LONG_MASK);
1035 rs = ((qrem & LONG_MASK) << 32) | nl;
1036 if (unsignedLongCompare(estProduct, rs))
1042 // D4 Multiply and subtract
1043 rem.value[j+rem.offset] = 0;
1044 int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
1046 // D5 Test remainder
1047 if (borrow + 0x80000000 > nh2) {
1049 divadd(divisor, rem.value, j+1+rem.offset);
1053 // Store the quotient digit
1059 rem.rightShift(shift);
1061 quotient.normalize();
1067 * Compare two longs as if they were unsigned.
1068 * Returns true iff one is bigger than two.
1070 private boolean unsignedLongCompare(long one, long two) {
1071 return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
1075 * This method divides a long quantity by an int to estimate
1076 * qhat for two multi precision numbers. It is used when
1077 * the signed value of n is less than zero.
1079 private void divWord(int[] result, long n, int d) {
1080 long dLong = d & LONG_MASK;
1088 // Approximate the quotient and remainder
1089 long q = (n >>> 1) / (dLong >>> 1);
1090 long r = n - q*dLong;
1092 // Correct the approximation
1097 while (r >= dLong) {
1102 // n - q*dlong == r && 0 <= r <dLong, hence we're done.
1108 * Calculate GCD of this and b. This and b are changed by the computation.
1110 MutableBigInteger hybridGCD(MutableBigInteger b) {
1111 // Use Euclid's algorithm until the numbers are approximately the
1112 // same length, then use the binary GCD algorithm to find the GCD.
1113 MutableBigInteger a = this;
1114 MutableBigInteger q = new MutableBigInteger();
1116 while (b.intLen != 0) {
1117 if (Math.abs(a.intLen - b.intLen) < 2)
1118 return a.binaryGCD(b);
1120 MutableBigInteger r = a.divide(b, q);
1128 * Calculate GCD of this and v.
1129 * Assumes that this and v are not zero.
1131 private MutableBigInteger binaryGCD(MutableBigInteger v) {
1132 // Algorithm B from Knuth section 4.5.2
1133 MutableBigInteger u = this;
1134 MutableBigInteger r = new MutableBigInteger();
1137 int s1 = u.getLowestSetBit();
1138 int s2 = v.getLowestSetBit();
1139 int k = (s1 < s2) ? s1 : s2;
1146 boolean uOdd = (k==s1);
1147 MutableBigInteger t = uOdd ? v: u;
1148 int tsign = uOdd ? -1 : 1;
1151 while ((lb = t.getLowestSetBit()) >= 0) {
1160 // Special case one word numbers
1161 if (u.intLen < 2 && v.intLen < 2) {
1162 int x = u.value[u.offset];
1163 int y = v.value[v.offset];
1164 x = binaryGcd(x, y);
1174 if ((tsign = u.difference(v)) == 0)
1176 t = (tsign >= 0) ? u : v;
1185 * Calculate GCD of a and b interpreted as unsigned integers.
1187 static int binaryGcd(int a, int b) {
1193 // Right shift a & b till their last bits equal to 1.
1194 int aZeros = Integer.numberOfTrailingZeros(a);
1195 int bZeros = Integer.numberOfTrailingZeros(b);
1199 int t = (aZeros < bZeros ? aZeros : bZeros);
1202 if ((a+0x80000000) > (b+0x80000000)) { // a > b as unsigned
1204 a >>>= Integer.numberOfTrailingZeros(a);
1207 b >>>= Integer.numberOfTrailingZeros(b);
1214 * Returns the modInverse of this mod p. This and p are not affected by
1217 MutableBigInteger mutableModInverse(MutableBigInteger p) {
1218 // Modulus is odd, use Schroeppel's algorithm
1220 return modInverse(p);
1222 // Base and modulus are even, throw exception
1224 throw new ArithmeticException("BigInteger not invertible.");
1226 // Get even part of modulus expressed as a power of 2
1227 int powersOf2 = p.getLowestSetBit();
1229 // Construct odd part of modulus
1230 MutableBigInteger oddMod = new MutableBigInteger(p);
1231 oddMod.rightShift(powersOf2);
1234 return modInverseMP2(powersOf2);
1236 // Calculate 1/a mod oddMod
1237 MutableBigInteger oddPart = modInverse(oddMod);
1239 // Calculate 1/a mod evenMod
1240 MutableBigInteger evenPart = modInverseMP2(powersOf2);
1242 // Combine the results using Chinese Remainder Theorem
1243 MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
1244 MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);
1246 MutableBigInteger temp1 = new MutableBigInteger();
1247 MutableBigInteger temp2 = new MutableBigInteger();
1248 MutableBigInteger result = new MutableBigInteger();
1250 oddPart.leftShift(powersOf2);
1251 oddPart.multiply(y1, result);
1253 evenPart.multiply(oddMod, temp1);
1254 temp1.multiply(y2, temp2);
1257 return result.divide(p, temp1);
1261 * Calculate the multiplicative inverse of this mod 2^k.
1263 MutableBigInteger modInverseMP2(int k) {
1265 throw new ArithmeticException("Non-invertible. (GCD != 1)");
1268 return euclidModInverse(k);
1270 int t = inverseMod32(value[offset+intLen-1]);
1273 t = (k == 32 ? t : t & ((1 << k) - 1));
1274 return new MutableBigInteger(t);
1277 long pLong = (value[offset+intLen-1] & LONG_MASK);
1279 pLong |= ((long)value[offset+intLen-2] << 32);
1280 long tLong = t & LONG_MASK;
1281 tLong = tLong * (2 - pLong * tLong); // 1 more Newton iter step
1282 tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
1284 MutableBigInteger result = new MutableBigInteger(new int[2]);
1285 result.value[0] = (int)(tLong >>> 32);
1286 result.value[1] = (int)tLong;
1293 * Returns the multiplicative inverse of val mod 2^32. Assumes val is odd.
1295 static int inverseMod32(int val) {
1296 // Newton's iteration!
1306 * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
1308 static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
1309 // Copy the mod to protect original
1310 return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);
1314 * Calculate the multiplicative inverse of this mod mod, where mod is odd.
1315 * This and mod are not changed by the calculation.
1317 * This method implements an algorithm due to Richard Schroeppel, that uses
1318 * the same intermediate representation as Montgomery Reduction
1319 * ("Montgomery Form"). The algorithm is described in an unpublished
1320 * manuscript entitled "Fast Modular Reciprocals."
1322 private MutableBigInteger modInverse(MutableBigInteger mod) {
1323 MutableBigInteger p = new MutableBigInteger(mod);
1324 MutableBigInteger f = new MutableBigInteger(this);
1325 MutableBigInteger g = new MutableBigInteger(p);
1326 SignedMutableBigInteger c = new SignedMutableBigInteger(1);
1327 SignedMutableBigInteger d = new SignedMutableBigInteger();
1328 MutableBigInteger temp = null;
1329 SignedMutableBigInteger sTemp = null;
1332 // Right shift f k times until odd, left shift d k times
1334 int trailingZeros = f.getLowestSetBit();
1335 f.rightShift(trailingZeros);
1336 d.leftShift(trailingZeros);
1340 // The Almost Inverse Algorithm
1342 // If gcd(f, g) != 1, number is not invertible modulo mod
1344 throw new ArithmeticException("BigInteger not invertible.");
1346 // If f < g exchange f, g and c, d
1347 if (f.compare(g) < 0) {
1348 temp = f; f = g; g = temp;
1349 sTemp = d; d = c; c = sTemp;
1352 // If f == g (mod 4)
1353 if (((f.value[f.offset + f.intLen - 1] ^
1354 g.value[g.offset + g.intLen - 1]) & 3) == 0) {
1356 c.signedSubtract(d);
1357 } else { // If f != g (mod 4)
1362 // Right shift f k times until odd, left shift d k times
1363 int trailingZeros = f.getLowestSetBit();
1364 f.rightShift(trailingZeros);
1365 d.leftShift(trailingZeros);
1372 return fixup(c, p, k);
1376 * The Fixup Algorithm
1377 * Calculates X such that X = C * 2^(-k) (mod P)
1378 * Assumes C<P and P is odd.
1380 static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,
1382 MutableBigInteger temp = new MutableBigInteger();
1383 // Set r to the multiplicative inverse of p mod 2^32
1384 int r = -inverseMod32(p.value[p.offset+p.intLen-1]);
1386 for(int i=0, numWords = k >> 5; i<numWords; i++) {
1387 // V = R * c (mod 2^j)
1388 int v = r * c.value[c.offset + c.intLen-1];
1395 int numBits = k & 0x1f;
1397 // V = R * c (mod 2^j)
1398 int v = r * c.value[c.offset + c.intLen-1];
1399 v &= ((1<<numBits) - 1);
1404 c.rightShift(numBits);
1407 // In theory, c may be greater than p at this point (Very rare!)
1408 while (c.compare(p) >= 0)
1415 * Uses the extended Euclidean algorithm to compute the modInverse of base
1416 * mod a modulus that is a power of 2. The modulus is 2^k.
1418 MutableBigInteger euclidModInverse(int k) {
1419 MutableBigInteger b = new MutableBigInteger(1);
1421 MutableBigInteger mod = new MutableBigInteger(b);
1423 MutableBigInteger a = new MutableBigInteger(this);
1424 MutableBigInteger q = new MutableBigInteger();
1425 MutableBigInteger r = b.divide(a, q);
1427 MutableBigInteger swapper = b;
1432 MutableBigInteger t1 = new MutableBigInteger(q);
1433 MutableBigInteger t0 = new MutableBigInteger(1);
1434 MutableBigInteger temp = new MutableBigInteger();
1436 while (!b.isOne()) {
1440 throw new ArithmeticException("BigInteger not invertible.");
1446 t1.mul(q.value[q.offset], temp);
1448 q.multiply(t1, temp);
1460 throw new ArithmeticException("BigInteger not invertible.");
1466 t0.mul(q.value[q.offset], temp);
1468 q.multiply(t0, temp);
1469 swapper = q; q = temp; temp = swapper;