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