0001: /*
0002: * Copyright 1999-2007 Sun Microsystems, Inc. All Rights Reserved.
0003: * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
0004: *
0005: * This code is free software; you can redistribute it and/or modify it
0006: * under the terms of the GNU General Public License version 2 only, as
0007: * published by the Free Software Foundation. Sun designates this
0008: * particular file as subject to the "Classpath" exception as provided
0009: * by Sun in the LICENSE file that accompanied this code.
0010: *
0011: * This code is distributed in the hope that it will be useful, but WITHOUT
0012: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
0013: * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
0014: * version 2 for more details (a copy is included in the LICENSE file that
0015: * accompanied this code).
0016: *
0017: * You should have received a copy of the GNU General Public License version
0018: * 2 along with this work; if not, write to the Free Software Foundation,
0019: * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
0020: *
0021: * Please contact Sun Microsystems, Inc., 4150 Network Circle, Santa Clara,
0022: * CA 95054 USA or visit www.sun.com if you need additional information or
0023: * have any questions.
0024: */
0025:
0026: package java.math;
0027:
0028: /**
0029: * A class used to represent multiprecision integers that makes efficient
0030: * use of allocated space by allowing a number to occupy only part of
0031: * an array so that the arrays do not have to be reallocated as often.
0032: * When performing an operation with many iterations the array used to
0033: * hold a number is only reallocated when necessary and does not have to
0034: * be the same size as the number it represents. A mutable number allows
0035: * calculations to occur on the same number without having to create
0036: * a new number for every step of the calculation as occurs with
0037: * BigIntegers.
0038: *
0039: * @see BigInteger
0040: * @author Michael McCloskey
0041: * @since 1.3
0042: */
0043:
0044: class MutableBigInteger {
0045: /**
0046: * Holds the magnitude of this MutableBigInteger in big endian order.
0047: * The magnitude may start at an offset into the value array, and it may
0048: * end before the length of the value array.
0049: */
0050: int[] value;
0051:
0052: /**
0053: * The number of ints of the value array that are currently used
0054: * to hold the magnitude of this MutableBigInteger. The magnitude starts
0055: * at an offset and offset + intLen may be less than value.length.
0056: */
0057: int intLen;
0058:
0059: /**
0060: * The offset into the value array where the magnitude of this
0061: * MutableBigInteger begins.
0062: */
0063: int offset = 0;
0064:
0065: /**
0066: * This mask is used to obtain the value of an int as if it were unsigned.
0067: */
0068: private final static long LONG_MASK = 0xffffffffL;
0069:
0070: // Constructors
0071:
0072: /**
0073: * The default constructor. An empty MutableBigInteger is created with
0074: * a one word capacity.
0075: */
0076: MutableBigInteger() {
0077: value = new int[1];
0078: intLen = 0;
0079: }
0080:
0081: /**
0082: * Construct a new MutableBigInteger with a magnitude specified by
0083: * the int val.
0084: */
0085: MutableBigInteger(int val) {
0086: value = new int[1];
0087: intLen = 1;
0088: value[0] = val;
0089: }
0090:
0091: /**
0092: * Construct a new MutableBigInteger with the specified value array
0093: * up to the specified length.
0094: */
0095: MutableBigInteger(int[] val, int len) {
0096: value = val;
0097: intLen = len;
0098: }
0099:
0100: /**
0101: * Construct a new MutableBigInteger with the specified value array
0102: * up to the length of the array supplied.
0103: */
0104: MutableBigInteger(int[] val) {
0105: value = val;
0106: intLen = val.length;
0107: }
0108:
0109: /**
0110: * Construct a new MutableBigInteger with a magnitude equal to the
0111: * specified BigInteger.
0112: */
0113: MutableBigInteger(BigInteger b) {
0114: value = (int[]) b.mag.clone();
0115: intLen = value.length;
0116: }
0117:
0118: /**
0119: * Construct a new MutableBigInteger with a magnitude equal to the
0120: * specified MutableBigInteger.
0121: */
0122: MutableBigInteger(MutableBigInteger val) {
0123: intLen = val.intLen;
0124: value = new int[intLen];
0125:
0126: for (int i = 0; i < intLen; i++)
0127: value[i] = val.value[val.offset + i];
0128: }
0129:
0130: /**
0131: * Clear out a MutableBigInteger for reuse.
0132: */
0133: void clear() {
0134: offset = intLen = 0;
0135: for (int index = 0, n = value.length; index < n; index++)
0136: value[index] = 0;
0137: }
0138:
0139: /**
0140: * Set a MutableBigInteger to zero, removing its offset.
0141: */
0142: void reset() {
0143: offset = intLen = 0;
0144: }
0145:
0146: /**
0147: * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
0148: * as this MutableBigInteger is numerically less than, equal to, or
0149: * greater than {@code b}.
0150: */
0151: final int compare(MutableBigInteger b) {
0152: if (intLen < b.intLen)
0153: return -1;
0154: if (intLen > b.intLen)
0155: return 1;
0156:
0157: for (int i = 0; i < intLen; i++) {
0158: int b1 = value[offset + i] + 0x80000000;
0159: int b2 = b.value[b.offset + i] + 0x80000000;
0160: if (b1 < b2)
0161: return -1;
0162: if (b1 > b2)
0163: return 1;
0164: }
0165: return 0;
0166: }
0167:
0168: /**
0169: * Return the index of the lowest set bit in this MutableBigInteger. If the
0170: * magnitude of this MutableBigInteger is zero, -1 is returned.
0171: */
0172: private final int getLowestSetBit() {
0173: if (intLen == 0)
0174: return -1;
0175: int j, b;
0176: for (j = intLen - 1; (j > 0) && (value[j + offset] == 0); j--)
0177: ;
0178: b = value[j + offset];
0179: if (b == 0)
0180: return -1;
0181: return ((intLen - 1 - j) << 5) + BigInteger.trailingZeroCnt(b);
0182: }
0183:
0184: /**
0185: * Return the int in use in this MutableBigInteger at the specified
0186: * index. This method is not used because it is not inlined on all
0187: * platforms.
0188: */
0189: private final int getInt(int index) {
0190: return value[offset + index];
0191: }
0192:
0193: /**
0194: * Return a long which is equal to the unsigned value of the int in
0195: * use in this MutableBigInteger at the specified index. This method is
0196: * not used because it is not inlined on all platforms.
0197: */
0198: private final long getLong(int index) {
0199: return value[offset + index] & LONG_MASK;
0200: }
0201:
0202: /**
0203: * Ensure that the MutableBigInteger is in normal form, specifically
0204: * making sure that there are no leading zeros, and that if the
0205: * magnitude is zero, then intLen is zero.
0206: */
0207: final void normalize() {
0208: if (intLen == 0) {
0209: offset = 0;
0210: return;
0211: }
0212:
0213: int index = offset;
0214: if (value[index] != 0)
0215: return;
0216:
0217: int indexBound = index + intLen;
0218: do {
0219: index++;
0220: } while (index < indexBound && value[index] == 0);
0221:
0222: int numZeros = index - offset;
0223: intLen -= numZeros;
0224: offset = (intLen == 0 ? 0 : offset + numZeros);
0225: }
0226:
0227: /**
0228: * If this MutableBigInteger cannot hold len words, increase the size
0229: * of the value array to len words.
0230: */
0231: private final void ensureCapacity(int len) {
0232: if (value.length < len) {
0233: value = new int[len];
0234: offset = 0;
0235: intLen = len;
0236: }
0237: }
0238:
0239: /**
0240: * Convert this MutableBigInteger into an int array with no leading
0241: * zeros, of a length that is equal to this MutableBigInteger's intLen.
0242: */
0243: int[] toIntArray() {
0244: int[] result = new int[intLen];
0245: for (int i = 0; i < intLen; i++)
0246: result[i] = value[offset + i];
0247: return result;
0248: }
0249:
0250: /**
0251: * Sets the int at index+offset in this MutableBigInteger to val.
0252: * This does not get inlined on all platforms so it is not used
0253: * as often as originally intended.
0254: */
0255: void setInt(int index, int val) {
0256: value[offset + index] = val;
0257: }
0258:
0259: /**
0260: * Sets this MutableBigInteger's value array to the specified array.
0261: * The intLen is set to the specified length.
0262: */
0263: void setValue(int[] val, int length) {
0264: value = val;
0265: intLen = length;
0266: offset = 0;
0267: }
0268:
0269: /**
0270: * Sets this MutableBigInteger's value array to a copy of the specified
0271: * array. The intLen is set to the length of the new array.
0272: */
0273: void copyValue(MutableBigInteger val) {
0274: int len = val.intLen;
0275: if (value.length < len)
0276: value = new int[len];
0277:
0278: for (int i = 0; i < len; i++)
0279: value[i] = val.value[val.offset + i];
0280: intLen = len;
0281: offset = 0;
0282: }
0283:
0284: /**
0285: * Sets this MutableBigInteger's value array to a copy of the specified
0286: * array. The intLen is set to the length of the specified array.
0287: */
0288: void copyValue(int[] val) {
0289: int len = val.length;
0290: if (value.length < len)
0291: value = new int[len];
0292: for (int i = 0; i < len; i++)
0293: value[i] = val[i];
0294: intLen = len;
0295: offset = 0;
0296: }
0297:
0298: /**
0299: * Returns true iff this MutableBigInteger has a value of one.
0300: */
0301: boolean isOne() {
0302: return (intLen == 1) && (value[offset] == 1);
0303: }
0304:
0305: /**
0306: * Returns true iff this MutableBigInteger has a value of zero.
0307: */
0308: boolean isZero() {
0309: return (intLen == 0);
0310: }
0311:
0312: /**
0313: * Returns true iff this MutableBigInteger is even.
0314: */
0315: boolean isEven() {
0316: return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0);
0317: }
0318:
0319: /**
0320: * Returns true iff this MutableBigInteger is odd.
0321: */
0322: boolean isOdd() {
0323: return ((value[offset + intLen - 1] & 1) == 1);
0324: }
0325:
0326: /**
0327: * Returns true iff this MutableBigInteger is in normal form. A
0328: * MutableBigInteger is in normal form if it has no leading zeros
0329: * after the offset, and intLen + offset <= value.length.
0330: */
0331: boolean isNormal() {
0332: if (intLen + offset > value.length)
0333: return false;
0334: if (intLen == 0)
0335: return true;
0336: return (value[offset] != 0);
0337: }
0338:
0339: /**
0340: * Returns a String representation of this MutableBigInteger in radix 10.
0341: */
0342: public String toString() {
0343: BigInteger b = new BigInteger(this , 1);
0344: return b.toString();
0345: }
0346:
0347: /**
0348: * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
0349: * in normal form.
0350: */
0351: void rightShift(int n) {
0352: if (intLen == 0)
0353: return;
0354: int nInts = n >>> 5;
0355: int nBits = n & 0x1F;
0356: this .intLen -= nInts;
0357: if (nBits == 0)
0358: return;
0359: int bitsInHighWord = BigInteger.bitLen(value[offset]);
0360: if (nBits >= bitsInHighWord) {
0361: this .primitiveLeftShift(32 - nBits);
0362: this .intLen--;
0363: } else {
0364: primitiveRightShift(nBits);
0365: }
0366: }
0367:
0368: /**
0369: * Left shift this MutableBigInteger n bits.
0370: */
0371: void leftShift(int n) {
0372: /*
0373: * If there is enough storage space in this MutableBigInteger already
0374: * the available space will be used. Space to the right of the used
0375: * ints in the value array is faster to utilize, so the extra space
0376: * will be taken from the right if possible.
0377: */
0378: if (intLen == 0)
0379: return;
0380: int nInts = n >>> 5;
0381: int nBits = n & 0x1F;
0382: int bitsInHighWord = BigInteger.bitLen(value[offset]);
0383:
0384: // If shift can be done without moving words, do so
0385: if (n <= (32 - bitsInHighWord)) {
0386: primitiveLeftShift(nBits);
0387: return;
0388: }
0389:
0390: int newLen = intLen + nInts + 1;
0391: if (nBits <= (32 - bitsInHighWord))
0392: newLen--;
0393: if (value.length < newLen) {
0394: // The array must grow
0395: int[] result = new int[newLen];
0396: for (int i = 0; i < intLen; i++)
0397: result[i] = value[offset + i];
0398: setValue(result, newLen);
0399: } else if (value.length - offset >= newLen) {
0400: // Use space on right
0401: for (int i = 0; i < newLen - intLen; i++)
0402: value[offset + intLen + i] = 0;
0403: } else {
0404: // Must use space on left
0405: for (int i = 0; i < intLen; i++)
0406: value[i] = value[offset + i];
0407: for (int i = intLen; i < newLen; i++)
0408: value[i] = 0;
0409: offset = 0;
0410: }
0411: intLen = newLen;
0412: if (nBits == 0)
0413: return;
0414: if (nBits <= (32 - bitsInHighWord))
0415: primitiveLeftShift(nBits);
0416: else
0417: primitiveRightShift(32 - nBits);
0418: }
0419:
0420: /**
0421: * A primitive used for division. This method adds in one multiple of the
0422: * divisor a back to the dividend result at a specified offset. It is used
0423: * when qhat was estimated too large, and must be adjusted.
0424: */
0425: private int divadd(int[] a, int[] result, int offset) {
0426: long carry = 0;
0427:
0428: for (int j = a.length - 1; j >= 0; j--) {
0429: long sum = (a[j] & LONG_MASK)
0430: + (result[j + offset] & LONG_MASK) + carry;
0431: result[j + offset] = (int) sum;
0432: carry = sum >>> 32;
0433: }
0434: return (int) carry;
0435: }
0436:
0437: /**
0438: * This method is used for division. It multiplies an n word input a by one
0439: * word input x, and subtracts the n word product from q. This is needed
0440: * when subtracting qhat*divisor from dividend.
0441: */
0442: private int mulsub(int[] q, int[] a, int x, int len, int offset) {
0443: long xLong = x & LONG_MASK;
0444: long carry = 0;
0445: offset += len;
0446:
0447: for (int j = len - 1; j >= 0; j--) {
0448: long product = (a[j] & LONG_MASK) * xLong + carry;
0449: long difference = q[offset] - product;
0450: q[offset--] = (int) difference;
0451: carry = (product >>> 32)
0452: + (((difference & LONG_MASK) > (((~(int) product) & LONG_MASK))) ? 1
0453: : 0);
0454: }
0455: return (int) carry;
0456: }
0457:
0458: /**
0459: * Right shift this MutableBigInteger n bits, where n is
0460: * less than 32.
0461: * Assumes that intLen > 0, n > 0 for speed
0462: */
0463: private final void primitiveRightShift(int n) {
0464: int[] val = value;
0465: int n2 = 32 - n;
0466: for (int i = offset + intLen - 1, c = val[i]; i > offset; i--) {
0467: int b = c;
0468: c = val[i - 1];
0469: val[i] = (c << n2) | (b >>> n);
0470: }
0471: val[offset] >>>= n;
0472: }
0473:
0474: /**
0475: * Left shift this MutableBigInteger n bits, where n is
0476: * less than 32.
0477: * Assumes that intLen > 0, n > 0 for speed
0478: */
0479: private final void primitiveLeftShift(int n) {
0480: int[] val = value;
0481: int n2 = 32 - n;
0482: for (int i = offset, c = val[i], m = i + intLen - 1; i < m; i++) {
0483: int b = c;
0484: c = val[i + 1];
0485: val[i] = (b << n) | (c >>> n2);
0486: }
0487: val[offset + intLen - 1] <<= n;
0488: }
0489:
0490: /**
0491: * Adds the contents of two MutableBigInteger objects.The result
0492: * is placed within this MutableBigInteger.
0493: * The contents of the addend are not changed.
0494: */
0495: void add(MutableBigInteger addend) {
0496: int x = intLen;
0497: int y = addend.intLen;
0498: int resultLen = (intLen > addend.intLen ? intLen
0499: : addend.intLen);
0500: int[] result = (value.length < resultLen ? new int[resultLen]
0501: : value);
0502:
0503: int rstart = result.length - 1;
0504: long sum = 0;
0505:
0506: // Add common parts of both numbers
0507: while (x > 0 && y > 0) {
0508: x--;
0509: y--;
0510: sum = (value[x + offset] & LONG_MASK)
0511: + (addend.value[y + addend.offset] & LONG_MASK)
0512: + (sum >>> 32);
0513: result[rstart--] = (int) sum;
0514: }
0515:
0516: // Add remainder of the longer number
0517: while (x > 0) {
0518: x--;
0519: sum = (value[x + offset] & LONG_MASK) + (sum >>> 32);
0520: result[rstart--] = (int) sum;
0521: }
0522: while (y > 0) {
0523: y--;
0524: sum = (addend.value[y + addend.offset] & LONG_MASK)
0525: + (sum >>> 32);
0526: result[rstart--] = (int) sum;
0527: }
0528:
0529: if ((sum >>> 32) > 0) { // Result must grow in length
0530: resultLen++;
0531: if (result.length < resultLen) {
0532: int temp[] = new int[resultLen];
0533: for (int i = resultLen - 1; i > 0; i--)
0534: temp[i] = result[i - 1];
0535: temp[0] = 1;
0536: result = temp;
0537: } else {
0538: result[rstart--] = 1;
0539: }
0540: }
0541:
0542: value = result;
0543: intLen = resultLen;
0544: offset = result.length - resultLen;
0545: }
0546:
0547: /**
0548: * Subtracts the smaller of this and b from the larger and places the
0549: * result into this MutableBigInteger.
0550: */
0551: int subtract(MutableBigInteger b) {
0552: MutableBigInteger a = this ;
0553:
0554: int[] result = value;
0555: int sign = a.compare(b);
0556:
0557: if (sign == 0) {
0558: reset();
0559: return 0;
0560: }
0561: if (sign < 0) {
0562: MutableBigInteger tmp = a;
0563: a = b;
0564: b = tmp;
0565: }
0566:
0567: int resultLen = a.intLen;
0568: if (result.length < resultLen)
0569: result = new int[resultLen];
0570:
0571: long diff = 0;
0572: int x = a.intLen;
0573: int y = b.intLen;
0574: int rstart = result.length - 1;
0575:
0576: // Subtract common parts of both numbers
0577: while (y > 0) {
0578: x--;
0579: y--;
0580:
0581: diff = (a.value[x + a.offset] & LONG_MASK)
0582: - (b.value[y + b.offset] & LONG_MASK)
0583: - ((int) -(diff >> 32));
0584: result[rstart--] = (int) diff;
0585: }
0586: // Subtract remainder of longer number
0587: while (x > 0) {
0588: x--;
0589: diff = (a.value[x + a.offset] & LONG_MASK)
0590: - ((int) -(diff >> 32));
0591: result[rstart--] = (int) diff;
0592: }
0593:
0594: value = result;
0595: intLen = resultLen;
0596: offset = value.length - resultLen;
0597: normalize();
0598: return sign;
0599: }
0600:
0601: /**
0602: * Subtracts the smaller of a and b from the larger and places the result
0603: * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
0604: * operation was performed.
0605: */
0606: private int difference(MutableBigInteger b) {
0607: MutableBigInteger a = this ;
0608: int sign = a.compare(b);
0609: if (sign == 0)
0610: return 0;
0611: if (sign < 0) {
0612: MutableBigInteger tmp = a;
0613: a = b;
0614: b = tmp;
0615: }
0616:
0617: long diff = 0;
0618: int x = a.intLen;
0619: int y = b.intLen;
0620:
0621: // Subtract common parts of both numbers
0622: while (y > 0) {
0623: x--;
0624: y--;
0625: diff = (a.value[a.offset + x] & LONG_MASK)
0626: - (b.value[b.offset + y] & LONG_MASK)
0627: - ((int) -(diff >> 32));
0628: a.value[a.offset + x] = (int) diff;
0629: }
0630: // Subtract remainder of longer number
0631: while (x > 0) {
0632: x--;
0633: diff = (a.value[a.offset + x] & LONG_MASK)
0634: - ((int) -(diff >> 32));
0635: a.value[a.offset + x] = (int) diff;
0636: }
0637:
0638: a.normalize();
0639: return sign;
0640: }
0641:
0642: /**
0643: * Multiply the contents of two MutableBigInteger objects. The result is
0644: * placed into MutableBigInteger z. The contents of y are not changed.
0645: */
0646: void multiply(MutableBigInteger y, MutableBigInteger z) {
0647: int xLen = intLen;
0648: int yLen = y.intLen;
0649: int newLen = xLen + yLen;
0650:
0651: // Put z into an appropriate state to receive product
0652: if (z.value.length < newLen)
0653: z.value = new int[newLen];
0654: z.offset = 0;
0655: z.intLen = newLen;
0656:
0657: // The first iteration is hoisted out of the loop to avoid extra add
0658: long carry = 0;
0659: for (int j = yLen - 1, k = yLen + xLen - 1; j >= 0; j--, k--) {
0660: long product = (y.value[j + y.offset] & LONG_MASK)
0661: * (value[xLen - 1 + offset] & LONG_MASK) + carry;
0662: z.value[k] = (int) product;
0663: carry = product >>> 32;
0664: }
0665: z.value[xLen - 1] = (int) carry;
0666:
0667: // Perform the multiplication word by word
0668: for (int i = xLen - 2; i >= 0; i--) {
0669: carry = 0;
0670: for (int j = yLen - 1, k = yLen + i; j >= 0; j--, k--) {
0671: long product = (y.value[j + y.offset] & LONG_MASK)
0672: * (value[i + offset] & LONG_MASK)
0673: + (z.value[k] & LONG_MASK) + carry;
0674: z.value[k] = (int) product;
0675: carry = product >>> 32;
0676: }
0677: z.value[i] = (int) carry;
0678: }
0679:
0680: // Remove leading zeros from product
0681: z.normalize();
0682: }
0683:
0684: /**
0685: * Multiply the contents of this MutableBigInteger by the word y. The
0686: * result is placed into z.
0687: */
0688: void mul(int y, MutableBigInteger z) {
0689: if (y == 1) {
0690: z.copyValue(this );
0691: return;
0692: }
0693:
0694: if (y == 0) {
0695: z.clear();
0696: return;
0697: }
0698:
0699: // Perform the multiplication word by word
0700: long ylong = y & LONG_MASK;
0701: int[] zval = (z.value.length < intLen + 1 ? new int[intLen + 1]
0702: : z.value);
0703: long carry = 0;
0704: for (int i = intLen - 1; i >= 0; i--) {
0705: long product = ylong * (value[i + offset] & LONG_MASK)
0706: + carry;
0707: zval[i + 1] = (int) product;
0708: carry = product >>> 32;
0709: }
0710:
0711: if (carry == 0) {
0712: z.offset = 1;
0713: z.intLen = intLen;
0714: } else {
0715: z.offset = 0;
0716: z.intLen = intLen + 1;
0717: zval[0] = (int) carry;
0718: }
0719: z.value = zval;
0720: }
0721:
0722: /**
0723: * This method is used for division of an n word dividend by a one word
0724: * divisor. The quotient is placed into quotient. The one word divisor is
0725: * specified by divisor. The value of this MutableBigInteger is the
0726: * dividend at invocation but is replaced by the remainder.
0727: *
0728: * NOTE: The value of this MutableBigInteger is modified by this method.
0729: */
0730: void divideOneWord(int divisor, MutableBigInteger quotient) {
0731: long divLong = divisor & LONG_MASK;
0732:
0733: // Special case of one word dividend
0734: if (intLen == 1) {
0735: long remValue = value[offset] & LONG_MASK;
0736: quotient.value[0] = (int) (remValue / divLong);
0737: quotient.intLen = (quotient.value[0] == 0) ? 0 : 1;
0738: quotient.offset = 0;
0739:
0740: value[0] = (int) (remValue - (quotient.value[0] * divLong));
0741: offset = 0;
0742: intLen = (value[0] == 0) ? 0 : 1;
0743:
0744: return;
0745: }
0746:
0747: if (quotient.value.length < intLen)
0748: quotient.value = new int[intLen];
0749: quotient.offset = 0;
0750: quotient.intLen = intLen;
0751:
0752: // Normalize the divisor
0753: int shift = 32 - BigInteger.bitLen(divisor);
0754:
0755: int rem = value[offset];
0756: long remLong = rem & LONG_MASK;
0757: if (remLong < divLong) {
0758: quotient.value[0] = 0;
0759: } else {
0760: quotient.value[0] = (int) (remLong / divLong);
0761: rem = (int) (remLong - (quotient.value[0] * divLong));
0762: remLong = rem & LONG_MASK;
0763: }
0764:
0765: int xlen = intLen;
0766: int[] qWord = new int[2];
0767: while (--xlen > 0) {
0768: long dividendEstimate = (remLong << 32)
0769: | (value[offset + intLen - xlen] & LONG_MASK);
0770: if (dividendEstimate >= 0) {
0771: qWord[0] = (int) (dividendEstimate / divLong);
0772: qWord[1] = (int) (dividendEstimate - (qWord[0] * divLong));
0773: } else {
0774: divWord(qWord, dividendEstimate, divisor);
0775: }
0776: quotient.value[intLen - xlen] = (int) qWord[0];
0777: rem = (int) qWord[1];
0778: remLong = rem & LONG_MASK;
0779: }
0780:
0781: // Unnormalize
0782: if (shift > 0)
0783: value[0] = rem %= divisor;
0784: else
0785: value[0] = rem;
0786: intLen = (value[0] == 0) ? 0 : 1;
0787:
0788: quotient.normalize();
0789: }
0790:
0791: /**
0792: * Calculates the quotient and remainder of this div b and places them
0793: * in the MutableBigInteger objects provided.
0794: *
0795: * Uses Algorithm D in Knuth section 4.3.1.
0796: * Many optimizations to that algorithm have been adapted from the Colin
0797: * Plumb C library.
0798: * It special cases one word divisors for speed.
0799: * The contents of a and b are not changed.
0800: *
0801: */
0802: void divide(MutableBigInteger b, MutableBigInteger quotient,
0803: MutableBigInteger rem) {
0804: if (b.intLen == 0)
0805: throw new ArithmeticException("BigInteger divide by zero");
0806:
0807: // Dividend is zero
0808: if (intLen == 0) {
0809: quotient.intLen = quotient.offset = rem.intLen = rem.offset = 0;
0810: return;
0811: }
0812:
0813: int cmp = compare(b);
0814:
0815: // Dividend less than divisor
0816: if (cmp < 0) {
0817: quotient.intLen = quotient.offset = 0;
0818: rem.copyValue(this );
0819: return;
0820: }
0821: // Dividend equal to divisor
0822: if (cmp == 0) {
0823: quotient.value[0] = quotient.intLen = 1;
0824: quotient.offset = rem.intLen = rem.offset = 0;
0825: return;
0826: }
0827:
0828: quotient.clear();
0829:
0830: // Special case one word divisor
0831: if (b.intLen == 1) {
0832: rem.copyValue(this );
0833: rem.divideOneWord(b.value[b.offset], quotient);
0834: return;
0835: }
0836:
0837: // Copy divisor value to protect divisor
0838: int[] d = new int[b.intLen];
0839: for (int i = 0; i < b.intLen; i++)
0840: d[i] = b.value[b.offset + i];
0841: int dlen = b.intLen;
0842:
0843: // Remainder starts as dividend with space for a leading zero
0844: if (rem.value.length < intLen + 1)
0845: rem.value = new int[intLen + 1];
0846:
0847: for (int i = 0; i < intLen; i++)
0848: rem.value[i + 1] = value[i + offset];
0849: rem.intLen = intLen;
0850: rem.offset = 1;
0851:
0852: int nlen = rem.intLen;
0853:
0854: // Set the quotient size
0855: int limit = nlen - dlen + 1;
0856: if (quotient.value.length < limit) {
0857: quotient.value = new int[limit];
0858: quotient.offset = 0;
0859: }
0860: quotient.intLen = limit;
0861: int[] q = quotient.value;
0862:
0863: // D1 normalize the divisor
0864: int shift = 32 - BigInteger.bitLen(d[0]);
0865: if (shift > 0) {
0866: // First shift will not grow array
0867: BigInteger.primitiveLeftShift(d, dlen, shift);
0868: // But this one might
0869: rem.leftShift(shift);
0870: }
0871:
0872: // Must insert leading 0 in rem if its length did not change
0873: if (rem.intLen == nlen) {
0874: rem.offset = 0;
0875: rem.value[0] = 0;
0876: rem.intLen++;
0877: }
0878:
0879: int dh = d[0];
0880: long dhLong = dh & LONG_MASK;
0881: int dl = d[1];
0882: int[] qWord = new int[2];
0883:
0884: // D2 Initialize j
0885: for (int j = 0; j < limit; j++) {
0886: // D3 Calculate qhat
0887: // estimate qhat
0888: int qhat = 0;
0889: int qrem = 0;
0890: boolean skipCorrection = false;
0891: int nh = rem.value[j + rem.offset];
0892: int nh2 = nh + 0x80000000;
0893: int nm = rem.value[j + 1 + rem.offset];
0894:
0895: if (nh == dh) {
0896: qhat = ~0;
0897: qrem = nh + nm;
0898: skipCorrection = qrem + 0x80000000 < nh2;
0899: } else {
0900: long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
0901: if (nChunk >= 0) {
0902: qhat = (int) (nChunk / dhLong);
0903: qrem = (int) (nChunk - (qhat * dhLong));
0904: } else {
0905: divWord(qWord, nChunk, dh);
0906: qhat = qWord[0];
0907: qrem = qWord[1];
0908: }
0909: }
0910:
0911: if (qhat == 0)
0912: continue;
0913:
0914: if (!skipCorrection) { // Correct qhat
0915: long nl = rem.value[j + 2 + rem.offset] & LONG_MASK;
0916: long rs = ((qrem & LONG_MASK) << 32) | nl;
0917: long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
0918:
0919: if (unsignedLongCompare(estProduct, rs)) {
0920: qhat--;
0921: qrem = (int) ((qrem & LONG_MASK) + dhLong);
0922: if ((qrem & LONG_MASK) >= dhLong) {
0923: estProduct = (dl & LONG_MASK)
0924: * (qhat & LONG_MASK);
0925: rs = ((qrem & LONG_MASK) << 32) | nl;
0926: if (unsignedLongCompare(estProduct, rs))
0927: qhat--;
0928: }
0929: }
0930: }
0931:
0932: // D4 Multiply and subtract
0933: rem.value[j + rem.offset] = 0;
0934: int borrow = mulsub(rem.value, d, qhat, dlen, j
0935: + rem.offset);
0936:
0937: // D5 Test remainder
0938: if (borrow + 0x80000000 > nh2) {
0939: // D6 Add back
0940: divadd(d, rem.value, j + 1 + rem.offset);
0941: qhat--;
0942: }
0943:
0944: // Store the quotient digit
0945: q[j] = qhat;
0946: } // D7 loop on j
0947:
0948: // D8 Unnormalize
0949: if (shift > 0)
0950: rem.rightShift(shift);
0951:
0952: rem.normalize();
0953: quotient.normalize();
0954: }
0955:
0956: /**
0957: * Compare two longs as if they were unsigned.
0958: * Returns true iff one is bigger than two.
0959: */
0960: private boolean unsignedLongCompare(long one, long two) {
0961: return (one + Long.MIN_VALUE) > (two + Long.MIN_VALUE);
0962: }
0963:
0964: /**
0965: * This method divides a long quantity by an int to estimate
0966: * qhat for two multi precision numbers. It is used when
0967: * the signed value of n is less than zero.
0968: */
0969: private void divWord(int[] result, long n, int d) {
0970: long dLong = d & LONG_MASK;
0971:
0972: if (dLong == 1) {
0973: result[0] = (int) n;
0974: result[1] = 0;
0975: return;
0976: }
0977:
0978: // Approximate the quotient and remainder
0979: long q = (n >>> 1) / (dLong >>> 1);
0980: long r = n - q * dLong;
0981:
0982: // Correct the approximation
0983: while (r < 0) {
0984: r += dLong;
0985: q--;
0986: }
0987: while (r >= dLong) {
0988: r -= dLong;
0989: q++;
0990: }
0991:
0992: // n - q*dlong == r && 0 <= r <dLong, hence we're done.
0993: result[0] = (int) q;
0994: result[1] = (int) r;
0995: }
0996:
0997: /**
0998: * Calculate GCD of this and b. This and b are changed by the computation.
0999: */
1000: MutableBigInteger hybridGCD(MutableBigInteger b) {
1001: // Use Euclid's algorithm until the numbers are approximately the
1002: // same length, then use the binary GCD algorithm to find the GCD.
1003: MutableBigInteger a = this ;
1004: MutableBigInteger q = new MutableBigInteger(), r = new MutableBigInteger();
1005:
1006: while (b.intLen != 0) {
1007: if (Math.abs(a.intLen - b.intLen) < 2)
1008: return a.binaryGCD(b);
1009:
1010: a.divide(b, q, r);
1011: MutableBigInteger swapper = a;
1012: a = b;
1013: b = r;
1014: r = swapper;
1015: }
1016: return a;
1017: }
1018:
1019: /**
1020: * Calculate GCD of this and v.
1021: * Assumes that this and v are not zero.
1022: */
1023: private MutableBigInteger binaryGCD(MutableBigInteger v) {
1024: // Algorithm B from Knuth section 4.5.2
1025: MutableBigInteger u = this ;
1026: MutableBigInteger r = new MutableBigInteger();
1027:
1028: // step B1
1029: int s1 = u.getLowestSetBit();
1030: int s2 = v.getLowestSetBit();
1031: int k = (s1 < s2) ? s1 : s2;
1032: if (k != 0) {
1033: u.rightShift(k);
1034: v.rightShift(k);
1035: }
1036:
1037: // step B2
1038: boolean uOdd = (k == s1);
1039: MutableBigInteger t = uOdd ? v : u;
1040: int tsign = uOdd ? -1 : 1;
1041:
1042: int lb;
1043: while ((lb = t.getLowestSetBit()) >= 0) {
1044: // steps B3 and B4
1045: t.rightShift(lb);
1046: // step B5
1047: if (tsign > 0)
1048: u = t;
1049: else
1050: v = t;
1051:
1052: // Special case one word numbers
1053: if (u.intLen < 2 && v.intLen < 2) {
1054: int x = u.value[u.offset];
1055: int y = v.value[v.offset];
1056: x = binaryGcd(x, y);
1057: r.value[0] = x;
1058: r.intLen = 1;
1059: r.offset = 0;
1060: if (k > 0)
1061: r.leftShift(k);
1062: return r;
1063: }
1064:
1065: // step B6
1066: if ((tsign = u.difference(v)) == 0)
1067: break;
1068: t = (tsign >= 0) ? u : v;
1069: }
1070:
1071: if (k > 0)
1072: u.leftShift(k);
1073: return u;
1074: }
1075:
1076: /**
1077: * Calculate GCD of a and b interpreted as unsigned integers.
1078: */
1079: static int binaryGcd(int a, int b) {
1080: if (b == 0)
1081: return a;
1082: if (a == 0)
1083: return b;
1084:
1085: int x;
1086: int aZeros = 0;
1087: while ((x = (int) a & 0xff) == 0) {
1088: a >>>= 8;
1089: aZeros += 8;
1090: }
1091: int y = BigInteger.trailingZeroTable[x];
1092: aZeros += y;
1093: a >>>= y;
1094:
1095: int bZeros = 0;
1096: while ((x = (int) b & 0xff) == 0) {
1097: b >>>= 8;
1098: bZeros += 8;
1099: }
1100: y = BigInteger.trailingZeroTable[x];
1101: bZeros += y;
1102: b >>>= y;
1103:
1104: int t = (aZeros < bZeros ? aZeros : bZeros);
1105:
1106: while (a != b) {
1107: if ((a + 0x80000000) > (b + 0x80000000)) { // a > b as unsigned
1108: a -= b;
1109:
1110: while ((x = (int) a & 0xff) == 0)
1111: a >>>= 8;
1112: a >>>= BigInteger.trailingZeroTable[x];
1113: } else {
1114: b -= a;
1115:
1116: while ((x = (int) b & 0xff) == 0)
1117: b >>>= 8;
1118: b >>>= BigInteger.trailingZeroTable[x];
1119: }
1120: }
1121: return a << t;
1122: }
1123:
1124: /**
1125: * Returns the modInverse of this mod p. This and p are not affected by
1126: * the operation.
1127: */
1128: MutableBigInteger mutableModInverse(MutableBigInteger p) {
1129: // Modulus is odd, use Schroeppel's algorithm
1130: if (p.isOdd())
1131: return modInverse(p);
1132:
1133: // Base and modulus are even, throw exception
1134: if (isEven())
1135: throw new ArithmeticException("BigInteger not invertible.");
1136:
1137: // Get even part of modulus expressed as a power of 2
1138: int powersOf2 = p.getLowestSetBit();
1139:
1140: // Construct odd part of modulus
1141: MutableBigInteger oddMod = new MutableBigInteger(p);
1142: oddMod.rightShift(powersOf2);
1143:
1144: if (oddMod.isOne())
1145: return modInverseMP2(powersOf2);
1146:
1147: // Calculate 1/a mod oddMod
1148: MutableBigInteger oddPart = modInverse(oddMod);
1149:
1150: // Calculate 1/a mod evenMod
1151: MutableBigInteger evenPart = modInverseMP2(powersOf2);
1152:
1153: // Combine the results using Chinese Remainder Theorem
1154: MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
1155: MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);
1156:
1157: MutableBigInteger temp1 = new MutableBigInteger();
1158: MutableBigInteger temp2 = new MutableBigInteger();
1159: MutableBigInteger result = new MutableBigInteger();
1160:
1161: oddPart.leftShift(powersOf2);
1162: oddPart.multiply(y1, result);
1163:
1164: evenPart.multiply(oddMod, temp1);
1165: temp1.multiply(y2, temp2);
1166:
1167: result.add(temp2);
1168: result.divide(p, temp1, temp2);
1169: return temp2;
1170: }
1171:
1172: /*
1173: * Calculate the multiplicative inverse of this mod 2^k.
1174: */
1175: MutableBigInteger modInverseMP2(int k) {
1176: if (isEven())
1177: throw new ArithmeticException("Non-invertible. (GCD != 1)");
1178:
1179: if (k > 64)
1180: return euclidModInverse(k);
1181:
1182: int t = inverseMod32(value[offset + intLen - 1]);
1183:
1184: if (k < 33) {
1185: t = (k == 32 ? t : t & ((1 << k) - 1));
1186: return new MutableBigInteger(t);
1187: }
1188:
1189: long pLong = (value[offset + intLen - 1] & LONG_MASK);
1190: if (intLen > 1)
1191: pLong |= ((long) value[offset + intLen - 2] << 32);
1192: long tLong = t & LONG_MASK;
1193: tLong = tLong * (2 - pLong * tLong); // 1 more Newton iter step
1194: tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
1195:
1196: MutableBigInteger result = new MutableBigInteger(new int[2]);
1197: result.value[0] = (int) (tLong >>> 32);
1198: result.value[1] = (int) tLong;
1199: result.intLen = 2;
1200: result.normalize();
1201: return result;
1202: }
1203:
1204: /*
1205: * Returns the multiplicative inverse of val mod 2^32. Assumes val is odd.
1206: */
1207: static int inverseMod32(int val) {
1208: // Newton's iteration!
1209: int t = val;
1210: t *= 2 - val * t;
1211: t *= 2 - val * t;
1212: t *= 2 - val * t;
1213: t *= 2 - val * t;
1214: return t;
1215: }
1216:
1217: /*
1218: * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
1219: */
1220: static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
1221: // Copy the mod to protect original
1222: return fixup(new MutableBigInteger(1), new MutableBigInteger(
1223: mod), k);
1224: }
1225:
1226: /**
1227: * Calculate the multiplicative inverse of this mod mod, where mod is odd.
1228: * This and mod are not changed by the calculation.
1229: *
1230: * This method implements an algorithm due to Richard Schroeppel, that uses
1231: * the same intermediate representation as Montgomery Reduction
1232: * ("Montgomery Form"). The algorithm is described in an unpublished
1233: * manuscript entitled "Fast Modular Reciprocals."
1234: */
1235: private MutableBigInteger modInverse(MutableBigInteger mod) {
1236: MutableBigInteger p = new MutableBigInteger(mod);
1237: MutableBigInteger f = new MutableBigInteger(this );
1238: MutableBigInteger g = new MutableBigInteger(p);
1239: SignedMutableBigInteger c = new SignedMutableBigInteger(1);
1240: SignedMutableBigInteger d = new SignedMutableBigInteger();
1241: MutableBigInteger temp = null;
1242: SignedMutableBigInteger sTemp = null;
1243:
1244: int k = 0;
1245: // Right shift f k times until odd, left shift d k times
1246: if (f.isEven()) {
1247: int trailingZeros = f.getLowestSetBit();
1248: f.rightShift(trailingZeros);
1249: d.leftShift(trailingZeros);
1250: k = trailingZeros;
1251: }
1252:
1253: // The Almost Inverse Algorithm
1254: while (!f.isOne()) {
1255: // If gcd(f, g) != 1, number is not invertible modulo mod
1256: if (f.isZero())
1257: throw new ArithmeticException(
1258: "BigInteger not invertible.");
1259:
1260: // If f < g exchange f, g and c, d
1261: if (f.compare(g) < 0) {
1262: temp = f;
1263: f = g;
1264: g = temp;
1265: sTemp = d;
1266: d = c;
1267: c = sTemp;
1268: }
1269:
1270: // If f == g (mod 4)
1271: if (((f.value[f.offset + f.intLen - 1] ^ g.value[g.offset
1272: + g.intLen - 1]) & 3) == 0) {
1273: f.subtract(g);
1274: c.signedSubtract(d);
1275: } else { // If f != g (mod 4)
1276: f.add(g);
1277: c.signedAdd(d);
1278: }
1279:
1280: // Right shift f k times until odd, left shift d k times
1281: int trailingZeros = f.getLowestSetBit();
1282: f.rightShift(trailingZeros);
1283: d.leftShift(trailingZeros);
1284: k += trailingZeros;
1285: }
1286:
1287: while (c.sign < 0)
1288: c.signedAdd(p);
1289:
1290: return fixup(c, p, k);
1291: }
1292:
1293: /*
1294: * The Fixup Algorithm
1295: * Calculates X such that X = C * 2^(-k) (mod P)
1296: * Assumes C<P and P is odd.
1297: */
1298: static MutableBigInteger fixup(MutableBigInteger c,
1299: MutableBigInteger p, int k) {
1300: MutableBigInteger temp = new MutableBigInteger();
1301: // Set r to the multiplicative inverse of p mod 2^32
1302: int r = -inverseMod32(p.value[p.offset + p.intLen - 1]);
1303:
1304: for (int i = 0, numWords = k >> 5; i < numWords; i++) {
1305: // V = R * c (mod 2^j)
1306: int v = r * c.value[c.offset + c.intLen - 1];
1307: // c = c + (v * p)
1308: p.mul(v, temp);
1309: c.add(temp);
1310: // c = c / 2^j
1311: c.intLen--;
1312: }
1313: int numBits = k & 0x1f;
1314: if (numBits != 0) {
1315: // V = R * c (mod 2^j)
1316: int v = r * c.value[c.offset + c.intLen - 1];
1317: v &= ((1 << numBits) - 1);
1318: // c = c + (v * p)
1319: p.mul(v, temp);
1320: c.add(temp);
1321: // c = c / 2^j
1322: c.rightShift(numBits);
1323: }
1324:
1325: // In theory, c may be greater than p at this point (Very rare!)
1326: while (c.compare(p) >= 0)
1327: c.subtract(p);
1328:
1329: return c;
1330: }
1331:
1332: /**
1333: * Uses the extended Euclidean algorithm to compute the modInverse of base
1334: * mod a modulus that is a power of 2. The modulus is 2^k.
1335: */
1336: MutableBigInteger euclidModInverse(int k) {
1337: MutableBigInteger b = new MutableBigInteger(1);
1338: b.leftShift(k);
1339: MutableBigInteger mod = new MutableBigInteger(b);
1340:
1341: MutableBigInteger a = new MutableBigInteger(this );
1342: MutableBigInteger q = new MutableBigInteger();
1343: MutableBigInteger r = new MutableBigInteger();
1344:
1345: b.divide(a, q, r);
1346: MutableBigInteger swapper = b;
1347: b = r;
1348: r = swapper;
1349:
1350: MutableBigInteger t1 = new MutableBigInteger(q);
1351: MutableBigInteger t0 = new MutableBigInteger(1);
1352: MutableBigInteger temp = new MutableBigInteger();
1353:
1354: while (!b.isOne()) {
1355: a.divide(b, q, r);
1356:
1357: if (r.intLen == 0)
1358: throw new ArithmeticException(
1359: "BigInteger not invertible.");
1360:
1361: swapper = r;
1362: r = a;
1363: a = swapper;
1364:
1365: if (q.intLen == 1)
1366: t1.mul(q.value[q.offset], temp);
1367: else
1368: q.multiply(t1, temp);
1369: swapper = q;
1370: q = temp;
1371: temp = swapper;
1372:
1373: t0.add(q);
1374:
1375: if (a.isOne())
1376: return t0;
1377:
1378: b.divide(a, q, r);
1379:
1380: if (r.intLen == 0)
1381: throw new ArithmeticException(
1382: "BigInteger not invertible.");
1383:
1384: swapper = b;
1385: b = r;
1386: r = swapper;
1387:
1388: if (q.intLen == 1)
1389: t0.mul(q.value[q.offset], temp);
1390: else
1391: q.multiply(t0, temp);
1392: swapper = q;
1393: q = temp;
1394: temp = swapper;
1395:
1396: t1.add(q);
1397: }
1398: mod.subtract(t1);
1399: return mod;
1400: }
1401:
1402: }
|
