1 /*
2 * Copyright (c) 1999, 2007, Oracle and/or its affiliates. All rights reserved.
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * This code is free software; you can redistribute it and/or modify it
6 * under the terms of the GNU General Public License version 2 only, as
7 * published by the Free Software Foundation. Oracle designates this
8 * particular file as subject to the "Classpath" exception as provided
9 * by Oracle in the LICENSE file that accompanied this code.
10 *
11 * This code is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 * version 2 for more details (a copy is included in the LICENSE file that
15 * accompanied this code).
16 *
17 * You should have received a copy of the GNU General Public License version
18 * 2 along with this work; if not, write to the Free Software Foundation,
19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20 *
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22 * or visit www.oracle.com if you need additional information or have any
143 return 0;
144 long d = value[offset] & LONG_MASK;
145 return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d;
146 }
147
148 /**
149 * Convert this MutableBigInteger to a BigInteger object.
150 */
151 BigInteger toBigInteger(int sign) {
152 if (intLen == 0 || sign == 0)
153 return BigInteger.ZERO;
154 return new BigInteger(getMagnitudeArray(), sign);
155 }
156
157 /**
158 * Convert this MutableBigInteger to BigDecimal object with the specified sign
159 * and scale.
160 */
161 BigDecimal toBigDecimal(int sign, int scale) {
162 if (intLen == 0 || sign == 0)
163 return BigDecimal.valueOf(0, scale);
164 int[] mag = getMagnitudeArray();
165 int len = mag.length;
166 int d = mag[0];
167 // If this MutableBigInteger can't be fit into long, we need to
168 // make a BigInteger object for the resultant BigDecimal object.
169 if (len > 2 || (d < 0 && len == 2))
170 return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0);
171 long v = (len == 2) ?
172 ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
173 d & LONG_MASK;
174 return new BigDecimal(null, sign == -1 ? -v : v, scale, 0);
175 }
176
177 /**
178 * Clear out a MutableBigInteger for reuse.
179 */
180 void clear() {
181 offset = intLen = 0;
182 for (int index=0, n=value.length; index < n; index++)
183 value[index] = 0;
184 }
185
186 /**
187 * Set a MutableBigInteger to zero, removing its offset.
188 */
189 void reset() {
190 offset = intLen = 0;
191 }
192
193 /**
194 * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
527 * word input x, and subtracts the n word product from q. This is needed
528 * when subtracting qhat*divisor from dividend.
529 */
530 private int mulsub(int[] q, int[] a, int x, int len, int offset) {
531 long xLong = x & LONG_MASK;
532 long carry = 0;
533 offset += len;
534
535 for (int j=len-1; j >= 0; j--) {
536 long product = (a[j] & LONG_MASK) * xLong + carry;
537 long difference = q[offset] - product;
538 q[offset--] = (int)difference;
539 carry = (product >>> 32)
540 + (((difference & LONG_MASK) >
541 (((~(int)product) & LONG_MASK))) ? 1:0);
542 }
543 return (int)carry;
544 }
545
546 /**
547 * Right shift this MutableBigInteger n bits, where n is
548 * less than 32.
549 * Assumes that intLen > 0, n > 0 for speed
550 */
551 private final void primitiveRightShift(int n) {
552 int[] val = value;
553 int n2 = 32 - n;
554 for (int i=offset+intLen-1, c=val[i]; i>offset; i--) {
555 int b = c;
556 c = val[i-1];
557 val[i] = (c << n2) | (b >>> n);
558 }
559 val[offset] >>>= n;
560 }
561
562 /**
563 * Left shift this MutableBigInteger n bits, where n is
564 * less than 32.
565 * Assumes that intLen > 0, n > 0 for speed
566 */
825 return r;
826 }
827
828 if (quotient.value.length < intLen)
829 quotient.value = new int[intLen];
830 quotient.offset = 0;
831 quotient.intLen = intLen;
832
833 // Normalize the divisor
834 int shift = Integer.numberOfLeadingZeros(divisor);
835
836 int rem = value[offset];
837 long remLong = rem & LONG_MASK;
838 if (remLong < divisorLong) {
839 quotient.value[0] = 0;
840 } else {
841 quotient.value[0] = (int)(remLong / divisorLong);
842 rem = (int) (remLong - (quotient.value[0] * divisorLong));
843 remLong = rem & LONG_MASK;
844 }
845
846 int xlen = intLen;
847 int[] qWord = new int[2];
848 while (--xlen > 0) {
849 long dividendEstimate = (remLong<<32) |
850 (value[offset + intLen - xlen] & LONG_MASK);
851 if (dividendEstimate >= 0) {
852 qWord[0] = (int) (dividendEstimate / divisorLong);
853 qWord[1] = (int) (dividendEstimate - qWord[0] * divisorLong);
854 } else {
855 divWord(qWord, dividendEstimate, divisor);
856 }
857 quotient.value[intLen - xlen] = qWord[0];
858 rem = qWord[1];
859 remLong = rem & LONG_MASK;
860 }
861
862 quotient.normalize();
863 // Unnormalize
864 if (shift > 0)
865 return rem % divisor;
866 else
867 return rem;
868 }
869
870 /**
871 * Calculates the quotient of this div b and places the quotient in the
872 * provided MutableBigInteger objects and the remainder object is returned.
873 *
874 * Uses Algorithm D in Knuth section 4.3.1.
875 * Many optimizations to that algorithm have been adapted from the Colin
876 * Plumb C library.
877 * It special cases one word divisors for speed. The content of b is not
878 * changed.
879 *
880 */
881 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
882 if (b.intLen == 0)
883 throw new ArithmeticException("BigInteger divide by zero");
884
885 // Dividend is zero
886 if (intLen == 0) {
887 quotient.intLen = quotient.offset;
888 return new MutableBigInteger();
889 }
890
891 int cmp = compare(b);
892 // Dividend less than divisor
893 if (cmp < 0) {
894 quotient.intLen = quotient.offset = 0;
895 return new MutableBigInteger(this);
896 }
897 // Dividend equal to divisor
898 if (cmp == 0) {
899 quotient.value[0] = quotient.intLen = 1;
900 quotient.offset = 0;
901 return new MutableBigInteger();
902 }
903
904 quotient.clear();
905 // Special case one word divisor
906 if (b.intLen == 1) {
907 int r = divideOneWord(b.value[b.offset], quotient);
908 if (r == 0)
909 return new MutableBigInteger();
910 return new MutableBigInteger(r);
911 }
912
913 // Copy divisor value to protect divisor
914 int[] div = Arrays.copyOfRange(b.value, b.offset, b.offset + b.intLen);
915 return divideMagnitude(div, quotient);
916 }
917
918 /**
919 * Internally used to calculate the quotient of this div v and places the
920 * quotient in the provided MutableBigInteger object and the remainder is
921 * returned.
922 *
923 * @return the remainder of the division will be returned.
924 */
925 long divide(long v, MutableBigInteger quotient) {
926 if (v == 0)
927 throw new ArithmeticException("BigInteger divide by zero");
928
929 // Dividend is zero
930 if (intLen == 0) {
931 quotient.intLen = quotient.offset = 0;
932 return 0;
933 }
934 if (v < 0)
935 v = -v;
936
937 int d = (int)(v >>> 32);
938 quotient.clear();
939 // Special case on word divisor
940 if (d == 0)
941 return divideOneWord((int)v, quotient) & LONG_MASK;
942 else {
943 int[] div = new int[]{ d, (int)(v & LONG_MASK) };
944 return divideMagnitude(div, quotient).toLong();
945 }
946 }
947
948 /**
949 * Divide this MutableBigInteger by the divisor represented by its magnitude
950 * array. The quotient will be placed into the provided quotient object &
951 * the remainder object is returned.
952 */
953 private MutableBigInteger divideMagnitude(int[] divisor,
954 MutableBigInteger quotient) {
955
956 // Remainder starts as dividend with space for a leading zero
957 MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
958 System.arraycopy(value, offset, rem.value, 1, intLen);
959 rem.intLen = intLen;
960 rem.offset = 1;
961
962 int nlen = rem.intLen;
963
964 // Set the quotient size
965 int dlen = divisor.length;
966 int limit = nlen - dlen + 1;
967 if (quotient.value.length < limit) {
968 quotient.value = new int[limit];
969 quotient.offset = 0;
970 }
971 quotient.intLen = limit;
972 int[] q = quotient.value;
973
974 // D1 normalize the divisor
975 int shift = Integer.numberOfLeadingZeros(divisor[0]);
976 if (shift > 0) {
977 // First shift will not grow array
978 BigInteger.primitiveLeftShift(divisor, dlen, shift);
979 // But this one might
980 rem.leftShift(shift);
981 }
982
983 // Must insert leading 0 in rem if its length did not change
984 if (rem.intLen == nlen) {
985 rem.offset = 0;
986 rem.value[0] = 0;
987 rem.intLen++;
988 }
989
990 int dh = divisor[0];
991 long dhLong = dh & LONG_MASK;
992 int dl = divisor[1];
993 int[] qWord = new int[2];
994
995 // D2 Initialize j
996 for(int j=0; j<limit; j++) {
997 // D3 Calculate qhat
998 // estimate qhat
999 int qhat = 0;
1000 int qrem = 0;
1001 boolean skipCorrection = false;
1002 int nh = rem.value[j+rem.offset];
1003 int nh2 = nh + 0x80000000;
1004 int nm = rem.value[j+1+rem.offset];
1005
1006 if (nh == dh) {
1007 qhat = ~0;
1008 qrem = nh + nm;
1009 skipCorrection = qrem + 0x80000000 < nh2;
1010 } else {
1011 long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
1012 if (nChunk >= 0) {
1013 qhat = (int) (nChunk / dhLong);
1014 qrem = (int) (nChunk - (qhat * dhLong));
1015 } else {
1016 divWord(qWord, nChunk, dh);
1017 qhat = qWord[0];
1018 qrem = qWord[1];
1019 }
1020 }
1021
1022 if (qhat == 0)
1023 continue;
1024
1025 if (!skipCorrection) { // Correct qhat
1026 long nl = rem.value[j+2+rem.offset] & LONG_MASK;
1027 long rs = ((qrem & LONG_MASK) << 32) | nl;
1028 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1029
1030 if (unsignedLongCompare(estProduct, rs)) {
1031 qhat--;
1032 qrem = (int)((qrem & LONG_MASK) + dhLong);
1033 if ((qrem & LONG_MASK) >= dhLong) {
1034 estProduct -= (dl & LONG_MASK);
1035 rs = ((qrem & LONG_MASK) << 32) | nl;
1036 if (unsignedLongCompare(estProduct, rs))
1037 qhat--;
1038 }
1039 }
1040 }
1041
1042 // D4 Multiply and subtract
1043 rem.value[j+rem.offset] = 0;
1044 int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
1045
1046 // D5 Test remainder
1047 if (borrow + 0x80000000 > nh2) {
1048 // D6 Add back
1049 divadd(divisor, rem.value, j+1+rem.offset);
1050 qhat--;
1051 }
1052
1053 // Store the quotient digit
1054 q[j] = qhat;
1055 } // D7 loop on j
1056
1057 // D8 Unnormalize
1058 if (shift > 0)
1059 rem.rightShift(shift);
1060
1061 quotient.normalize();
1062 rem.normalize();
1063 return rem;
1064 }
1065
1066 /**
1067 * Compare two longs as if they were unsigned.
1068 * Returns true iff one is bigger than two.
1069 */
1070 private boolean unsignedLongCompare(long one, long two) {
1071 return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
1072 }
1073
1074 /**
1075 * This method divides a long quantity by an int to estimate
1076 * qhat for two multi precision numbers. It is used when
1077 * the signed value of n is less than zero.
1078 */
1079 private void divWord(int[] result, long n, int d) {
1080 long dLong = d & LONG_MASK;
1081
1082 if (dLong == 1) {
1083 result[0] = (int)n;
1084 result[1] = 0;
1085 return;
1086 }
1087
1088 // Approximate the quotient and remainder
1089 long q = (n >>> 1) / (dLong >>> 1);
1090 long r = n - q*dLong;
1091
1092 // Correct the approximation
1093 while (r < 0) {
1094 r += dLong;
1095 q--;
1096 }
1097 while (r >= dLong) {
1098 r -= dLong;
1099 q++;
1100 }
1101
1102 // n - q*dlong == r && 0 <= r <dLong, hence we're done.
1103 result[0] = (int)q;
1104 result[1] = (int)r;
1105 }
1106
1107 /**
1108 * Calculate GCD of this and b. This and b are changed by the computation.
1109 */
1110 MutableBigInteger hybridGCD(MutableBigInteger b) {
1111 // Use Euclid's algorithm until the numbers are approximately the
1112 // same length, then use the binary GCD algorithm to find the GCD.
1113 MutableBigInteger a = this;
1114 MutableBigInteger q = new MutableBigInteger();
1115
1116 while (b.intLen != 0) {
1117 if (Math.abs(a.intLen - b.intLen) < 2)
1118 return a.binaryGCD(b);
1119
1120 MutableBigInteger r = a.divide(b, q);
1121 a = b;
1122 b = r;
1123 }
1124 return a;
1456
1457 r = b.divide(a, q);
1458
1459 if (r.intLen == 0)
1460 throw new ArithmeticException("BigInteger not invertible.");
1461
1462 swapper = b;
1463 b = r;
1464
1465 if (q.intLen == 1)
1466 t0.mul(q.value[q.offset], temp);
1467 else
1468 q.multiply(t0, temp);
1469 swapper = q; q = temp; temp = swapper;
1470
1471 t1.add(q);
1472 }
1473 mod.subtract(t1);
1474 return mod;
1475 }
1476
1477 }
|
1 /*
2 * Copyright (c) 1999, 2011, Oracle and/or its affiliates. All rights reserved.
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * This code is free software; you can redistribute it and/or modify it
6 * under the terms of the GNU General Public License version 2 only, as
7 * published by the Free Software Foundation. Oracle designates this
8 * particular file as subject to the "Classpath" exception as provided
9 * by Oracle in the LICENSE file that accompanied this code.
10 *
11 * This code is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 * version 2 for more details (a copy is included in the LICENSE file that
15 * accompanied this code).
16 *
17 * You should have received a copy of the GNU General Public License version
18 * 2 along with this work; if not, write to the Free Software Foundation,
19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20 *
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22 * or visit www.oracle.com if you need additional information or have any
143 return 0;
144 long d = value[offset] & LONG_MASK;
145 return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d;
146 }
147
148 /**
149 * Convert this MutableBigInteger to a BigInteger object.
150 */
151 BigInteger toBigInteger(int sign) {
152 if (intLen == 0 || sign == 0)
153 return BigInteger.ZERO;
154 return new BigInteger(getMagnitudeArray(), sign);
155 }
156
157 /**
158 * Convert this MutableBigInteger to BigDecimal object with the specified sign
159 * and scale.
160 */
161 BigDecimal toBigDecimal(int sign, int scale) {
162 if (intLen == 0 || sign == 0)
163 return BigDecimal.zeroValueOf(scale);
164 int[] mag = getMagnitudeArray();
165 int len = mag.length;
166 int d = mag[0];
167 // If this MutableBigInteger can't be fit into long, we need to
168 // make a BigInteger object for the resultant BigDecimal object.
169 if (len > 2 || (d < 0 && len == 2))
170 return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0);
171 long v = (len == 2) ?
172 ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
173 d & LONG_MASK;
174 return BigDecimal.valueOf(sign == -1 ? -v : v, scale);
175 }
176
177 /**
178 * This is for internal use in converting from a MutableBigInteger
179 * object into a long value given a specified sign.
180 * returns INFLATED if value is not fit into long
181 */
182 long toCompactValue(int sign) {
183 if (intLen == 0 || sign == 0)
184 return 0L;
185 int[] mag = getMagnitudeArray();
186 int len = mag.length;
187 int d = mag[0];
188 // If this MutableBigInteger can not be fitted into long, we need to
189 // make a BigInteger object for the resultant BigDecimal object.
190 if (len > 2 || (d < 0 && len == 2))
191 return INFLATED;
192 long v = (len == 2) ?
193 ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
194 d & LONG_MASK;
195 return sign == -1 ? -v : v;
196 }
197
198 /**
199 * Clear out a MutableBigInteger for reuse.
200 */
201 void clear() {
202 offset = intLen = 0;
203 for (int index=0, n=value.length; index < n; index++)
204 value[index] = 0;
205 }
206
207 /**
208 * Set a MutableBigInteger to zero, removing its offset.
209 */
210 void reset() {
211 offset = intLen = 0;
212 }
213
214 /**
215 * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
548 * word input x, and subtracts the n word product from q. This is needed
549 * when subtracting qhat*divisor from dividend.
550 */
551 private int mulsub(int[] q, int[] a, int x, int len, int offset) {
552 long xLong = x & LONG_MASK;
553 long carry = 0;
554 offset += len;
555
556 for (int j=len-1; j >= 0; j--) {
557 long product = (a[j] & LONG_MASK) * xLong + carry;
558 long difference = q[offset] - product;
559 q[offset--] = (int)difference;
560 carry = (product >>> 32)
561 + (((difference & LONG_MASK) >
562 (((~(int)product) & LONG_MASK))) ? 1:0);
563 }
564 return (int)carry;
565 }
566
567 /**
568 * The method is the same as mulsun, except the fact that q array is not
569 * updated, the only result of the method is borrow flag.
570 */
571 private int mulsubBorrow(int[] q, int[] a, int x, int len, int offset) {
572 long xLong = x & LONG_MASK;
573 long carry = 0;
574 offset += len;
575 for (int j=len-1; j >= 0; j--) {
576 long product = (a[j] & LONG_MASK) * xLong + carry;
577 long difference = q[offset--] - product;
578 carry = (product >>> 32)
579 + (((difference & LONG_MASK) >
580 (((~(int)product) & LONG_MASK))) ? 1:0);
581 }
582 return (int)carry;
583 }
584
585 /**
586 * Right shift this MutableBigInteger n bits, where n is
587 * less than 32.
588 * Assumes that intLen > 0, n > 0 for speed
589 */
590 private final void primitiveRightShift(int n) {
591 int[] val = value;
592 int n2 = 32 - n;
593 for (int i=offset+intLen-1, c=val[i]; i>offset; i--) {
594 int b = c;
595 c = val[i-1];
596 val[i] = (c << n2) | (b >>> n);
597 }
598 val[offset] >>>= n;
599 }
600
601 /**
602 * Left shift this MutableBigInteger n bits, where n is
603 * less than 32.
604 * Assumes that intLen > 0, n > 0 for speed
605 */
864 return r;
865 }
866
867 if (quotient.value.length < intLen)
868 quotient.value = new int[intLen];
869 quotient.offset = 0;
870 quotient.intLen = intLen;
871
872 // Normalize the divisor
873 int shift = Integer.numberOfLeadingZeros(divisor);
874
875 int rem = value[offset];
876 long remLong = rem & LONG_MASK;
877 if (remLong < divisorLong) {
878 quotient.value[0] = 0;
879 } else {
880 quotient.value[0] = (int)(remLong / divisorLong);
881 rem = (int) (remLong - (quotient.value[0] * divisorLong));
882 remLong = rem & LONG_MASK;
883 }
884 int xlen = intLen;
885 while (--xlen > 0) {
886 long dividendEstimate = (remLong << 32) |
887 (value[offset + intLen - xlen] & LONG_MASK);
888 int q;
889 if (dividendEstimate >= 0) {
890 q = (int) (dividendEstimate / divisorLong);
891 rem = (int) (dividendEstimate - q * divisorLong);
892 } else {
893 long tmp = divWord(dividendEstimate, divisor);
894 q = (int) (tmp & LONG_MASK);
895 rem = (int) (tmp >>> 32);
896 }
897 quotient.value[intLen - xlen] = q;
898 remLong = rem & LONG_MASK;
899 }
900
901 quotient.normalize();
902 // Unnormalize
903 if (shift > 0)
904 return rem % divisor;
905 else
906 return rem;
907 }
908
909 /**
910 * Calculates the quotient of this div b and places the quotient in the
911 * provided MutableBigInteger objects and the remainder object is returned.
912 *
913 * Uses Algorithm D in Knuth section 4.3.1.
914 * Many optimizations to that algorithm have been adapted from the Colin
915 * Plumb C library.
916 * It special cases one word divisors for speed. The content of b is not
917 * changed.
918 *
919 */
920 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
921 return divide(b,quotient,true);
922 }
923
924 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needReminder) {
925 if (b.intLen == 0)
926 throw new ArithmeticException("BigInteger divide by zero");
927
928 // Dividend is zero
929 if (intLen == 0) {
930 quotient.intLen = quotient.offset;
931 return needReminder ? new MutableBigInteger() : null;
932 }
933
934 int cmp = compare(b);
935 // Dividend less than divisor
936 if (cmp < 0) {
937 quotient.intLen = quotient.offset = 0;
938 return needReminder ? new MutableBigInteger(this) : null;
939 }
940 // Dividend equal to divisor
941 if (cmp == 0) {
942 quotient.value[0] = quotient.intLen = 1;
943 quotient.offset = 0;
944 return needReminder ? new MutableBigInteger() : null;
945 }
946
947 quotient.clear();
948 // Special case one word divisor
949 if (b.intLen == 1) {
950 int r = divideOneWord(b.value[b.offset], quotient);
951 if(needReminder) {
952 if (r == 0)
953 return new MutableBigInteger();
954 return new MutableBigInteger(r);
955 } else {
956 return null;
957 }
958 }
959 return divideMagnitude(b, quotient, needReminder);
960 }
961
962 /**
963 * Internally used to calculate the quotient of this div v and places the
964 * quotient in the provided MutableBigInteger object and the remainder is
965 * returned.
966 *
967 * @return the remainder of the division will be returned.
968 */
969 long divide(long v, MutableBigInteger quotient) {
970 if (v == 0)
971 throw new ArithmeticException("BigInteger divide by zero");
972
973 // Dividend is zero
974 if (intLen == 0) {
975 quotient.intLen = quotient.offset = 0;
976 return 0;
977 }
978 if (v < 0)
979 v = -v;
980
981 int d = (int)(v >>> 32);
982 quotient.clear();
983 // Special case on word divisor
984 if (d == 0)
985 return divideOneWord((int)v, quotient) & LONG_MASK;
986 else {
987 return divideLongMagnitude(v, quotient).toLong();
988 }
989 }
990
991 private static void copyAndShift(int[] src, int srcFrom, int srcLen, int[] dst, int dstFrom, int shift) {
992 int n2 = 32 - shift;
993 int c=src[srcFrom];
994 for (int i=0; i < srcLen-1; i++) {
995 int b = c;
996 c = src[++srcFrom];
997 dst[dstFrom+i] = (b << shift) | (c >>> n2);
998 }
999 dst[dstFrom+srcLen-1] = c << shift;
1000 }
1001
1002 /**
1003 * Divide this MutableBigInteger by the divisor.
1004 * The quotient will be placed into the provided quotient object &
1005 * the remainder object is returned.
1006 */
1007 private MutableBigInteger divideMagnitude(MutableBigInteger div,
1008 MutableBigInteger quotient,
1009 boolean needReminder ) {
1010 // assert div.intLen > 1
1011 // D1 normalize the divisor
1012 int shift = Integer.numberOfLeadingZeros(div.value[div.offset]);
1013 // Copy divisor value to protect divisor
1014 final int dlen = div.intLen;
1015 int[] divisor;
1016 MutableBigInteger rem; // Remainder starts as dividend with space for a leading zero
1017 if (shift > 0) {
1018 divisor = new int[dlen];
1019 copyAndShift(div.value,div.offset,dlen,divisor,0,shift);
1020 if(Integer.numberOfLeadingZeros(value[offset])>=shift) {
1021 int[] remarr = new int[intLen + 1];
1022 rem = new MutableBigInteger(remarr);
1023 rem.intLen = intLen;
1024 rem.offset = 1;
1025 copyAndShift(value,offset,intLen,remarr,1,shift);
1026 } else {
1027 int[] remarr = new int[intLen + 2];
1028 rem = new MutableBigInteger(remarr);
1029 rem.intLen = intLen+1;
1030 rem.offset = 1;
1031 int rFrom = offset;
1032 int c=0;
1033 int n2 = 32 - shift;
1034 for (int i=1; i < intLen+1; i++,rFrom++) {
1035 int b = c;
1036 c = value[rFrom];
1037 remarr[i] = (b << shift) | (c >>> n2);
1038 }
1039 remarr[intLen+1] = c << shift;
1040 }
1041 } else {
1042 divisor = Arrays.copyOfRange(div.value, div.offset, div.offset + div.intLen);
1043 rem = new MutableBigInteger(new int[intLen + 1]);
1044 System.arraycopy(value, offset, rem.value, 1, intLen);
1045 rem.intLen = intLen;
1046 rem.offset = 1;
1047 }
1048
1049 int nlen = rem.intLen;
1050
1051 // Set the quotient size
1052 final int limit = nlen - dlen + 1;
1053 if (quotient.value.length < limit) {
1054 quotient.value = new int[limit];
1055 quotient.offset = 0;
1056 }
1057 quotient.intLen = limit;
1058 int[] q = quotient.value;
1059
1060
1061 // Must insert leading 0 in rem if its length did not change
1062 if (rem.intLen == nlen) {
1063 rem.offset = 0;
1064 rem.value[0] = 0;
1065 rem.intLen++;
1066 }
1067
1068 int dh = divisor[0];
1069 long dhLong = dh & LONG_MASK;
1070 int dl = divisor[1];
1071
1072 // D2 Initialize j
1073 for(int j=0; j<limit-1; j++) {
1074 // D3 Calculate qhat
1075 // estimate qhat
1076 int qhat = 0;
1077 int qrem = 0;
1078 boolean skipCorrection = false;
1079 int nh = rem.value[j+rem.offset];
1080 int nh2 = nh + 0x80000000;
1081 int nm = rem.value[j+1+rem.offset];
1082
1083 if (nh == dh) {
1084 qhat = ~0;
1085 qrem = nh + nm;
1086 skipCorrection = qrem + 0x80000000 < nh2;
1087 } else {
1088 long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
1089 if (nChunk >= 0) {
1090 qhat = (int) (nChunk / dhLong);
1091 qrem = (int) (nChunk - (qhat * dhLong));
1092 } else {
1093 long tmp = divWord(nChunk, dh);
1094 qhat = (int) (tmp & LONG_MASK);
1095 qrem = (int) (tmp >>> 32);
1096 }
1097 }
1098
1099 if (qhat == 0)
1100 continue;
1101
1102 if (!skipCorrection) { // Correct qhat
1103 long nl = rem.value[j+2+rem.offset] & LONG_MASK;
1104 long rs = ((qrem & LONG_MASK) << 32) | nl;
1105 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1106
1107 if (unsignedLongCompare(estProduct, rs)) {
1108 qhat--;
1109 qrem = (int)((qrem & LONG_MASK) + dhLong);
1110 if ((qrem & LONG_MASK) >= dhLong) {
1111 estProduct -= (dl & LONG_MASK);
1112 rs = ((qrem & LONG_MASK) << 32) | nl;
1113 if (unsignedLongCompare(estProduct, rs))
1114 qhat--;
1115 }
1116 }
1117 }
1118
1119 // D4 Multiply and subtract
1120 rem.value[j+rem.offset] = 0;
1121 int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
1122
1123 // D5 Test remainder
1124 if (borrow + 0x80000000 > nh2) {
1125 // D6 Add back
1126 divadd(divisor, rem.value, j+1+rem.offset);
1127 qhat--;
1128 }
1129
1130 // Store the quotient digit
1131 q[j] = qhat;
1132 } // D7 loop on j
1133 // D3 Calculate qhat
1134 // estimate qhat
1135 int qhat = 0;
1136 int qrem = 0;
1137 boolean skipCorrection = false;
1138 int nh = rem.value[limit - 1 + rem.offset];
1139 int nh2 = nh + 0x80000000;
1140 int nm = rem.value[limit + rem.offset];
1141
1142 if (nh == dh) {
1143 qhat = ~0;
1144 qrem = nh + nm;
1145 skipCorrection = qrem + 0x80000000 < nh2;
1146 } else {
1147 long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
1148 if (nChunk >= 0) {
1149 qhat = (int) (nChunk / dhLong);
1150 qrem = (int) (nChunk - (qhat * dhLong));
1151 } else {
1152 long tmp = divWord(nChunk, dh);
1153 qhat = (int) (tmp & LONG_MASK);
1154 qrem = (int) (tmp >>> 32);
1155 }
1156 }
1157 if (qhat != 0) {
1158 if (!skipCorrection) { // Correct qhat
1159 long nl = rem.value[limit + 1 + rem.offset] & LONG_MASK;
1160 long rs = ((qrem & LONG_MASK) << 32) | nl;
1161 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1162
1163 if (unsignedLongCompare(estProduct, rs)) {
1164 qhat--;
1165 qrem = (int) ((qrem & LONG_MASK) + dhLong);
1166 if ((qrem & LONG_MASK) >= dhLong) {
1167 estProduct -= (dl & LONG_MASK);
1168 rs = ((qrem & LONG_MASK) << 32) | nl;
1169 if (unsignedLongCompare(estProduct, rs))
1170 qhat--;
1171 }
1172 }
1173 }
1174
1175
1176 // D4 Multiply and subtract
1177 int borrow;
1178 rem.value[limit - 1 + rem.offset] = 0;
1179 if(needReminder)
1180 borrow = mulsub(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
1181 else
1182 borrow = mulsubBorrow(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
1183
1184 // D5 Test remainder
1185 if (borrow + 0x80000000 > nh2) {
1186 // D6 Add back
1187 if(needReminder)
1188 divadd(divisor, rem.value, limit - 1 + 1 + rem.offset);
1189 qhat--;
1190 }
1191
1192 // Store the quotient digit
1193 q[(limit - 1)] = qhat;
1194 }
1195
1196
1197 if(needReminder) {
1198 // D8 Unnormalize
1199 if (shift > 0)
1200 rem.rightShift(shift);
1201 rem.normalize();
1202 }
1203 quotient.normalize();
1204 return needReminder ? rem : null;
1205 }
1206
1207 /**
1208 * Divide this MutableBigInteger by the divisor represented by positive long
1209 * value. The quotient will be placed into the provided quotient object &
1210 * the remainder object is returned.
1211 */
1212 private MutableBigInteger divideLongMagnitude(long ldivisor, MutableBigInteger quotient) {
1213 // Remainder starts as dividend with space for a leading zero
1214 MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
1215 System.arraycopy(value, offset, rem.value, 1, intLen);
1216 rem.intLen = intLen;
1217 rem.offset = 1;
1218
1219 int nlen = rem.intLen;
1220
1221 int limit = nlen - 2 + 1;
1222 if (quotient.value.length < limit) {
1223 quotient.value = new int[limit];
1224 quotient.offset = 0;
1225 }
1226 quotient.intLen = limit;
1227 int[] q = quotient.value;
1228
1229 // D1 normalize the divisor
1230 int shift = Long.numberOfLeadingZeros(ldivisor);
1231 if (shift > 0) {
1232 ldivisor<<=shift;
1233 rem.leftShift(shift);
1234 }
1235
1236 // Must insert leading 0 in rem if its length did not change
1237 if (rem.intLen == nlen) {
1238 rem.offset = 0;
1239 rem.value[0] = 0;
1240 rem.intLen++;
1241 }
1242
1243 int dh = (int)(ldivisor >>> 32);
1244 long dhLong = dh & LONG_MASK;
1245 int dl = (int)(ldivisor & LONG_MASK);
1246
1247 // D2 Initialize j
1248 for (int j = 0; j < limit; j++) {
1249 // D3 Calculate qhat
1250 // estimate qhat
1251 int qhat = 0;
1252 int qrem = 0;
1253 boolean skipCorrection = false;
1254 int nh = rem.value[j + rem.offset];
1255 int nh2 = nh + 0x80000000;
1256 int nm = rem.value[j + 1 + rem.offset];
1257
1258 if (nh == dh) {
1259 qhat = ~0;
1260 qrem = nh + nm;
1261 skipCorrection = qrem + 0x80000000 < nh2;
1262 } else {
1263 long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
1264 if (nChunk >= 0) {
1265 qhat = (int) (nChunk / dhLong);
1266 qrem = (int) (nChunk - (qhat * dhLong));
1267 } else {
1268 long tmp = divWord(nChunk, dh);
1269 qhat =(int)(tmp & LONG_MASK);
1270 qrem = (int)(tmp>>>32);
1271 }
1272 }
1273
1274 if (qhat == 0)
1275 continue;
1276
1277 if (!skipCorrection) { // Correct qhat
1278 long nl = rem.value[j + 2 + rem.offset] & LONG_MASK;
1279 long rs = ((qrem & LONG_MASK) << 32) | nl;
1280 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1281
1282 if (unsignedLongCompare(estProduct, rs)) {
1283 qhat--;
1284 qrem = (int) ((qrem & LONG_MASK) + dhLong);
1285 if ((qrem & LONG_MASK) >= dhLong) {
1286 estProduct -= (dl & LONG_MASK);
1287 rs = ((qrem & LONG_MASK) << 32) | nl;
1288 if (unsignedLongCompare(estProduct, rs))
1289 qhat--;
1290 }
1291 }
1292 }
1293
1294 // D4 Multiply and subtract
1295 rem.value[j + rem.offset] = 0;
1296 int borrow = mulsubLong(rem.value, dh, dl, qhat, j + rem.offset);
1297
1298 // D5 Test remainder
1299 if (borrow + 0x80000000 > nh2) {
1300 // D6 Add back
1301 divaddLong(dh,dl, rem.value, j + 1 + rem.offset);
1302 qhat--;
1303 }
1304
1305 // Store the quotient digit
1306 q[j] = qhat;
1307 } // D7 loop on j
1308
1309 // D8 Unnormalize
1310 if (shift > 0)
1311 rem.rightShift(shift);
1312
1313 quotient.normalize();
1314 rem.normalize();
1315 return rem;
1316 }
1317
1318 /**
1319 * A primitive used for division by long.
1320 * Specialized version of the method divadd.
1321 * dh is a high part of the divisor, dl is a low part
1322 */
1323 private int divaddLong(int dh, int dl, int[] result, int offset) {
1324 long carry = 0;
1325
1326 long sum = (dl & LONG_MASK) + (result[1+offset] & LONG_MASK);
1327 result[1+offset] = (int)sum;
1328
1329 sum = (dh & LONG_MASK) + (result[offset] & LONG_MASK) + carry;
1330 result[offset] = (int)sum;
1331 carry = sum >>> 32;
1332 return (int)carry;
1333 }
1334
1335 /**
1336 * This method is used for division by long.
1337 * Specialized version of the method sulsub.
1338 * dh is a high part of the divisor, dl is a low part
1339 */
1340 private int mulsubLong(int[] q, int dh, int dl, int x, int offset) {
1341 long xLong = x & LONG_MASK;
1342 offset += 2;
1343 long product = (dl & LONG_MASK) * xLong;
1344 long difference = q[offset] - product;
1345 q[offset--] = (int)difference;
1346 long carry = (product >>> 32)
1347 + (((difference & LONG_MASK) >
1348 (((~(int)product) & LONG_MASK))) ? 1:0);
1349 product = (dh & LONG_MASK) * xLong + carry;
1350 difference = q[offset] - product;
1351 q[offset--] = (int)difference;
1352 carry = (product >>> 32)
1353 + (((difference & LONG_MASK) >
1354 (((~(int)product) & LONG_MASK))) ? 1:0);
1355 return (int)carry;
1356 }
1357
1358 /**
1359 * Compare two longs as if they were unsigned.
1360 * Returns true iff one is bigger than two.
1361 */
1362 private boolean unsignedLongCompare(long one, long two) {
1363 return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
1364 }
1365
1366 /**
1367 * This method divides a long quantity by an int to estimate
1368 * qhat for two multi precision numbers. It is used when
1369 * the signed value of n is less than zero.
1370 * Returns long value where high 32 bits contain reminder value and
1371 * low 32 bits contain quotient value.
1372 */
1373 static long divWord(long n, int d) {
1374 long dLong = d & LONG_MASK;
1375 long r;
1376 long q;
1377 if (dLong == 1) {
1378 q = (int)n;
1379 r = 0;
1380 return (r << 32) | (q & LONG_MASK);
1381 }
1382
1383 // Approximate the quotient and remainder
1384 q = (n >>> 1) / (dLong >>> 1);
1385 r = n - q*dLong;
1386
1387 // Correct the approximation
1388 while (r < 0) {
1389 r += dLong;
1390 q--;
1391 }
1392 while (r >= dLong) {
1393 r -= dLong;
1394 q++;
1395 }
1396 // n - q*dlong == r && 0 <= r <dLong, hence we're done.
1397 return (r << 32) | (q & LONG_MASK);
1398 }
1399
1400 /**
1401 * Calculate GCD of this and b. This and b are changed by the computation.
1402 */
1403 MutableBigInteger hybridGCD(MutableBigInteger b) {
1404 // Use Euclid's algorithm until the numbers are approximately the
1405 // same length, then use the binary GCD algorithm to find the GCD.
1406 MutableBigInteger a = this;
1407 MutableBigInteger q = new MutableBigInteger();
1408
1409 while (b.intLen != 0) {
1410 if (Math.abs(a.intLen - b.intLen) < 2)
1411 return a.binaryGCD(b);
1412
1413 MutableBigInteger r = a.divide(b, q);
1414 a = b;
1415 b = r;
1416 }
1417 return a;
1749
1750 r = b.divide(a, q);
1751
1752 if (r.intLen == 0)
1753 throw new ArithmeticException("BigInteger not invertible.");
1754
1755 swapper = b;
1756 b = r;
1757
1758 if (q.intLen == 1)
1759 t0.mul(q.value[q.offset], temp);
1760 else
1761 q.multiply(t0, temp);
1762 swapper = q; q = temp; temp = swapper;
1763
1764 t1.add(q);
1765 }
1766 mod.subtract(t1);
1767 return mod;
1768 }
1769 }
|