0001: // Copyright (c) 1997, 1998, 1999, 2000, 2006 Per M.A. Bothner.
0002: // This is free software; for terms and warranty disclaimer see ./COPYING.
0003:
0004: package gnu.math;
0005:
0006: import java.io.*;
0007: import java.math.BigDecimal;
0008: import java.math.BigInteger;
0009:
0010: /** A class for infinite-precision integers.
0011: * @author Per Bothner
0012: */
0013:
0014: public class IntNum extends RatNum implements Externalizable {
0015: /** All integers are stored in 2's-complement form.
0016: * If words == null, the ival is the value of this IntNum.
0017: * Otherwise, the first ival elements of words make the value
0018: * of this IntNum, stored in little-endian order, 2's-complement form. */
0019: public int ival;
0020: public int[] words;
0021:
0022: /** We pre-allocate integers in the range minFixNum..maxFixNum. */
0023: static final int minFixNum = -100;
0024: static final int maxFixNum = 1024;
0025: static final int numFixNum = maxFixNum - minFixNum + 1;
0026: static final IntNum[] smallFixNums = new IntNum[numFixNum];
0027:
0028: static {
0029: for (int i = numFixNum; --i >= 0;)
0030: smallFixNums[i] = new IntNum(i + minFixNum);
0031: }
0032:
0033: public IntNum() {
0034: }
0035:
0036: /** Create a new (non-shared) IntNum, and initialize to an int.
0037: * @param value the initial value */
0038: public IntNum(int value) {
0039: ival = value;
0040: }
0041:
0042: /** Return a (possibly-shared) IntNum with a given int value. */
0043: public static IntNum make(int value) {
0044: if (value >= minFixNum && value <= maxFixNum)
0045: return smallFixNums[(int) value - minFixNum];
0046: else
0047: return new IntNum(value);
0048: }
0049:
0050: public static final IntNum zero() {
0051: return smallFixNums[-minFixNum];
0052: }
0053:
0054: public static final IntNum one() {
0055: return smallFixNums[1 - minFixNum];
0056: }
0057:
0058: public static final IntNum ten() {
0059: return smallFixNums[10 - minFixNum];
0060: }
0061:
0062: /** Return the IntNum for -1. */
0063: public static IntNum minusOne() {
0064: return smallFixNums[-1 - minFixNum];
0065: }
0066:
0067: /** Return a (possibly-shared) IntNum with a given long value. */
0068: public static IntNum make(long value) {
0069: if (value >= minFixNum && value <= maxFixNum)
0070: return smallFixNums[(int) value - minFixNum];
0071: int i = (int) value;
0072: if ((long) i == value)
0073: return new IntNum(i);
0074: IntNum result = alloc(2);
0075: result.ival = 2;
0076: result.words[0] = i;
0077: result.words[1] = (int) (value >> 32);
0078: return result;
0079: }
0080:
0081: /** Make an IntNum from an unsigned 64-bit value. */
0082: public static IntNum makeU(long value) {
0083: if (value >= 0)
0084: return make(value);
0085: IntNum result = alloc(3);
0086: result.ival = 3;
0087: result.words[0] = (int) value;
0088: result.words[1] = (int) (value >> 32);
0089: result.words[2] = 0;
0090: return result;
0091: }
0092:
0093: /** Make a canonicalized IntNum from an array of words.
0094: * The array may be reused (without copying). */
0095: public static IntNum make(int[] words, int len) {
0096: if (words == null)
0097: return make(len);
0098: len = IntNum.wordsNeeded(words, len);
0099: if (len <= 1)
0100: return len == 0 ? zero() : make(words[0]);
0101: IntNum num = new IntNum();
0102: num.words = words;
0103: num.ival = len;
0104: return num;
0105: }
0106:
0107: public static IntNum make(int[] words) {
0108: return make(words, words.length);
0109: }
0110:
0111: /** Allocate a new non-shared IntNum.
0112: * @param nwords number of words to allocate
0113: */
0114: public static IntNum alloc(int nwords) {
0115: if (nwords <= 1)
0116: return new IntNum();
0117: IntNum result = new IntNum();
0118: result.words = new int[nwords];
0119: return result;
0120: }
0121:
0122: /** Change words.length to nwords.
0123: * We allow words.length to be upto nwords+2 without reallocating. */
0124: public void realloc(int nwords) {
0125: if (nwords == 0) {
0126: if (words != null) {
0127: if (ival > 0)
0128: ival = words[0];
0129: words = null;
0130: }
0131: } else if (words == null || words.length < nwords
0132: || words.length > nwords + 2) {
0133: int[] new_words = new int[nwords];
0134: if (words == null) {
0135: new_words[0] = ival;
0136: ival = 1;
0137: } else {
0138: if (nwords < ival)
0139: ival = nwords;
0140: System.arraycopy(words, 0, new_words, 0, ival);
0141: }
0142: words = new_words;
0143: }
0144: }
0145:
0146: public final IntNum numerator() {
0147: return this ;
0148: }
0149:
0150: public final IntNum denominator() {
0151: return one();
0152: }
0153:
0154: public final boolean isNegative() {
0155: return (words == null ? ival : words[ival - 1]) < 0;
0156: }
0157:
0158: public int sign() {
0159: int n = ival;
0160: int[] w = words;
0161: if (w == null)
0162: return n > 0 ? 1 : n < 0 ? -1 : 0;
0163: int i = w[--n];
0164: if (i > 0)
0165: return 1;
0166: if (i < 0)
0167: return -1;
0168: for (;;) {
0169: if (n == 0)
0170: return 0;
0171: if (w[--n] != 0)
0172: return 1;
0173: }
0174: }
0175:
0176: /** Return -1, 0, or 1, depending on which value is greater. */
0177: public static int compare(IntNum x, IntNum y) {
0178: if (x.words == null && y.words == null)
0179: return x.ival < y.ival ? -1 : x.ival > y.ival ? 1 : 0;
0180: boolean x_negative = x.isNegative();
0181: boolean y_negative = y.isNegative();
0182: if (x_negative != y_negative)
0183: return x_negative ? -1 : 1;
0184: int x_len = x.words == null ? 1 : x.ival;
0185: int y_len = y.words == null ? 1 : y.ival;
0186: if (x_len != y_len) // We assume x and y are canonicalized.
0187: return (x_len > y_len) != x_negative ? 1 : -1;
0188: return MPN.cmp(x.words, y.words, x_len);
0189: }
0190:
0191: /** Return -1, 0, or 1, depending on which value is greater. */
0192: public static int compare(IntNum x, long y) {
0193: long x_word;
0194: if (x.words == null)
0195: x_word = x.ival;
0196: else {
0197: boolean x_negative = x.isNegative();
0198: boolean y_negative = y < 0;
0199: if (x_negative != y_negative)
0200: return x_negative ? -1 : 1;
0201: int x_len = x.words == null ? 1 : x.ival;
0202: if (x_len == 1)
0203: x_word = x.words[0];
0204: else if (x_len == 2)
0205: x_word = x.longValue();
0206: else
0207: // We assume x is canonicalized.
0208: return x_negative ? -1 : 1;
0209: }
0210: return x_word < y ? -1 : x_word > y ? 1 : 0;
0211: }
0212:
0213: public int compare(Object obj) {
0214: if (obj instanceof IntNum)
0215: return compare(this , (IntNum) obj);
0216: return ((RealNum) obj).compareReversed(this );
0217: }
0218:
0219: public final boolean isOdd() {
0220: int low = words == null ? ival : words[0];
0221: return (low & 1) != 0;
0222: }
0223:
0224: public final boolean isZero() {
0225: return words == null && ival == 0;
0226: }
0227:
0228: public final boolean isOne() {
0229: return words == null && ival == 1;
0230: }
0231:
0232: public final boolean isMinusOne() {
0233: return words == null && ival == -1;
0234: }
0235:
0236: /** Calculate how many words are significant in words[0:len-1].
0237: * Returns the least value x such that x>0 && words[0:x-1]==words[0:len-1],
0238: * when words is viewed as a 2's complement integer.
0239: */
0240: public static int wordsNeeded(int[] words, int len) {
0241: int i = len;
0242: if (i > 0) {
0243: int word = words[--i];
0244: if (word == -1) {
0245: while (i > 0 && (word = words[i - 1]) < 0) {
0246: i--;
0247: if (word != -1)
0248: break;
0249: }
0250: } else {
0251: while (word == 0 && i > 0 && (word = words[i - 1]) >= 0)
0252: i--;
0253: }
0254: }
0255: return i + 1;
0256: }
0257:
0258: public IntNum canonicalize() {
0259: if (words != null
0260: && (ival = IntNum.wordsNeeded(words, ival)) <= 1) {
0261: if (ival == 1)
0262: ival = words[0];
0263: words = null;
0264: }
0265: if (words == null && ival >= minFixNum && ival <= maxFixNum)
0266: return smallFixNums[(int) ival - minFixNum];
0267: return this ;
0268: }
0269:
0270: /** Add two ints, yielding an IntNum. */
0271: public static final IntNum add(int x, int y) {
0272: return IntNum.make((long) x + (long) y);
0273: }
0274:
0275: /** Add an IntNum and an int, yielding a new IntNum. */
0276: public static IntNum add(IntNum x, int y) {
0277: if (x.words == null)
0278: return IntNum.add(x.ival, y);
0279: IntNum result = new IntNum(0);
0280: result.setAdd(x, y);
0281: return result.canonicalize();
0282: }
0283:
0284: /** Set this to the sum of x and y.
0285: * OK if x==this. */
0286: public void setAdd(IntNum x, int y) {
0287: if (x.words == null) {
0288: set((long) x.ival + (long) y);
0289: return;
0290: }
0291: int len = x.ival;
0292: realloc(len + 1);
0293: long carry = y;
0294: for (int i = 0; i < len; i++) {
0295: carry += ((long) x.words[i] & 0xffffffffL);
0296: words[i] = (int) carry;
0297: carry >>= 32;
0298: }
0299: if (x.words[len - 1] < 0)
0300: carry--;
0301: words[len] = (int) carry;
0302: ival = wordsNeeded(words, len + 1);
0303: }
0304:
0305: /** Destructively add an int to this. */
0306: public final void setAdd(int y) {
0307: setAdd(this , y);
0308: }
0309:
0310: /** Destructively set the value of this to an int. */
0311: public final void set(int y) {
0312: words = null;
0313: ival = y;
0314: }
0315:
0316: /** Destructively set the value of this to a long. */
0317: public final void set(long y) {
0318: int i = (int) y;
0319: if ((long) i == y) {
0320: ival = i;
0321: words = null;
0322: } else {
0323: realloc(2);
0324: words[0] = i;
0325: words[1] = (int) (y >> 32);
0326: ival = 2;
0327: }
0328: }
0329:
0330: /** Destructively set the value of this to the given words.
0331: * The words array is reused, not copied. */
0332: public final void set(int[] words, int length) {
0333: this .ival = length;
0334: this .words = words;
0335: }
0336:
0337: /** Destructively set the value of this to that of y. */
0338: public final void set(IntNum y) {
0339: if (y.words == null)
0340: set(y.ival);
0341: else if (this != y) {
0342: realloc(y.ival);
0343: System.arraycopy(y.words, 0, words, 0, y.ival);
0344: ival = y.ival;
0345: }
0346: }
0347:
0348: /** Add two IntNums, yielding their sum as another IntNum. */
0349: public static IntNum add(IntNum x, IntNum y) {
0350: return add(x, y, 1);
0351: }
0352:
0353: /** Subtract two IntNums, yielding their sum as another IntNum. */
0354: public static IntNum sub(IntNum x, IntNum y) {
0355: return add(x, y, -1);
0356: }
0357:
0358: /** Add two IntNums, yielding their sum as another IntNum. */
0359: public static IntNum add(IntNum x, IntNum y, int k) {
0360: if (x.words == null && y.words == null)
0361: return IntNum
0362: .make((long) k * (long) y.ival + (long) x.ival);
0363: if (k != 1) {
0364: if (k == -1)
0365: y = IntNum.neg(y);
0366: else
0367: y = IntNum.times(y, IntNum.make(k));
0368: }
0369: if (x.words == null)
0370: return IntNum.add(y, x.ival);
0371: if (y.words == null)
0372: return IntNum.add(x, y.ival);
0373: // Both are big
0374: int len;
0375: if (y.ival > x.ival) { // Swap so x is longer then y.
0376: IntNum tmp = x;
0377: x = y;
0378: y = tmp;
0379: }
0380: IntNum result = alloc(x.ival + 1);
0381: int i = y.ival;
0382: long carry = MPN.add_n(result.words, x.words, y.words, i);
0383: long y_ext = y.words[i - 1] < 0 ? 0xffffffffL : 0;
0384: for (; i < x.ival; i++) {
0385: carry += ((long) x.words[i] & 0xffffffffL) + y_ext;
0386: ;
0387: result.words[i] = (int) carry;
0388: carry >>>= 32;
0389: }
0390: if (x.words[i - 1] < 0)
0391: y_ext--;
0392: result.words[i] = (int) (carry + y_ext);
0393: result.ival = i + 1;
0394: return result.canonicalize();
0395: }
0396:
0397: /** Multiply two ints, yielding an IntNum. */
0398: public static final IntNum times(int x, int y) {
0399: return IntNum.make((long) x * (long) y);
0400: }
0401:
0402: public static final IntNum times(IntNum x, int y) {
0403: if (y == 0)
0404: return zero();
0405: if (y == 1)
0406: return x;
0407: int[] xwords = x.words;
0408: int xlen = x.ival;
0409: if (xwords == null)
0410: return IntNum.make((long) xlen * (long) y);
0411: boolean negative;
0412: IntNum result = IntNum.alloc(xlen + 1);
0413: if (xwords[xlen - 1] < 0) {
0414: negative = true;
0415: negate(result.words, xwords, xlen);
0416: xwords = result.words;
0417: } else
0418: negative = false;
0419: if (y < 0) {
0420: negative = !negative;
0421: y = -y;
0422: }
0423: result.words[xlen] = MPN.mul_1(result.words, xwords, xlen, y);
0424: result.ival = xlen + 1;
0425: if (negative)
0426: result.setNegative();
0427: return result.canonicalize();
0428: }
0429:
0430: public static final IntNum times(IntNum x, IntNum y) {
0431: if (y.words == null)
0432: return times(x, y.ival);
0433: if (x.words == null)
0434: return times(y, x.ival);
0435: boolean negative = false;
0436: int[] xwords;
0437: int[] ywords;
0438: int xlen = x.ival;
0439: int ylen = y.ival;
0440: if (x.isNegative()) {
0441: negative = true;
0442: xwords = new int[xlen];
0443: negate(xwords, x.words, xlen);
0444: } else {
0445: negative = false;
0446: xwords = x.words;
0447: }
0448: if (y.isNegative()) {
0449: negative = !negative;
0450: ywords = new int[ylen];
0451: negate(ywords, y.words, ylen);
0452: } else
0453: ywords = y.words;
0454: // Swap if x is shorter then y.
0455: if (xlen < ylen) {
0456: int[] twords = xwords;
0457: xwords = ywords;
0458: ywords = twords;
0459: int tlen = xlen;
0460: xlen = ylen;
0461: ylen = tlen;
0462: }
0463: IntNum result = IntNum.alloc(xlen + ylen);
0464: MPN.mul(result.words, xwords, xlen, ywords, ylen);
0465: result.ival = xlen + ylen;
0466: if (negative)
0467: result.setNegative();
0468: return result.canonicalize();
0469: }
0470:
0471: public static void divide(long x, long y, IntNum quotient,
0472: IntNum remainder, int rounding_mode) {
0473: boolean xNegative, yNegative;
0474: if (x < 0) {
0475: xNegative = true;
0476: if (x == Long.MIN_VALUE) {
0477: divide(IntNum.make(x), IntNum.make(y), quotient,
0478: remainder, rounding_mode);
0479: return;
0480: }
0481: x = -x;
0482: } else
0483: xNegative = false;
0484:
0485: if (y < 0) {
0486: yNegative = true;
0487: if (y == Long.MIN_VALUE) {
0488: if (rounding_mode == TRUNCATE) { // x != Long.Min_VALUE implies abs(x) < abs(y)
0489: if (quotient != null)
0490: quotient.set(0);
0491: if (remainder != null)
0492: remainder.set(x);
0493: } else
0494: divide(IntNum.make(x), IntNum.make(y), quotient,
0495: remainder, rounding_mode);
0496: return;
0497: }
0498: y = -y;
0499: } else
0500: yNegative = false;
0501:
0502: long q = x / y;
0503: long r = x % y;
0504: boolean qNegative = xNegative ^ yNegative;
0505:
0506: boolean add_one = false;
0507: if (r != 0) {
0508: switch (rounding_mode) {
0509: case TRUNCATE:
0510: break;
0511: case CEILING:
0512: case FLOOR:
0513: if (qNegative == (rounding_mode == FLOOR))
0514: add_one = true;
0515: break;
0516: case ROUND:
0517: add_one = r > ((y - (q & 1)) >> 1);
0518: break;
0519: }
0520: }
0521: if (quotient != null) {
0522: if (add_one)
0523: q++;
0524: if (qNegative)
0525: q = -q;
0526: quotient.set(q);
0527: }
0528: if (remainder != null) {
0529: // The remainder is by definition: X-Q*Y
0530: if (add_one) {
0531: // Subtract the remainder from Y.
0532: r = y - r;
0533: // In this case, abs(Q*Y) > abs(X).
0534: // So sign(remainder) = -sign(X).
0535: xNegative = !xNegative;
0536: } else {
0537: // If !add_one, then: abs(Q*Y) <= abs(X).
0538: // So sign(remainder) = sign(X).
0539: }
0540: if (xNegative)
0541: r = -r;
0542: remainder.set(r);
0543: }
0544: }
0545:
0546: /** Divide two integers, yielding quotient and remainder.
0547: * @param x the numerator in the division
0548: * @param y the denominator in the division
0549: * @param quotient is set to the quotient of the result (iff quotient!=null)
0550: * @param remainder is set to the remainder of the result
0551: * (iff remainder!=null)
0552: * @param rounding_mode one of FLOOR, CEILING, TRUNCATE, or ROUND.
0553: */
0554: public static void divide(IntNum x, IntNum y, IntNum quotient,
0555: IntNum remainder, int rounding_mode) {
0556: if ((x.words == null || x.ival <= 2)
0557: && (y.words == null || y.ival <= 2)) {
0558: long x_l = x.longValue();
0559: long y_l = y.longValue();
0560: if (x_l != Long.MIN_VALUE && y_l != Long.MIN_VALUE) {
0561: divide(x_l, y_l, quotient, remainder, rounding_mode);
0562: return;
0563: }
0564: }
0565:
0566: boolean xNegative = x.isNegative();
0567: boolean yNegative = y.isNegative();
0568: boolean qNegative = xNegative ^ yNegative;
0569:
0570: int ylen = y.words == null ? 1 : y.ival;
0571: int[] ywords = new int[ylen];
0572: y.getAbsolute(ywords);
0573: while (ylen > 1 && ywords[ylen - 1] == 0)
0574: ylen--;
0575:
0576: int xlen = x.words == null ? 1 : x.ival;
0577: int[] xwords = new int[xlen + 2];
0578: x.getAbsolute(xwords);
0579: while (xlen > 1 && xwords[xlen - 1] == 0)
0580: xlen--;
0581:
0582: int qlen, rlen;
0583:
0584: int cmpval = MPN.cmp(xwords, xlen, ywords, ylen);
0585: if (cmpval < 0) // abs(x) < abs(y)
0586: { // quotient = 0; remainder = num.
0587: int[] rwords = xwords;
0588: xwords = ywords;
0589: ywords = rwords;
0590: rlen = xlen;
0591: qlen = 1;
0592: xwords[0] = 0;
0593: } else if (cmpval == 0) // abs(x) == abs(y)
0594: {
0595: xwords[0] = 1;
0596: qlen = 1; // quotient = 1
0597: ywords[0] = 0;
0598: rlen = 1; // remainder = 0;
0599: } else if (ylen == 1) {
0600: qlen = xlen;
0601: rlen = 1;
0602: ywords[0] = MPN.divmod_1(xwords, xwords, xlen, ywords[0]);
0603: } else // abs(x) > abs(y)
0604: {
0605: // Normalize the denominator, i.e. make its most significant bit set by
0606: // shifting it normalization_steps bits to the left. Also shift the
0607: // numerator the same number of steps (to keep the quotient the same!).
0608:
0609: int nshift = MPN.count_leading_zeros(ywords[ylen - 1]);
0610: if (nshift != 0) {
0611: // Shift up the denominator setting the most significant bit of
0612: // the most significant word.
0613: MPN.lshift(ywords, 0, ywords, ylen, nshift);
0614:
0615: // Shift up the numerator, possibly introducing a new most
0616: // significant word.
0617: int x_high = MPN
0618: .lshift(xwords, 0, xwords, xlen, nshift);
0619: xwords[xlen++] = x_high;
0620: }
0621:
0622: if (xlen == ylen)
0623: xwords[xlen++] = 0;
0624: MPN.divide(xwords, xlen, ywords, ylen);
0625: rlen = ylen;
0626: MPN.rshift0(ywords, xwords, 0, rlen, nshift);
0627:
0628: qlen = xlen + 1 - ylen;
0629: if (quotient != null) {
0630: for (int i = 0; i < qlen; i++)
0631: xwords[i] = xwords[i + ylen];
0632: }
0633: }
0634:
0635: while (rlen > 1 && ywords[rlen - 1] == 0)
0636: rlen--;
0637: if (ywords[rlen - 1] < 0) {
0638: ywords[rlen] = 0;
0639: rlen++;
0640: }
0641:
0642: // Now the quotient is in xwords, and the remainder is in ywords.
0643:
0644: boolean add_one = false;
0645: if (rlen > 1 || ywords[0] != 0) { // Non-zero remainder i.e. in-exact quotient.
0646: switch (rounding_mode) {
0647: case TRUNCATE:
0648: break;
0649: case CEILING:
0650: case FLOOR:
0651: if (qNegative == (rounding_mode == FLOOR))
0652: add_one = true;
0653: break;
0654: case ROUND:
0655: // int cmp = compare (remainder<<1, abs(y));
0656: IntNum tmp = remainder == null ? new IntNum()
0657: : remainder;
0658: tmp.set(ywords, rlen);
0659: tmp = shift(tmp, 1);
0660: if (yNegative)
0661: tmp.setNegative();
0662: int cmp = compare(tmp, y);
0663: // Now cmp == compare(sign(y)*(remainder<<1), y)
0664: if (yNegative)
0665: cmp = -cmp;
0666: add_one = (cmp == 1)
0667: || (cmp == 0 && (xwords[0] & 1) != 0);
0668: }
0669: }
0670: if (quotient != null) {
0671: if (xwords[qlen - 1] < 0) {
0672: xwords[qlen] = 0;
0673: qlen++;
0674: }
0675: quotient.set(xwords, qlen);
0676: if (qNegative) {
0677: if (add_one) // -(quotient + 1) == ~(quotient)
0678: quotient.setInvert();
0679: else
0680: quotient.setNegative();
0681: } else if (add_one)
0682: quotient.setAdd(1);
0683: }
0684: if (remainder != null) {
0685: // The remainder is by definition: X-Q*Y
0686: remainder.set(ywords, rlen);
0687: if (add_one) {
0688: // Subtract the remainder from Y:
0689: // abs(R) = abs(Y) - abs(orig_rem) = -(abs(orig_rem) - abs(Y)).
0690: IntNum tmp;
0691: if (y.words == null) {
0692: tmp = remainder;
0693: tmp.set(yNegative ? ywords[0] + y.ival : ywords[0]
0694: - y.ival);
0695: } else
0696: tmp = IntNum.add(remainder, y, yNegative ? 1 : -1);
0697: // Now tmp <= 0.
0698: // In this case, abs(Q) = 1 + floor(abs(X)/abs(Y)).
0699: // Hence, abs(Q*Y) > abs(X).
0700: // So sign(remainder) = -sign(X).
0701: if (xNegative)
0702: remainder.setNegative(tmp);
0703: else
0704: remainder.set(tmp);
0705: } else {
0706: // If !add_one, then: abs(Q*Y) <= abs(X).
0707: // So sign(remainder) = sign(X).
0708: if (xNegative)
0709: remainder.setNegative();
0710: }
0711: }
0712: }
0713:
0714: public static IntNum quotient(IntNum x, IntNum y, int rounding_mode) {
0715: IntNum quotient = new IntNum();
0716: divide(x, y, quotient, null, rounding_mode);
0717: return quotient.canonicalize();
0718: }
0719:
0720: public static IntNum quotient(IntNum x, IntNum y) {
0721: return quotient(x, y, TRUNCATE);
0722: }
0723:
0724: public IntNum toExactInt(int rounding_mode) {
0725: return this ;
0726: }
0727:
0728: public RealNum toInt(int rounding_mode) {
0729: return this ;
0730: }
0731:
0732: public static IntNum remainder(IntNum x, IntNum y) {
0733: if (y.isZero())
0734: return x;
0735: IntNum rem = new IntNum();
0736: divide(x, y, null, rem, TRUNCATE);
0737: return rem.canonicalize();
0738: }
0739:
0740: public static IntNum modulo(IntNum x, IntNum y) {
0741: if (y.isZero())
0742: return x;
0743: IntNum rem = new IntNum();
0744: divide(x, y, null, rem, FLOOR);
0745: return rem.canonicalize();
0746: }
0747:
0748: public Numeric power(IntNum y) {
0749: if (isOne())
0750: return this ;
0751: if (isMinusOne())
0752: return y.isOdd() ? this : IntNum.one();
0753: if (y.words == null && y.ival >= 0)
0754: return power(this , y.ival);
0755: if (isZero())
0756: return y.isNegative() ? RatNum.infinity(-1) : (RatNum) this ;
0757: return super .power(y);
0758: }
0759:
0760: /** Calculate the integral power of an IntNum.
0761: * @param x the value (base) to exponentiate
0762: * @param y the exponent (must be non-negative)
0763: */
0764: public static IntNum power(IntNum x, int y) {
0765: if (y <= 0) {
0766: if (y == 0)
0767: return one();
0768: else
0769: throw new Error("negative exponent");
0770: }
0771: if (x.isZero())
0772: return x;
0773: int plen = x.words == null ? 1 : x.ival; // Length of pow2.
0774: int blen = ((x.intLength() * y) >> 5) + 2 * plen;
0775: boolean negative = x.isNegative() && (y & 1) != 0;
0776: int[] pow2 = new int[blen];
0777: int[] rwords = new int[blen];
0778: int[] work = new int[blen];
0779: x.getAbsolute(pow2); // pow2 = abs(x);
0780: int rlen = 1;
0781: rwords[0] = 1; // rwords = 1;
0782: for (;;) // for (i = 0; ; i++)
0783: {
0784: // pow2 == x**(2**i)
0785: // prod = x**(sum(j=0..i-1, (y>>j)&1))
0786: if ((y & 1) != 0) { // r *= pow2
0787: MPN.mul(work, pow2, plen, rwords, rlen);
0788: int[] temp = work;
0789: work = rwords;
0790: rwords = temp;
0791: rlen += plen;
0792: while (rwords[rlen - 1] == 0)
0793: rlen--;
0794: }
0795: y >>= 1;
0796: if (y == 0)
0797: break;
0798: // pow2 *= pow2;
0799: MPN.mul(work, pow2, plen, pow2, plen);
0800: int[] temp = work;
0801: work = pow2;
0802: pow2 = temp; // swap to avoid a copy
0803: plen *= 2;
0804: while (pow2[plen - 1] == 0)
0805: plen--;
0806: }
0807: if (rwords[rlen - 1] < 0)
0808: rlen++;
0809: if (negative)
0810: negate(rwords, rwords, rlen);
0811: return IntNum.make(rwords, rlen);
0812: }
0813:
0814: /** Calculate Greatest Common Divisor for non-negative ints. */
0815: public static final int gcd(int a, int b) {
0816: // Euclid's algorithm, copied from libg++.
0817: if (b > a) {
0818: int tmp = a;
0819: a = b;
0820: b = tmp;
0821: }
0822: for (;;) {
0823: if (b == 0)
0824: return a;
0825: else if (b == 1)
0826: return b;
0827: else {
0828: int tmp = b;
0829: b = a % b;
0830: a = tmp;
0831: }
0832: }
0833: }
0834:
0835: public static IntNum gcd(IntNum x, IntNum y) {
0836: int xval = x.ival;
0837: int yval = y.ival;
0838: if (x.words == null) {
0839: if (xval == 0)
0840: return IntNum.abs(y);
0841: if (y.words == null && xval != Integer.MIN_VALUE
0842: && yval != Integer.MIN_VALUE) {
0843: if (xval < 0)
0844: xval = -xval;
0845: if (yval < 0)
0846: yval = -yval;
0847: return IntNum.make(IntNum.gcd(xval, yval));
0848: }
0849: xval = 1;
0850: }
0851: if (y.words == null) {
0852: if (yval == 0)
0853: return IntNum.abs(x);
0854: yval = 1;
0855: }
0856: int len = (xval > yval ? xval : yval) + 1;
0857: int[] xwords = new int[len];
0858: int[] ywords = new int[len];
0859: x.getAbsolute(xwords);
0860: y.getAbsolute(ywords);
0861: len = MPN.gcd(xwords, ywords, len);
0862: IntNum result = new IntNum(0);
0863: result.ival = len;
0864: result.words = xwords;
0865: return result.canonicalize();
0866: }
0867:
0868: public static IntNum lcm(IntNum x, IntNum y) {
0869: if (x.isZero() || y.isZero())
0870: return IntNum.zero();
0871: x = IntNum.abs(x);
0872: y = IntNum.abs(y);
0873: IntNum quotient = new IntNum();
0874: divide(times(x, y), gcd(x, y), quotient, null, TRUNCATE);
0875: return quotient.canonicalize();
0876: }
0877:
0878: void setInvert() {
0879: if (words == null)
0880: ival = ~ival;
0881: else {
0882: for (int i = ival; --i >= 0;)
0883: words[i] = ~words[i];
0884: }
0885: }
0886:
0887: void setShiftLeft(IntNum x, int count) {
0888: int[] xwords;
0889: int xlen;
0890: if (x.words == null) {
0891: if (count < 32) {
0892: set((long) x.ival << count);
0893: return;
0894: }
0895: xwords = new int[1];
0896: xwords[0] = x.ival;
0897: xlen = 1;
0898: } else {
0899: xwords = x.words;
0900: xlen = x.ival;
0901: }
0902: int word_count = count >> 5;
0903: count &= 31;
0904: int new_len = xlen + word_count;
0905: if (count == 0) {
0906: realloc(new_len);
0907: for (int i = xlen; --i >= 0;)
0908: words[i + word_count] = xwords[i];
0909: } else {
0910: new_len++;
0911: realloc(new_len);
0912: int shift_out = MPN.lshift(words, word_count, xwords, xlen,
0913: count);
0914: count = 32 - count;
0915: words[new_len - 1] = (shift_out << count) >> count; // sign-extend.
0916: }
0917: ival = new_len;
0918: for (int i = word_count; --i >= 0;)
0919: words[i] = 0;
0920: }
0921:
0922: void setShiftRight(IntNum x, int count) {
0923: if (x.words == null)
0924: set(count < 32 ? x.ival >> count : x.ival < 0 ? -1 : 0);
0925: else if (count == 0)
0926: set(x);
0927: else {
0928: boolean neg = x.isNegative();
0929: int word_count = count >> 5;
0930: count &= 31;
0931: int d_len = x.ival - word_count;
0932: if (d_len <= 0)
0933: set(neg ? -1 : 0);
0934: else {
0935: if (words == null || words.length < d_len)
0936: realloc(d_len);
0937: MPN.rshift0(words, x.words, word_count, d_len, count);
0938: ival = d_len;
0939: if (neg)
0940: words[d_len - 1] |= -2 << (31 - count);
0941: }
0942: }
0943: }
0944:
0945: void setShift(IntNum x, int count) {
0946: if (count > 0)
0947: setShiftLeft(x, count);
0948: else
0949: setShiftRight(x, -count);
0950: }
0951:
0952: public static IntNum shift(IntNum x, int count) {
0953: if (x.words == null) {
0954: if (count <= 0)
0955: return make(count > -32 ? x.ival >> (-count)
0956: : x.ival < 0 ? -1 : 0);
0957: if (count < 32)
0958: return make((long) x.ival << count);
0959: }
0960: if (count == 0)
0961: return x;
0962: IntNum result = new IntNum(0);
0963: result.setShift(x, count);
0964: return result.canonicalize();
0965: }
0966:
0967: public void format(int radix, StringBuffer buffer) {
0968: if (words == null)
0969: buffer.append(Integer.toString(ival, radix));
0970: else if (ival <= 2)
0971: buffer.append(Long.toString(longValue(), radix));
0972: else {
0973: boolean neg = isNegative();
0974: int[] work;
0975: if (neg || radix != 16) {
0976: work = new int[ival];
0977: getAbsolute(work);
0978: } else
0979: work = words;
0980: int len = ival;
0981:
0982: if (radix == 16) {
0983: if (neg)
0984: buffer.append('-');
0985: int buf_start = buffer.length();
0986: for (int i = len; --i >= 0;) {
0987: int word = work[i];
0988: for (int j = 8; --j >= 0;) {
0989: int hex_digit = (word >> (4 * j)) & 0xF;
0990: // Suppress leading zeros:
0991: if (hex_digit > 0
0992: || buffer.length() > buf_start)
0993: buffer.append(Character.forDigit(hex_digit,
0994: 16));
0995: }
0996: }
0997: } else {
0998: int i = buffer.length();
0999: for (;;) {
1000: int digit = MPN.divmod_1(work, work, len, radix);
1001: buffer.append(Character.forDigit(digit, radix));
1002: while (len > 0 && work[len - 1] == 0)
1003: len--;
1004: if (len == 0)
1005: break;
1006: }
1007: if (neg)
1008: buffer.append('-');
1009: /* Reverse buffer. */
1010: int j = buffer.length() - 1;
1011: while (i < j) {
1012: char tmp = buffer.charAt(i);
1013: buffer.setCharAt(i, buffer.charAt(j));
1014: buffer.setCharAt(j, tmp);
1015: i++;
1016: j--;
1017: }
1018: }
1019: }
1020: }
1021:
1022: public String toString(int radix) {
1023: if (words == null)
1024: return Integer.toString(ival, radix);
1025: else if (ival <= 2)
1026: return Long.toString(longValue(), radix);
1027: int buf_size = ival * (MPN.chars_per_word(radix) + 1);
1028: StringBuffer buffer = new StringBuffer(buf_size);
1029: format(radix, buffer);
1030: return buffer.toString();
1031: }
1032:
1033: public int intValue() {
1034: if (words == null)
1035: return ival;
1036: return words[0];
1037: }
1038:
1039: /** Cast an Object to an int. The Object must (currently) be an IntNum. */
1040: public static int intValue(Object obj) {
1041: IntNum inum = (IntNum) obj;
1042: if (inum.words != null)
1043: // This is not quite appropriate, but will do.
1044: throw new ClassCastException("integer too large");
1045: return inum.ival;
1046: }
1047:
1048: public long longValue() {
1049: if (words == null)
1050: return ival;
1051: if (ival == 1)
1052: return words[0];
1053: return ((long) words[1] << 32)
1054: + ((long) words[0] & 0xffffffffL);
1055: }
1056:
1057: public int hashCode() {
1058: return words == null ? ival : (words[0] + words[ival - 1]);
1059: }
1060:
1061: /* Assumes x and y are both canonicalized. */
1062: public static boolean equals(IntNum x, IntNum y) {
1063: if (x.words == null && y.words == null)
1064: return x.ival == y.ival;
1065: if (x.words == null || y.words == null || x.ival != y.ival)
1066: return false;
1067: for (int i = x.ival; --i >= 0;) {
1068: if (x.words[i] != y.words[i])
1069: return false;
1070: }
1071: return true;
1072: }
1073:
1074: /* Assumes this and obj are both canonicalized. */
1075: public boolean equals(Object obj) {
1076: if (obj == null || !(obj instanceof IntNum))
1077: return false;
1078: return IntNum.equals(this , (IntNum) obj);
1079: }
1080:
1081: public static IntNum valueOf(char[] buf, int offset, int length,
1082: int radix, boolean negative) {
1083: int byte_len = 0;
1084: byte[] bytes = new byte[length];
1085: for (int i = 0; i < length; i++) {
1086: char ch = buf[offset + i];
1087: if (ch == '-')
1088: negative = true;
1089: else if (ch == '_'
1090: || (byte_len == 0 && (ch == ' ' || ch == '\t')))
1091: continue;
1092: else {
1093: int digit = Character.digit(ch, radix);
1094: if (digit < 0)
1095: break;
1096: bytes[byte_len++] = (byte) digit;
1097: }
1098: }
1099: return valueOf(bytes, byte_len, negative, radix);
1100: }
1101:
1102: public static IntNum valueOf(String s, int radix)
1103: throws NumberFormatException {
1104: int len = s.length();
1105: // Testing (len < 2 * MPN.chars_per_word(radix)) would be more accurate,
1106: // but slightly more expensive, for little practical gain.
1107: if (len + radix <= 28)
1108: return IntNum.make(Long.parseLong(s, radix));
1109:
1110: int byte_len = 0;
1111: byte[] bytes = new byte[len];
1112: boolean negative = false;
1113: for (int i = 0; i < len; i++) {
1114: char ch = s.charAt(i);
1115: if (ch == '-')
1116: negative = true;
1117: else if (ch == '_'
1118: || (byte_len == 0 && (ch == ' ' || ch == '\t')))
1119: continue;
1120: else {
1121: int digit = Character.digit(ch, radix);
1122: if (digit < 0)
1123: throw new NumberFormatException(
1124: "For input string: \"" + s + '"');
1125: bytes[byte_len++] = (byte) digit;
1126: }
1127: }
1128: return valueOf(bytes, byte_len, negative, radix);
1129: }
1130:
1131: public static IntNum valueOf(byte[] digits, int byte_len,
1132: boolean negative, int radix) {
1133: int chars_per_word = MPN.chars_per_word(radix);
1134: int[] words = new int[byte_len / chars_per_word + 1];
1135: int size = MPN.set_str(words, digits, byte_len, radix);
1136: if (size == 0)
1137: return zero();
1138: if (words[size - 1] < 0)
1139: words[size++] = 0;
1140: if (negative)
1141: negate(words, words, size);
1142: return make(words, size);
1143: }
1144:
1145: public static IntNum valueOf(String s) throws NumberFormatException {
1146: return IntNum.valueOf(s, 10);
1147: }
1148:
1149: public double doubleValue() {
1150: if (words == null)
1151: return (double) ival;
1152: if (ival <= 2)
1153: return (double) longValue();
1154: if (isNegative())
1155: return IntNum.neg(this ).roundToDouble(0, true, false);
1156: else
1157: return roundToDouble(0, false, false);
1158: }
1159:
1160: /** Return true if any of the lowest n bits are one.
1161: * (false if n is negative). */
1162: boolean checkBits(int n) {
1163: if (n <= 0)
1164: return false;
1165: if (words == null)
1166: return n > 31 || ((ival & ((1 << n) - 1)) != 0);
1167: int i;
1168: for (i = 0; i < (n >> 5); i++)
1169: if (words[i] != 0)
1170: return true;
1171: return (n & 31) != 0 && (words[i] & ((1 << (n & 31)) - 1)) != 0;
1172: }
1173:
1174: /** Convert a semi-processed IntNum to double.
1175: * Number must be non-negative. Multiplies by a power of two, applies sign,
1176: * and converts to double, with the usual java rounding.
1177: * @param exp power of two, positive or negative, by which to multiply
1178: * @param neg true if negative
1179: * @param remainder true if the IntNum is the result of a truncating
1180: * division that had non-zero remainder. To ensure proper rounding in
1181: * this case, the IntNum must have at least 54 bits. */
1182: public double roundToDouble(int exp, boolean neg, boolean remainder) {
1183: // Compute length.
1184: int il = intLength();
1185:
1186: // Exponent when normalized to have decimal point directly after
1187: // leading one. This is stored excess 1023 in the exponent bit field.
1188: exp += il - 1;
1189:
1190: // Gross underflow. If exp == -1075, we let the rounding
1191: // computation determine whether it is minval or 0 (which are just
1192: // 0x0000 0000 0000 0001 and 0x0000 0000 0000 0000 as bit
1193: // patterns).
1194: if (exp < -1075)
1195: return neg ? -0.0 : 0.0;
1196:
1197: // gross overflow
1198: if (exp > 1023)
1199: return neg ? Double.NEGATIVE_INFINITY
1200: : Double.POSITIVE_INFINITY;
1201:
1202: // number of bits in mantissa, including the leading one.
1203: // 53 unless it's denormalized
1204: int ml = (exp >= -1022 ? 53 : 53 + exp + 1022);
1205:
1206: // Get top ml + 1 bits. The extra one is for rounding.
1207: long m;
1208: int excess_bits = il - (ml + 1);
1209: if (excess_bits > 0)
1210: m = ((words == null) ? ival >> excess_bits : MPN
1211: .rshift_long(words, ival, excess_bits));
1212: else
1213: m = longValue() << (-excess_bits);
1214:
1215: // Special rounding for maxval. If the number exceeds maxval by
1216: // any amount, even if it's less than half a step, it overflows.
1217: if (exp == 1023 && ((m >> 1) == (1L << 53) - 1)) {
1218: if (remainder || checkBits(il - ml))
1219: return neg ? Double.NEGATIVE_INFINITY
1220: : Double.POSITIVE_INFINITY;
1221: else
1222: return neg ? -Double.MAX_VALUE : Double.MAX_VALUE;
1223: }
1224:
1225: // Normal round-to-even rule: round up if the bit dropped is a one, and
1226: // the bit above it or any of the bits below it is a one.
1227: if ((m & 1) == 1
1228: && ((m & 2) == 2 || remainder || checkBits(excess_bits))) {
1229: m += 2;
1230: // Check if we overflowed the mantissa
1231: if ((m & (1L << 54)) != 0) {
1232: exp++;
1233: // renormalize
1234: m >>= 1;
1235: }
1236: // Check if a denormalized mantissa was just rounded up to a
1237: // normalized one.
1238: else if (ml == 52 && (m & (1L << 53)) != 0)
1239: exp++;
1240: }
1241:
1242: // Discard the rounding bit
1243: m >>= 1;
1244:
1245: long bits_sign = neg ? (1L << 63) : 0;
1246: exp += 1023;
1247: long bits_exp = (exp <= 0) ? 0 : ((long) exp) << 52;
1248: long bits_mant = m & ~(1L << 52);
1249: return Double
1250: .longBitsToDouble(bits_sign | bits_exp | bits_mant);
1251: }
1252:
1253: public Numeric add(Object y, int k) {
1254: if (y instanceof IntNum)
1255: return IntNum.add(this , (IntNum) y, k);
1256: if (!(y instanceof Numeric))
1257: throw new IllegalArgumentException();
1258: return ((Numeric) y).addReversed(this , k);
1259: }
1260:
1261: public Numeric mul(Object y) {
1262: if (y instanceof IntNum)
1263: return IntNum.times(this , (IntNum) y);
1264: if (!(y instanceof Numeric))
1265: throw new IllegalArgumentException();
1266: return ((Numeric) y).mulReversed(this );
1267: }
1268:
1269: public Numeric div(Object y) {
1270: if (y instanceof RatNum) {
1271: RatNum r = (RatNum) y;
1272: return RatNum.make(IntNum.times(this , r.denominator()), r
1273: .numerator());
1274: }
1275: if (!(y instanceof Numeric))
1276: throw new IllegalArgumentException();
1277: return ((Numeric) y).divReversed(this );
1278: }
1279:
1280: /** Copy the abolute value of this into an array of words.
1281: * Assumes words.length >= (this.words == null ? 1 : this.ival).
1282: * Result is zero-extended, but need not be a valid 2's complement number.
1283: */
1284:
1285: public void getAbsolute(int[] words) {
1286: int len;
1287: if (this .words == null) {
1288: len = 1;
1289: words[0] = this .ival;
1290: } else {
1291: len = this .ival;
1292: for (int i = len; --i >= 0;)
1293: words[i] = this .words[i];
1294: }
1295: if (words[len - 1] < 0)
1296: negate(words, words, len);
1297: for (int i = words.length; --i > len;)
1298: words[i] = 0;
1299: }
1300:
1301: /** Set dest[0:len-1] to the negation of src[0:len-1].
1302: * Return true if overflow (i.e. if src is -2**(32*len-1)).
1303: * Ok for src==dest. */
1304: public static boolean negate(int[] dest, int[] src, int len) {
1305: long carry = 1;
1306: boolean negative = src[len - 1] < 0;
1307: for (int i = 0; i < len; i++) {
1308: carry += ((long) (~src[i]) & 0xffffffffL);
1309: dest[i] = (int) carry;
1310: carry >>= 32;
1311: }
1312: return (negative && dest[len - 1] < 0);
1313: }
1314:
1315: /** Destructively set this to the negative of x.
1316: * It is OK if x==this.*/
1317: public void setNegative(IntNum x) {
1318: int len = x.ival;
1319: if (x.words == null) {
1320: if (len == Integer.MIN_VALUE)
1321: set(-(long) len);
1322: else
1323: set(-len);
1324: return;
1325: }
1326: realloc(len + 1);
1327: if (IntNum.negate(words, x.words, len))
1328: words[len++] = 0;
1329: ival = len;
1330: }
1331:
1332: /** Destructively negate this. */
1333: public final void setNegative() {
1334: setNegative(this );
1335: }
1336:
1337: public static IntNum abs(IntNum x) {
1338: return x.isNegative() ? neg(x) : x;
1339: }
1340:
1341: public static IntNum neg(IntNum x) {
1342: if (x.words == null && x.ival != Integer.MIN_VALUE)
1343: return make(-x.ival);
1344: IntNum result = new IntNum(0);
1345: result.setNegative(x);
1346: return result.canonicalize();
1347: }
1348:
1349: public Numeric neg() {
1350: return IntNum.neg(this );
1351: }
1352:
1353: /** Calculates {@code ceiling(log2(this < 0 ? -this : this+1))}.
1354: * See Common Lisp: the Language, 2nd ed, p. 361.
1355: */
1356: public int intLength() {
1357: if (words == null)
1358: return MPN.intLength(ival);
1359: else
1360: return MPN.intLength(words, ival);
1361: }
1362:
1363: /**
1364: * @serialData If the value is in the range (int)0xC000000 .. 0x7fffffff
1365: * (inclusive) write out the value (using writeInt).
1366: * Otherwise, write (using writeInt) (0x80000000|nwords), where nwords is
1367: * the number of words following. The words are the minimal
1368: * 2's complement big-endian representation of the value, written using
1369: * writeint.
1370: * (Even if the current value is not canonicalized, the output is).
1371: */
1372: public void writeExternal(ObjectOutput out) throws IOException {
1373: int nwords = words == null ? 1 : wordsNeeded(words, ival);
1374: if (nwords <= 1) {
1375: int i = words == null ? ival : words.length == 0 ? 0
1376: : words[0];
1377: if (i >= (int) 0xC0000000)
1378: out.writeInt(i);
1379: else {
1380: out.writeInt(0x80000001);
1381: out.writeInt(i);
1382: }
1383: } else {
1384: out.writeInt(0x80000000 | nwords);
1385: while (--nwords >= 0)
1386: out.writeInt(words[nwords]);
1387: }
1388:
1389: }
1390:
1391: public void readExternal(ObjectInput in) throws IOException,
1392: ClassNotFoundException {
1393: int i = in.readInt();
1394: if (ival <= (int) 0xC0000000) {
1395: i &= 0x7fffffff;
1396: if (i == 1)
1397: i = in.readInt();
1398: else {
1399: int[] w = new int[i];
1400: for (int j = i; --j >= 0;)
1401: w[j] = in.readInt();
1402: words = w;
1403: }
1404: }
1405: ival = i;
1406: }
1407:
1408: public Object readResolve() throws ObjectStreamException {
1409: return canonicalize();
1410: }
1411:
1412: public BigInteger asBigInteger() {
1413: if (words == null || ival <= 2)
1414: return BigInteger.valueOf(longValue());
1415: return new BigInteger(toString());
1416: }
1417:
1418: public BigDecimal asBigDecimal() {
1419: if (words == null)
1420: return new BigDecimal(ival);
1421: if (ival <= 2)
1422: return BigDecimal.valueOf(longValue());
1423: return new BigDecimal(toString());
1424: }
1425: }
|