0001: /*
0002: *
0003: *
0004: * Copyright 1990-2007 Sun Microsystems, Inc. All Rights Reserved.
0005: * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER
0006: *
0007: * This program is free software; you can redistribute it and/or
0008: * modify it under the terms of the GNU General Public License version
0009: * 2 only, as published by the Free Software Foundation.
0010: *
0011: * This program is distributed in the hope that it will be useful, but
0012: * WITHOUT ANY WARRANTY; without even the implied warranty of
0013: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
0014: * General Public License version 2 for more details (a copy is
0015: * included at /legal/license.txt).
0016: *
0017: * You should have received a copy of the GNU General Public License
0018: * version 2 along with this work; if not, write to the Free Software
0019: * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
0020: * 02110-1301 USA
0021: *
0022: * Please contact Sun Microsystems, Inc., 4150 Network Circle, Santa
0023: * Clara, CA 95054 or visit www.sun.com if you need additional
0024: * information or have any questions.
0025: */
0026:
0027: package java.lang;
0028:
0029: import com.sun.cldchi.jvm.JVM;
0030:
0031: class FloatingDecimal {
0032: boolean isExceptional;
0033: boolean isNegative;
0034: int decExponent;
0035: char digits[];
0036: int nDigits;
0037: int bigIntExp;
0038: int bigIntNBits;
0039: boolean mustSetRoundDir = false;
0040: int roundDir; // set by doubleValue
0041:
0042: private FloatingDecimal(boolean negSign, int decExponent,
0043: char[] digits, int n, boolean e) {
0044: isNegative = negSign;
0045: isExceptional = e;
0046: this .decExponent = decExponent;
0047: this .digits = digits;
0048: this .nDigits = n;
0049: }
0050:
0051: /*
0052: * Constants of the implementation
0053: * Most are IEEE-754 related.
0054: * (There are more really boring constants at the end.)
0055: */
0056: static final long signMask = 0x8000000000000000L;
0057: static final long expMask = 0x7ff0000000000000L;
0058: static final long fractMask = ~(signMask | expMask);
0059: static final int expShift = 52;
0060: static final int expBias = 1023;
0061: static final long fractHOB = (1L << expShift); // assumed High-Order bit
0062: static final long expOne = ((long) expBias) << expShift; // exponent of 1.0
0063: static final int maxSmallBinExp = 62;
0064: static final int minSmallBinExp = -(63 / 3);
0065: static final int maxDecimalDigits = 15;
0066: static final int maxDecimalExponent = 308;
0067: static final int minDecimalExponent = -324;
0068: static final int bigDecimalExponent = 324; // i.e. abs(minDecimalExponent)
0069:
0070: static final long highbyte = 0xff00000000000000L;
0071: static final long highbit = 0x8000000000000000L;
0072: static final long lowbytes = ~highbyte;
0073:
0074: static final int singleSignMask = 0x80000000;
0075: static final int singleExpMask = 0x7f800000;
0076: static final int singleFractMask = ~(singleSignMask | singleExpMask);
0077: static final int singleExpShift = 23;
0078: static final int singleFractHOB = 1 << singleExpShift;
0079: static final int singleExpBias = 127;
0080: static final int singleMaxDecimalDigits = 7;
0081: static final int singleMaxDecimalExponent = 38;
0082: static final int singleMinDecimalExponent = -45;
0083:
0084: static final int intDecimalDigits = 9;
0085:
0086: /*
0087: * count number of bits from high-order 1 bit to low-order 1 bit,
0088: * inclusive.
0089: */
0090: private static int countBits(long v) {
0091: //
0092: // the strategy is to shift until we get a non-zero sign bit
0093: // then shift until we have no bits left, counting the difference.
0094: // we do byte shifting as a special fix.
0095: //
0096: if (v == 0L)
0097: return 0;
0098:
0099: while ((v & highbyte) == 0L) {
0100: v <<= 8;
0101: }
0102: while (v > 0L) { // i.e. while ((v&highbit) == 0L )
0103: v <<= 1;
0104: }
0105:
0106: int n = 0;
0107: while ((v & lowbytes) != 0L) {
0108: v <<= 8;
0109: n += 8;
0110: }
0111: while (v != 0L) {
0112: v <<= 1;
0113: n += 1;
0114: }
0115: return n;
0116: }
0117:
0118: /*
0119: * Keep big powers of 5 handy for future reference.
0120: */
0121: private static FDBigInt b5p[];
0122:
0123: private static synchronized FDBigInt big5pow(int p) {
0124: if (p < 0) {
0125: throw new RuntimeException(
0126: /* #ifdef VERBOSE_EXCEPTIONS */
0127: /// skipped "Assertion botch: negative power of 5"
0128: /* #endif */
0129: );
0130: }
0131: if (b5p == null) {
0132: b5p = new FDBigInt[p + 1];
0133: } else if (b5p.length <= p) {
0134: FDBigInt t[] = new FDBigInt[p + 1];
0135: JVM.unchecked_obj_arraycopy(b5p, 0, t, 0, b5p.length);
0136: b5p = t;
0137: }
0138: if (b5p[p] != null)
0139: return b5p[p];
0140: else if (p < small5pow.length)
0141: return b5p[p] = new FDBigInt(small5pow[p]);
0142: else if (p < long5pow.length)
0143: return b5p[p] = new FDBigInt(long5pow[p]);
0144: else {
0145: // construct the value.
0146: // recursively.
0147: int q, r;
0148: // in order to compute 5^p,
0149: // compute its square root, 5^(p/2) and square.
0150: // or, let q = p / 2, r = p -q, then
0151: // 5^p = 5^(q+r) = 5^q * 5^r
0152: q = p >> 1;
0153: r = p - q;
0154: FDBigInt bigq = b5p[q];
0155: if (bigq == null)
0156: bigq = big5pow(q);
0157: if (r < small5pow.length) {
0158: return (b5p[p] = bigq.mult(small5pow[r]));
0159: } else {
0160: FDBigInt bigr = b5p[r];
0161: if (bigr == null)
0162: bigr = big5pow(r);
0163: return (b5p[p] = bigq.mult(bigr));
0164: }
0165: }
0166: }
0167:
0168: //
0169: // a common operation
0170: //
0171: private static FDBigInt multPow52(FDBigInt v, int p5, int p2) {
0172: if (p5 != 0) {
0173: if (p5 < small5pow.length) {
0174: v = v.mult(small5pow[p5]);
0175: } else {
0176: v = v.mult(big5pow(p5));
0177: }
0178: }
0179: if (p2 != 0) {
0180: v.lshiftMe(p2);
0181: }
0182: return v;
0183: }
0184:
0185: //
0186: // another common operation
0187: //
0188: private static FDBigInt constructPow52(int p5, int p2) {
0189: FDBigInt v = new FDBigInt(big5pow(p5));
0190: if (p2 != 0) {
0191: v.lshiftMe(p2);
0192: }
0193: return v;
0194: }
0195:
0196: /*
0197: * Make a floating double into a FDBigInt.
0198: * This could also be structured as a FDBigInt
0199: * constructor, but we'd have to build a lot of knowledge
0200: * about floating-point representation into it, and we don't want to.
0201: *
0202: * AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES
0203: * bigIntExp and bigIntNBits
0204: *
0205: */
0206: private FDBigInt doubleToBigInt(double dval) {
0207: long lbits = Double.doubleToLongBits(dval) & ~signMask;
0208: int binexp = (int) (lbits >>> expShift);
0209: lbits &= fractMask;
0210: if (binexp > 0) {
0211: lbits |= fractHOB;
0212: } else {
0213: if (lbits == 0L) {
0214: throw new RuntimeException(
0215: /* #ifdef VERBOSE_EXCEPTIONS */
0216: /// skipped "Assertion botch: doubleToBigInt(0.0)"
0217: /* #endif */
0218: );
0219: }
0220: binexp += 1;
0221: while ((lbits & fractHOB) == 0L) {
0222: lbits <<= 1;
0223: binexp -= 1;
0224: }
0225: }
0226: binexp -= expBias;
0227: int nbits = countBits(lbits);
0228: /*
0229: * We now know where the high-order 1 bit is,
0230: * and we know how many there are.
0231: */
0232: int lowOrderZeros = expShift + 1 - nbits;
0233: lbits >>>= lowOrderZeros;
0234:
0235: bigIntExp = binexp + 1 - nbits;
0236: bigIntNBits = nbits;
0237: return new FDBigInt(lbits);
0238: }
0239:
0240: /*
0241: * Compute a number that is the ULP of the given value,
0242: * for purposes of addition/subtraction. Generally easy.
0243: * More difficult if subtracting and the argument
0244: * is a normalized a power of 2, as the ULP changes at these points.
0245: */
0246: private static double ulp(double dval, boolean subtracting) {
0247: long lbits = Double.doubleToLongBits(dval) & ~signMask;
0248: int binexp = (int) (lbits >>> expShift);
0249: double ulpval;
0250: if (subtracting && (binexp >= expShift)
0251: && ((lbits & fractMask) == 0L)) {
0252: // for subtraction from normalized, powers of 2,
0253: // use next-smaller exponent
0254: binexp -= 1;
0255: }
0256: if (binexp > expShift) {
0257: ulpval = Double
0258: .longBitsToDouble(((long) (binexp - expShift)) << expShift);
0259: } else if (binexp == 0) {
0260: ulpval = Double.MIN_VALUE;
0261: } else {
0262: ulpval = Double.longBitsToDouble(1L << (binexp - 1));
0263: }
0264: if (subtracting)
0265: ulpval = -ulpval;
0266:
0267: return ulpval;
0268: }
0269:
0270: /*
0271: * Round a double to a float.
0272: * In addition to the fraction bits of the double,
0273: * look at the class instance variable roundDir,
0274: * which should help us avoid double-rounding error.
0275: * roundDir was set in hardValueOf if the estimate was
0276: * close enough, but not exact. It tells us which direction
0277: * of rounding is preferred.
0278: */
0279: float stickyRound(double dval) {
0280: long lbits = Double.doubleToLongBits(dval);
0281: long binexp = lbits & expMask;
0282: if (binexp == 0L || binexp == expMask) {
0283: // what we have here is special.
0284: // don't worry, the right thing will happen.
0285: return (float) dval;
0286: }
0287: lbits += (long) roundDir;
0288: return (float) Double.longBitsToDouble(lbits);
0289: }
0290:
0291: /*
0292: * This is the easy subcase --
0293: * all the significant bits, after scaling, are held in lvalue.
0294: * negSign and decExponent tell us what processing and scaling
0295: * has already been done. Exceptional cases have already been
0296: * stripped out.
0297: * In particular:
0298: * lvalue is a finite number (not Inf, nor NaN)
0299: * lvalue > 0L (not zero, nor negative).
0300: *
0301: * The only reason that we develop the digits here, rather than
0302: * calling on Long.toString() is that we can do it a little faster,
0303: * and besides want to treat trailing 0s specially. If Long.toString
0304: * changes, we should re-evaluate this strategy!
0305: */
0306: private void developLongDigits(int decExponent, long lvalue,
0307: long insignificant) {
0308: char digits[];
0309: int ndigits;
0310: int digitno;
0311: int c;
0312: //
0313: // Discard non-significant low-order bits, while rounding,
0314: // up to insignificant value.
0315: int i;
0316: for (i = 0; insignificant >= 10L; i++)
0317: insignificant /= 10L;
0318: if (i != 0) {
0319: long pow10 = long5pow[i] << i; // 10^i == 5^i * 2^i;
0320: long residue = lvalue % pow10;
0321: lvalue /= pow10;
0322: decExponent += i;
0323: if (residue >= (pow10 >> 1)) {
0324: // round up based on the low-order bits we're discarding
0325: lvalue++;
0326: }
0327: }
0328: if (lvalue <= Integer.MAX_VALUE) {
0329: if (lvalue <= 0L) {
0330: throw new RuntimeException(
0331: /* #ifdef VERBOSE_EXCEPTIONS */
0332: /// skipped "Assertion botch: value "+lvalue+" <= 0"
0333: /* #endif */
0334: );
0335: }
0336: // even easier subcase!
0337: // can do int arithmetic rather than long!
0338: int ivalue = (int) lvalue;
0339: digits = new char[ndigits = 10];
0340: digitno = ndigits - 1;
0341: c = ivalue % 10;
0342: ivalue /= 10;
0343: while (c == 0) {
0344: decExponent++;
0345: c = ivalue % 10;
0346: ivalue /= 10;
0347: }
0348: while (ivalue != 0) {
0349: digits[digitno--] = (char) (c + '0');
0350: decExponent++;
0351: c = ivalue % 10;
0352: ivalue /= 10;
0353: }
0354: digits[digitno] = (char) (c + '0');
0355: } else {
0356: // same algorithm as above (same bugs, too )
0357: // but using long arithmetic.
0358: digits = new char[ndigits = 20];
0359: digitno = ndigits - 1;
0360: c = (int) (lvalue % 10L);
0361: lvalue /= 10L;
0362: while (c == 0) {
0363: decExponent++;
0364: c = (int) (lvalue % 10L);
0365: lvalue /= 10L;
0366: }
0367: while (lvalue != 0L) {
0368: digits[digitno--] = (char) (c + '0');
0369: decExponent++;
0370: c = (int) (lvalue % 10L);
0371: lvalue /= 10;
0372: }
0373: digits[digitno] = (char) (c + '0');
0374: }
0375: char result[];
0376: ndigits -= digitno;
0377: if (digitno == 0)
0378: result = digits;
0379: else {
0380: result = new char[ndigits];
0381: JVM.unchecked_char_arraycopy(digits, digitno, result, 0,
0382: ndigits);
0383: }
0384: this .digits = result;
0385: this .decExponent = decExponent + 1;
0386: this .nDigits = ndigits;
0387: }
0388:
0389: //
0390: // add one to the least significant digit.
0391: // deal with carry out.
0392: // assert that this will only happen where there
0393: // is only one digit, e.g. (float)1e-44 seems to do it.
0394: //
0395: private void roundup() {
0396: int i;
0397: int q = digits[i = (nDigits - 1)];
0398: if (q == '9') {
0399: while (q == '9' && i > 0) {
0400: digits[i] = '0';
0401: q = digits[--i];
0402: }
0403: if (q == '9') {
0404: // carryout! High-order 1, rest 0s, larger exp.
0405: decExponent += 1;
0406: digits[0] = '1';
0407: return;
0408: }
0409: // else fall through.
0410: }
0411: digits[i] = (char) (q + 1);
0412: }
0413:
0414: /*
0415: * FIRST IMPORTANT CONSTRUCTOR: DOUBLE
0416: */
0417: public FloatingDecimal(double d) {
0418: long dBits = Double.doubleToLongBits(d);
0419: long fractBits;
0420: int binExp;
0421: int nSignificantBits;
0422:
0423: // discover and delete sign
0424: if ((dBits & signMask) != 0) {
0425: isNegative = true;
0426: dBits ^= signMask;
0427: } else {
0428: isNegative = false;
0429: }
0430: // Begin to unpack
0431: // Discover obvious special cases of NaN and Infinity.
0432: binExp = (int) ((dBits & expMask) >> expShift);
0433: fractBits = dBits & fractMask;
0434: if (binExp == (int) (expMask >> expShift)) {
0435: isExceptional = true;
0436: if (fractBits == 0L) {
0437: digits = infinity;
0438: } else {
0439: digits = notANumber;
0440: isNegative = false; // NaN has no sign!
0441: }
0442: nDigits = digits.length;
0443: return;
0444: }
0445: isExceptional = false;
0446: // Finish unpacking
0447: // Normalize denormalized numbers.
0448: // Insert assumed high-order bit for normalized numbers.
0449: // Subtract exponent bias.
0450: if (binExp == 0) {
0451: if (fractBits == 0L) {
0452: // not a denorm, just a 0!
0453: decExponent = 0;
0454: digits = zero;
0455: nDigits = 1;
0456: return;
0457: }
0458: while ((fractBits & fractHOB) == 0L) {
0459: fractBits <<= 1;
0460: binExp -= 1;
0461: }
0462: nSignificantBits = expShift + binExp + 1; // recall binExp is - shift count.
0463: binExp += 1;
0464: } else {
0465: fractBits |= fractHOB;
0466: nSignificantBits = expShift + 1;
0467: }
0468: binExp -= expBias;
0469: // call the routine that actually does all the hard work.
0470: dtoa(binExp, fractBits, nSignificantBits);
0471: }
0472:
0473: /*
0474: * SECOND IMPORTANT CONSTRUCTOR: SINGLE
0475: */
0476: /* public */FloatingDecimal(float f) {
0477: int fBits = Float.floatToIntBits(f);
0478: int fractBits;
0479: int binExp;
0480: int nSignificantBits;
0481:
0482: // discover and delete sign
0483: if ((fBits & singleSignMask) != 0) {
0484: isNegative = true;
0485: fBits ^= singleSignMask;
0486: } else {
0487: isNegative = false;
0488: }
0489: // Begin to unpack
0490: // Discover obvious special cases of NaN and Infinity.
0491: binExp = (int) ((fBits & singleExpMask) >> singleExpShift);
0492: fractBits = fBits & singleFractMask;
0493: if (binExp == (int) (singleExpMask >> singleExpShift)) {
0494: isExceptional = true;
0495: if (fractBits == 0L) {
0496: digits = infinity;
0497: } else {
0498: digits = notANumber;
0499: isNegative = false; // NaN has no sign!
0500: }
0501: nDigits = digits.length;
0502: return;
0503: }
0504: isExceptional = false;
0505: // Finish unpacking
0506: // Normalize denormalized numbers.
0507: // Insert assumed high-order bit for normalized numbers.
0508: // Subtract exponent bias.
0509: if (binExp == 0) {
0510: if (fractBits == 0) {
0511: // not a denorm, just a 0!
0512: decExponent = 0;
0513: digits = zero;
0514: nDigits = 1;
0515: return;
0516: }
0517: while ((fractBits & singleFractHOB) == 0) {
0518: fractBits <<= 1;
0519: binExp -= 1;
0520: }
0521: nSignificantBits = singleExpShift + binExp + 1; // recall binExp is - shift count.
0522: binExp += 1;
0523: } else {
0524: fractBits |= singleFractHOB;
0525: nSignificantBits = singleExpShift + 1;
0526: }
0527: binExp -= singleExpBias;
0528: // call the routine that actually does all the hard work.
0529: dtoa(binExp, ((long) fractBits) << (expShift - singleExpShift),
0530: nSignificantBits);
0531: }
0532:
0533: private void dtoa(int binExp, long fractBits, int nSignificantBits) {
0534: int nFractBits; // number of significant bits of fractBits;
0535: int nTinyBits; // number of these to the right of the point.
0536: int decExp;
0537:
0538: // Examine number. Determine if it is an easy case,
0539: // which we can do pretty trivially using float/long conversion,
0540: // or whether we must do real work.
0541: nFractBits = countBits(fractBits);
0542: nTinyBits = Math.max(0, nFractBits - binExp - 1);
0543: if (binExp <= maxSmallBinExp && binExp >= minSmallBinExp) {
0544: // Look more closely at the number to decide if,
0545: // with scaling by 10^nTinyBits, the result will fit in
0546: // a long.
0547: if ((nTinyBits < long5pow.length)
0548: && ((nFractBits + n5bits[nTinyBits]) < 64)) {
0549: /*
0550: * We can do this:
0551: * take the fraction bits, which are normalized.
0552: * (a) nTinyBits == 0: Shift left or right appropriately
0553: * to align the binary point at the extreme right, i.e.
0554: * where a long int point is expected to be. The integer
0555: * result is easily converted to a string.
0556: * (b) nTinyBits > 0: Shift right by expShift-nFractBits,
0557: * which effectively converts to long and scales by
0558: * 2^nTinyBits. Then multiply by 5^nTinyBits to
0559: * complete the scaling. We know this won't overflow
0560: * because we just counted the number of bits necessary
0561: * in the result. The integer you get from this can
0562: * then be converted to a string pretty easily.
0563: */
0564: long halfULP;
0565: if (nTinyBits == 0) {
0566: if (binExp > nSignificantBits) {
0567: halfULP = 1L << (binExp - nSignificantBits - 1);
0568: } else {
0569: halfULP = 0L;
0570: }
0571: if (binExp >= expShift) {
0572: fractBits <<= (binExp - expShift);
0573: } else {
0574: fractBits >>>= (expShift - binExp);
0575: }
0576: developLongDigits(0, fractBits, halfULP);
0577: return;
0578: }
0579: /*
0580: * The following causes excess digits to be printed
0581: * out in the single-float case. Our manipulation of
0582: * halfULP here is apparently not correct. If we
0583: * better understand how this works, perhaps we can
0584: * use this special case again. But for the time being,
0585: * we do not.
0586: * else {
0587: * fractBits >>>= expShift+1-nFractBits;
0588: * fractBits *= long5pow[ nTinyBits ];
0589: * halfULP = long5pow[ nTinyBits ] >> (1+nSignificantBits-nFractBits);
0590: * developLongDigits( -nTinyBits, fractBits, halfULP );
0591: * return;
0592: * }
0593: */
0594: }
0595: }
0596: /*
0597: * This is the hard case. We are going to compute large positive
0598: * integers B and S and integer decExp, s.t.
0599: * d = ( B / S ) * 10^decExp
0600: * 1 <= B / S < 10
0601: * Obvious choices are:
0602: * decExp = floor( log10(d) )
0603: * B = d * 2^nTinyBits * 10^max( 0, -decExp )
0604: * S = 10^max( 0, decExp) * 2^nTinyBits
0605: * (noting that nTinyBits has already been forced to non-negative)
0606: * I am also going to compute a large positive integer
0607: * M = (1/2^nSignificantBits) * 2^nTinyBits * 10^max( 0, -decExp )
0608: * i.e. M is (1/2) of the ULP of d, scaled like B.
0609: * When we iterate through dividing B/S and picking off the
0610: * quotient bits, we will know when to stop when the remainder
0611: * is <= M.
0612: *
0613: * We keep track of powers of 2 and powers of 5.
0614: */
0615:
0616: /*
0617: * Estimate decimal exponent. (If it is small-ish,
0618: * we could double-check.)
0619: *
0620: * First, scale the mantissa bits such that 1 <= d2 < 2.
0621: * We are then going to estimate
0622: * log10(d2) ~=~ (d2-1.5)/1.5 + log(1.5)
0623: * and so we can estimate
0624: * log10(d) ~=~ log10(d2) + binExp * log10(2)
0625: * take the floor and call it decExp.
0626: * IMPL_NOTE -- use more precise constants here. It costs no more.
0627: */
0628: double d2 = Double.longBitsToDouble(expOne
0629: | (fractBits & ~fractHOB));
0630: decExp = (int) Math.floor((d2 - 1.5D) * 0.289529654D
0631: + 0.176091259 + (double) binExp * 0.301029995663981);
0632: int B2, B5; // powers of 2 and powers of 5, respectively, in B
0633: int S2, S5; // powers of 2 and powers of 5, respectively, in S
0634: int M2, M5; // powers of 2 and powers of 5, respectively, in M
0635: int Bbits; // binary digits needed to represent B, approx.
0636: int tenSbits; // binary digits needed to represent 10*S, approx.
0637: FDBigInt Sval, Bval, Mval;
0638:
0639: B5 = Math.max(0, -decExp);
0640: B2 = B5 + nTinyBits + binExp;
0641:
0642: S5 = Math.max(0, decExp);
0643: S2 = S5 + nTinyBits;
0644:
0645: M5 = B5;
0646: M2 = B2 - nSignificantBits;
0647:
0648: /*
0649: * the long integer fractBits contains the (nFractBits) interesting
0650: * bits from the mantissa of d ( hidden 1 added if necessary) followed
0651: * by (expShift+1-nFractBits) zeros. In the interest of compactness,
0652: * I will shift out those zeros before turning fractBits into a
0653: * FDBigInt. The resulting whole number will be
0654: * d * 2^(nFractBits-1-binExp).
0655: */
0656: fractBits >>>= (expShift + 1 - nFractBits);
0657: B2 -= nFractBits - 1;
0658: int common2factor = Math.min(B2, S2);
0659: B2 -= common2factor;
0660: S2 -= common2factor;
0661: M2 -= common2factor;
0662:
0663: /*
0664: * IMPL_NOTE: For exact powers of two, the next smallest number
0665: * is only half as far away as we think (because the meaning of
0666: * ULP changes at power-of-two bounds) for this reason, we
0667: * fix M2.
0668: */
0669: if (nFractBits == 1)
0670: M2 -= 1;
0671:
0672: if (M2 < 0) {
0673: // oops.
0674: // since we cannot scale M down far enough,
0675: // we must scale the other values up.
0676: B2 -= M2;
0677: S2 -= M2;
0678: M2 = 0;
0679: }
0680: /*
0681: * Construct, Scale, iterate.
0682: * We use a symmetric test.
0683: */
0684: char digits[] = this .digits = new char[18];
0685: int ndigit = 0;
0686: boolean low, high;
0687: long lowDigitDifference;
0688: int q;
0689:
0690: /*
0691: * Detect the special cases where all the numbers we are about
0692: * to compute will fit in int or long integers.
0693: * In these cases, we will avoid doing FDBigInt arithmetic.
0694: * We use the same algorithms, except that we "normalize"
0695: * our FDBigInts before iterating. This is to make division easier,
0696: * as it makes our fist guess (quotient of high-order words)
0697: * more accurate!
0698: *
0699: * We use a symmetric test.
0700: */
0701: Bbits = nFractBits + B2
0702: + ((B5 < n5bits.length) ? n5bits[B5] : (B5 * 3));
0703: tenSbits = S2
0704: + 1
0705: + (((S5 + 1) < n5bits.length) ? n5bits[(S5 + 1)]
0706: : ((S5 + 1) * 3));
0707: if (Bbits < 64 && tenSbits < 64) {
0708: if (Bbits < 32 && tenSbits < 32) {
0709: // wa-hoo! They're all ints!
0710: int b = ((int) fractBits * small5pow[B5]) << B2;
0711: int s = small5pow[S5] << S2;
0712: int m = small5pow[M5] << M2;
0713: int tens = s * 10;
0714: /*
0715: * Unroll the first iteration. If our decExp estimate
0716: * was too high, our first quotient will be zero. In this
0717: * case, we discard it and decrement decExp.
0718: */
0719: ndigit = 0;
0720: q = (int) (b / s);
0721: b = 10 * (b % s);
0722: m *= 10;
0723: low = (b < m);
0724: high = (b + m > tens);
0725: if (q >= 10) {
0726: // bummer, dude
0727: throw new RuntimeException(
0728: /* #ifdef VERBOSE_EXCEPTIONS */
0729: /// skipped "Assertion botch: excessivly large digit "+q
0730: /* #endif */
0731: );
0732: } else if ((q == 0) && !high) {
0733: // oops. Usually ignore leading zero.
0734: decExp--;
0735: } else {
0736: digits[ndigit++] = (char) ('0' + q);
0737: }
0738: /*
0739: * IMPL_NOTE: Java spec sez that we always have at least
0740: * one digit after the . in either F- or E-form output.
0741: * Thus we will need more than one digit if we're using
0742: * E-form
0743: */
0744: if (decExp <= -3 || decExp >= 8) {
0745: high = low = false;
0746: }
0747: while (!low && !high) {
0748: q = (int) (b / s);
0749: b = 10 * (b % s);
0750: m *= 10;
0751: if (q >= 10) {
0752: // bummer, dude
0753: throw new RuntimeException(
0754: /* #ifdef VERBOSE_EXCEPTIONS */
0755: /// skipped "Assertion botch: excessivly large digit "+q
0756: /* #endif */
0757: );
0758: }
0759: if (m > 0L) {
0760: low = (b < m);
0761: high = (b + m > tens);
0762: } else {
0763: // m might overflow!
0764: // in this case, it is certainly > b,
0765: // which won't
0766: // and b+m > tens, too, since that has overflowed
0767: // either!
0768: low = true;
0769: high = true;
0770: }
0771: digits[ndigit++] = (char) ('0' + q);
0772: }
0773: lowDigitDifference = (b << 1) - tens;
0774: } else {
0775: // still good! they're all longs!
0776: long b = (fractBits * long5pow[B5]) << B2;
0777: long s = long5pow[S5] << S2;
0778: long m = long5pow[M5] << M2;
0779: long tens = s * 10L;
0780: /*
0781: * Unroll the first iteration. If our decExp estimate
0782: * was too high, our first quotient will be zero. In this
0783: * case, we discard it and decrement decExp.
0784: */
0785: ndigit = 0;
0786: q = (int) (b / s);
0787: b = 10L * (b % s);
0788: m *= 10L;
0789: low = (b < m);
0790: high = (b + m > tens);
0791: if (q >= 10) {
0792: // bummer, dude
0793: throw new RuntimeException(
0794: /* #ifdef VERBOSE_EXCEPTIONS */
0795: /// skipped "Assertion botch: excessivly large digit "+q
0796: /* #endif */
0797: );
0798: } else if ((q == 0) && !high) {
0799: // oops. Usually ignore leading zero.
0800: decExp--;
0801: } else {
0802: digits[ndigit++] = (char) ('0' + q);
0803: }
0804: /*
0805: * IMPL_NOTE: Java spec sez that we always have at least
0806: * one digit after the . in either F- or E-form output.
0807: * Thus we will need more than one digit if we're using
0808: * E-form
0809: */
0810: if (decExp <= -3 || decExp >= 8) {
0811: high = low = false;
0812: }
0813: while (!low && !high) {
0814: q = (int) (b / s);
0815: b = 10 * (b % s);
0816: m *= 10;
0817: if (q >= 10) {
0818: // bummer, dude
0819: throw new RuntimeException(
0820: /* #ifdef VERBOSE_EXCEPTIONS */
0821: /// skipped "Assertion botch: excessivly large digit "+q
0822: /* #endif */
0823: );
0824: }
0825: if (m > 0L) {
0826: low = (b < m);
0827: high = (b + m > tens);
0828: } else {
0829: // m might overflow!
0830: // in this case, it is certainly > b,
0831: // which won't
0832: // and b+m > tens, too, since that has overflowed
0833: // either!
0834: low = true;
0835: high = true;
0836: }
0837: digits[ndigit++] = (char) ('0' + q);
0838: }
0839: lowDigitDifference = (b << 1) - tens;
0840: }
0841: } else {
0842: FDBigInt tenSval;
0843: int shiftBias;
0844:
0845: /*
0846: * We really must do FDBigInt arithmetic.
0847: * Fist, construct our FDBigInt initial values.
0848: */
0849: Bval = multPow52(new FDBigInt(fractBits), B5, B2);
0850: Sval = constructPow52(S5, S2);
0851: Mval = constructPow52(M5, M2);
0852:
0853: // normalize so that division works better
0854: Bval.lshiftMe(shiftBias = Sval.normalizeMe());
0855: Mval.lshiftMe(shiftBias);
0856: tenSval = Sval.mult(10);
0857: /*
0858: * Unroll the first iteration. If our decExp estimate
0859: * was too high, our first quotient will be zero. In this
0860: * case, we discard it and decrement decExp.
0861: */
0862: ndigit = 0;
0863: q = Bval.quoRemIteration(Sval);
0864: Mval = Mval.mult(10);
0865: low = (Bval.cmp(Mval) < 0);
0866: high = (Bval.add(Mval).cmp(tenSval) > 0);
0867: if (q >= 10) {
0868: // bummer, dude
0869: throw new RuntimeException(
0870: /* #ifdef VERBOSE_EXCEPTIONS */
0871: /// skipped "Assertion botch: excessivly large digit "+q
0872: /* #endif */
0873: );
0874: } else if ((q == 0) && !high) {
0875: // oops. Usually ignore leading zero.
0876: decExp--;
0877: } else {
0878: digits[ndigit++] = (char) ('0' + q);
0879: }
0880: /*
0881: * IMPL_NOTE: Java spec sez that we always have at least
0882: * one digit after the . in either F- or E-form output.
0883: * Thus we will need more than one digit if we're using
0884: * E-form
0885: */
0886: if (decExp <= -3 || decExp >= 8) {
0887: high = low = false;
0888: }
0889: while (!low && !high) {
0890: q = Bval.quoRemIteration(Sval);
0891: Mval = Mval.mult(10);
0892: if (q >= 10) {
0893: // bummer, dude
0894: throw new RuntimeException(
0895: /* #ifdef VERBOSE_EXCEPTIONS */
0896: /// skipped "Assertion botch: excessivly large digit "+q
0897: /* #endif */
0898: );
0899: }
0900: low = (Bval.cmp(Mval) < 0);
0901: high = (Bval.add(Mval).cmp(tenSval) > 0);
0902: digits[ndigit++] = (char) ('0' + q);
0903: }
0904: if (high && low) {
0905: Bval.lshiftMe(1);
0906: lowDigitDifference = Bval.cmp(tenSval);
0907: } else
0908: lowDigitDifference = 0L; // this here only for flow analysis!
0909: }
0910: this .decExponent = decExp + 1;
0911: this .digits = digits;
0912: this .nDigits = ndigit;
0913: /*
0914: * Last digit gets rounded based on stopping condition.
0915: */
0916: if (high) {
0917: if (low) {
0918: if (lowDigitDifference == 0L) {
0919: // it's a tie!
0920: // choose based on which digits we like.
0921: if ((digits[nDigits - 1] & 1) != 0)
0922: roundup();
0923: } else if (lowDigitDifference > 0) {
0924: roundup();
0925: }
0926: } else {
0927: roundup();
0928: }
0929: }
0930: }
0931:
0932: public String toString() {
0933: StringBuffer result = new StringBuffer(nDigits + 8);
0934: if (isNegative) {
0935: result.append('-');
0936: }
0937: if (isExceptional) {
0938: result.append(digits, 0, nDigits);
0939: } else {
0940: result.append("0.");
0941: result.append(digits, 0, nDigits);
0942: result.append('e');
0943: result.append(decExponent);
0944: }
0945: return new String(result);
0946: }
0947:
0948: public String toJavaFormatString() {
0949: char result[] = new char[nDigits + 10];
0950: int i = 0;
0951: if (isNegative) {
0952: result[0] = '-';
0953: i = 1;
0954: }
0955: if (isExceptional) {
0956: JVM.unchecked_char_arraycopy(digits, 0, result, i, nDigits);
0957: i += nDigits;
0958: } else {
0959: if (decExponent > 0 && decExponent < 8) {
0960: // print digits.digits.
0961: int charLength = Math.min(nDigits, decExponent);
0962: JVM.unchecked_char_arraycopy(digits, 0, result, i,
0963: charLength);
0964: i += charLength;
0965: if (charLength < decExponent) {
0966: charLength = decExponent - charLength;
0967: JVM.unchecked_char_arraycopy(zero, 0, result, i,
0968: charLength);
0969: i += charLength;
0970: result[i++] = '.';
0971: result[i++] = '0';
0972: } else {
0973: result[i++] = '.';
0974: if (charLength < nDigits) {
0975: int t = nDigits - charLength;
0976: JVM.unchecked_char_arraycopy(digits,
0977: charLength, result, i, t);
0978: i += t;
0979: } else {
0980: result[i++] = '0';
0981: }
0982: }
0983: } else if (decExponent <= 0 && decExponent > -3) {
0984: result[i++] = '0';
0985: result[i++] = '.';
0986: if (decExponent != 0) {
0987: JVM.unchecked_char_arraycopy(zero, 0, result, i,
0988: -decExponent);
0989: i -= decExponent;
0990: }
0991: JVM.unchecked_char_arraycopy(digits, 0, result, i,
0992: nDigits);
0993: i += nDigits;
0994: } else {
0995: result[i++] = digits[0];
0996: result[i++] = '.';
0997: if (nDigits > 1) {
0998: JVM.unchecked_char_arraycopy(digits, 1, result, i,
0999: nDigits - 1);
1000: i += nDigits - 1;
1001: } else {
1002: result[i++] = '0';
1003: }
1004: result[i++] = 'E';
1005: int e;
1006: if (decExponent <= 0) {
1007: result[i++] = '-';
1008: e = -decExponent + 1;
1009: } else {
1010: e = decExponent - 1;
1011: }
1012: // decExponent has 1, 2, or 3, digits
1013: if (e <= 9) {
1014: result[i++] = (char) (e + '0');
1015: } else if (e <= 99) {
1016: result[i++] = (char) (e / 10 + '0');
1017: result[i++] = (char) (e % 10 + '0');
1018: } else {
1019: result[i++] = (char) (e / 100 + '0');
1020: e %= 100;
1021: result[i++] = (char) (e / 10 + '0');
1022: result[i++] = (char) (e % 10 + '0');
1023: }
1024: }
1025: }
1026: return new String(result, 0, i);
1027: }
1028:
1029: public static FloatingDecimal readJavaFormatString(String in)
1030: throws NumberFormatException {
1031: boolean isNegative = false;
1032: boolean signSeen = false;
1033: int decExp;
1034: char c;
1035:
1036: parseNumber: try {
1037: in = in.trim(); // throws NullPointerException if null
1038: int l = in.length();
1039: if (l == 0) {
1040: throw new NumberFormatException(
1041: /* #ifdef VERBOSE_EXCEPTIONS */
1042: /// skipped "empty String"
1043: /* #endif */
1044: );
1045: }
1046: int i = 0;
1047: switch (c = in.charAt(i)) {
1048: case '-':
1049: isNegative = true;
1050: //FALLTHROUGH
1051: case '+':
1052: i++;
1053: signSeen = true;
1054: }
1055: // Would handle NaN and Infinity here, but it isn't
1056: // part of the spec!
1057: //
1058: char[] digits = new char[l];
1059: int nDigits = 0;
1060: boolean decSeen = false;
1061: int decPt = 0;
1062: int nLeadZero = 0;
1063: int nTrailZero = 0;
1064: digitLoop: while (i < l) {
1065: switch (c = in.charAt(i)) {
1066: case '0':
1067: if (nDigits > 0) {
1068: nTrailZero += 1;
1069: } else {
1070: nLeadZero += 1;
1071: }
1072: break; // out of switch.
1073: case '1':
1074: case '2':
1075: case '3':
1076: case '4':
1077: case '5':
1078: case '6':
1079: case '7':
1080: case '8':
1081: case '9':
1082: while (nTrailZero > 0) {
1083: digits[nDigits++] = '0';
1084: nTrailZero -= 1;
1085: }
1086: digits[nDigits++] = c;
1087: break; // out of switch.
1088: case '.':
1089: if (decSeen) {
1090: // already saw one ., this is the 2nd.
1091: throw new NumberFormatException(
1092: /* #ifdef VERBOSE_EXCEPTIONS */
1093: /// skipped "multiple points"
1094: /* #endif */
1095: );
1096: }
1097: decPt = i;
1098: if (signSeen) {
1099: decPt -= 1;
1100: }
1101: decSeen = true;
1102: break; // out of switch.
1103: default:
1104: break digitLoop;
1105: }
1106: i++;
1107: }
1108: /*
1109: * At this point, we've scanned all the digits and decimal
1110: * point we're going to see. Trim off leading and trailing
1111: * zeros, which will just confuse us later, and adjust
1112: * our initial decimal exponent accordingly.
1113: * To review:
1114: * we have seen i total characters.
1115: * nLeadZero of them were zeros before any other digits.
1116: * nTrailZero of them were zeros after any other digits.
1117: * if ( decSeen ), then a . was seen after decPt characters
1118: * ( including leading zeros which have been discarded )
1119: * nDigits characters were neither lead nor trailing
1120: * zeros, nor point
1121: */
1122: /*
1123: * special fix: if we saw no non-zero digits, then the
1124: * answer is zero!
1125: * Unfortunately, we feel honor-bound to keep parsing!
1126: */
1127: if (nDigits == 0) {
1128: digits = zero;
1129: nDigits = 1;
1130: if (nLeadZero == 0) {
1131: // we saw NO DIGITS AT ALL,
1132: // not even a crummy 0!
1133: // this is not allowed.
1134: break parseNumber; // go throw exception
1135: }
1136:
1137: }
1138:
1139: /* Our initial exponent is decPt, adjusted by the number of
1140: * discarded zeros. Or, if there was no decPt,
1141: * then its just nDigits adjusted by discarded trailing zeros.
1142: */
1143: if (decSeen) {
1144: decExp = decPt - nLeadZero;
1145: } else {
1146: decExp = nDigits + nTrailZero;
1147: }
1148:
1149: /*
1150: * Look for 'e' or 'E' and an optionally signed integer.
1151: */
1152: if ((i < l) && ((c = in.charAt(i)) == 'e') || (c == 'E')) {
1153: int expSign = 1;
1154: int expVal = 0;
1155: int reallyBig = Integer.MAX_VALUE / 10;
1156: boolean expOverflow = false;
1157: switch (in.charAt(++i)) {
1158: case '-':
1159: expSign = -1;
1160: //FALLTHROUGH
1161: case '+':
1162: i++;
1163: }
1164: int expAt = i;
1165: expLoop: while (i < l) {
1166: if (expVal >= reallyBig) {
1167: // the next character will cause integer
1168: // overflow.
1169: expOverflow = true;
1170: }
1171: switch (c = in.charAt(i++)) {
1172: case '0':
1173: case '1':
1174: case '2':
1175: case '3':
1176: case '4':
1177: case '5':
1178: case '6':
1179: case '7':
1180: case '8':
1181: case '9':
1182: expVal = expVal * 10 + ((int) c - (int) '0');
1183: continue;
1184: default:
1185: i--; // back up.
1186: break expLoop; // stop parsing exponent.
1187: }
1188: }
1189: int expLimit = bigDecimalExponent + nDigits
1190: + nTrailZero;
1191: if (expOverflow || (expVal > expLimit)) {
1192: //
1193: // The intent here is to end up with
1194: // infinity or zero, as appropriate.
1195: // The reason for yielding such a small decExponent,
1196: // rather than something intuitive such as
1197: // expSign*Integer.MAX_VALUE, is that this value
1198: // is subject to further manipulation in
1199: // doubleValue() and floatValue(), and I don't want
1200: // it to be able to cause overflow there!
1201: // (The only way we can get into trouble here is for
1202: // really outrageous nDigits+nTrailZero, such as 2 billion. )
1203: //
1204: decExp = expSign * expLimit;
1205: } else {
1206: // this should not overflow, since we tested
1207: // for expVal > (MAX+N), where N >= abs(decExp)
1208: decExp = decExp + expSign * expVal;
1209: }
1210:
1211: // if we saw something not a digit ( or end of string )
1212: // after the [Ee][+-], without seeing any digits at all
1213: // this is certainly an error. If we saw some digits,
1214: // but then some trailing garbage, that might be ok.
1215: // so we just fall through in that case.
1216: // HUMBUG
1217: if (i == expAt)
1218: break parseNumber; // certainly bad
1219: }
1220: /*
1221: * We parsed everything we could.
1222: * If there are leftovers, then this is not good input!
1223: */
1224: if (i < l
1225: && ((i != l - 1) || (in.charAt(i) != 'f'
1226: && in.charAt(i) != 'F'
1227: && in.charAt(i) != 'd' && in.charAt(i) != 'D'))) {
1228: break parseNumber; // go throw exception
1229: }
1230:
1231: return new FloatingDecimal(isNegative, decExp, digits,
1232: nDigits, false);
1233: } catch (StringIndexOutOfBoundsException e) {
1234: }
1235:
1236: throw new NumberFormatException(
1237: /* #ifdef VERBOSE_EXCEPTIONS */
1238: /// skipped in
1239: /* #endif */
1240: );
1241: }
1242:
1243: /*
1244: * Take a FloatingDecimal, which we presumably just scanned in,
1245: * and find out what its value is, as a double.
1246: *
1247: * AS A SIDE EFFECT, SET roundDir TO INDICATE PREFERRED
1248: * ROUNDING DIRECTION in case the result is really destined
1249: * for a single-precision float.
1250: */
1251:
1252: public double doubleValue() {
1253: int kDigits = Math.min(nDigits, maxDecimalDigits + 1);
1254: long lValue;
1255: double dValue;
1256: double rValue, tValue;
1257:
1258: roundDir = 0;
1259:
1260: /*
1261: * convert the lead kDigits to a long integer.
1262: */
1263: // (special performance fix: start to do it using int)
1264: int iValue = (int) digits[0] - (int) '0';
1265: int iDigits = Math.min(kDigits, intDecimalDigits);
1266: for (int i = 1; i < iDigits; i++) {
1267: iValue = iValue * 10 + (int) digits[i] - (int) '0';
1268: }
1269: lValue = (long) iValue;
1270: for (int i = iDigits; i < kDigits; i++) {
1271: lValue = lValue * 10L
1272: + (long) ((int) digits[i] - (int) '0');
1273: }
1274: dValue = (double) lValue;
1275: int exp = decExponent - kDigits;
1276: /*
1277: * lValue now contains a long integer with the value of
1278: * the first kDigits digits of the number.
1279: * dValue contains the (double) of the same.
1280: */
1281:
1282: if (nDigits <= maxDecimalDigits) {
1283: /*
1284: * possibly an easy case.
1285: * We know that the digits can be represented
1286: * exactly. And if the exponent isn't too outrageous,
1287: * the whole thing can be done with one operation,
1288: * thus one rounding error.
1289: * Note that all our constructors trim all leading and
1290: * trailing zeros, so simple values (including zero)
1291: * will always end up here
1292: */
1293: if (exp == 0 || dValue == 0.0)
1294: return (isNegative) ? -dValue : dValue; // small floating integer
1295: else if (exp >= 0) {
1296: if (exp <= maxSmallTen) {
1297: /*
1298: * Can get the answer with one operation,
1299: * thus one roundoff.
1300: */
1301: rValue = dValue * small10pow[exp];
1302: if (mustSetRoundDir) {
1303: tValue = rValue / small10pow[exp];
1304: roundDir = (tValue == dValue) ? 0
1305: : (tValue < dValue) ? 1 : -1;
1306: }
1307: return (isNegative) ? -rValue : rValue;
1308: }
1309: int slop = maxDecimalDigits - kDigits;
1310: if (exp <= maxSmallTen + slop) {
1311: /*
1312: * We can multiply dValue by 10^(slop)
1313: * and it is still "small" and exact.
1314: * Then we can multiply by 10^(exp-slop)
1315: * with one rounding.
1316: */
1317: dValue *= small10pow[slop];
1318: rValue = dValue * small10pow[exp - slop];
1319:
1320: if (mustSetRoundDir) {
1321: tValue = rValue / small10pow[exp - slop];
1322: roundDir = (tValue == dValue) ? 0
1323: : (tValue < dValue) ? 1 : -1;
1324: }
1325: return (isNegative) ? -rValue : rValue;
1326: }
1327: /*
1328: * Else we have a hard case with a positive exp.
1329: */
1330: } else {
1331: if (exp >= -maxSmallTen) {
1332: /*
1333: * Can get the answer in one division.
1334: */
1335: rValue = dValue / small10pow[-exp];
1336: tValue = rValue * small10pow[-exp];
1337: if (mustSetRoundDir) {
1338: roundDir = (tValue == dValue) ? 0
1339: : (tValue < dValue) ? 1 : -1;
1340: }
1341: return (isNegative) ? -rValue : rValue;
1342: }
1343: /*
1344: * Else we have a hard case with a negative exp.
1345: */
1346: }
1347: }
1348:
1349: /*
1350: * Harder cases:
1351: * The sum of digits plus exponent is greater than
1352: * what we think we can do with one error.
1353: *
1354: * Start by approximating the right answer by,
1355: * naively, scaling by powers of 10.
1356: */
1357: if (exp > 0) {
1358: if (decExponent > maxDecimalExponent + 1) {
1359: /*
1360: * Lets face it. This is going to be
1361: * Infinity. Cut to the chase.
1362: */
1363: return (isNegative) ? Double.NEGATIVE_INFINITY
1364: : Double.POSITIVE_INFINITY;
1365: }
1366: if ((exp & 15) != 0) {
1367: dValue *= small10pow[exp & 15];
1368: }
1369: if ((exp >>= 4) != 0) {
1370: int j = 0;
1371: while (exp > 1) {
1372: if ((exp & 1) != 0)
1373: dValue *= big10pow[j];
1374: j++;
1375: exp >>= 1;
1376: }
1377: /*
1378: * The reason for the exp > 1 condition
1379: * in the above loop was so that the last multiply
1380: * would get unrolled. We handle it here.
1381: * It could overflow.
1382: */
1383: double t = dValue * big10pow[j];
1384: if (Double.isInfinite(t)) {
1385: /*
1386: * It did overflow.
1387: * Look more closely at the result.
1388: * If the exponent is just one too large,
1389: * then use the maximum finite as our estimate
1390: * value. Else call the result infinity
1391: * and punt it.
1392: * ( I presume this could happen because
1393: * rounding forces the result here to be
1394: * an ULP or two larger than
1395: * Double.MAX_VALUE ).
1396: */
1397: t = dValue / 2.0;
1398: t *= big10pow[j];
1399: if (Double.isInfinite(t)) {
1400: return (isNegative) ? Double.NEGATIVE_INFINITY
1401: : Double.POSITIVE_INFINITY;
1402: }
1403: t = Double.MAX_VALUE;
1404: }
1405: dValue = t;
1406: }
1407: } else if (exp < 0) {
1408: exp = -exp;
1409: if (decExponent < minDecimalExponent - 1) {
1410: /*
1411: * Lets face it. This is going to be
1412: * zero. Cut to the chase.
1413: */
1414: return (isNegative) ? -0.0 : 0.0;
1415: }
1416: if ((exp & 15) != 0) {
1417: dValue /= small10pow[exp & 15];
1418: }
1419: if ((exp >>= 4) != 0) {
1420: int j = 0;
1421: while (exp > 1) {
1422: if ((exp & 1) != 0)
1423: dValue *= tiny10pow[j];
1424: j++;
1425: exp >>= 1;
1426: }
1427: /*
1428: * The reason for the exp > 1 condition
1429: * in the above loop was so that the last multiply
1430: * would get unrolled. We handle it here.
1431: * It could underflow.
1432: */
1433: double t = dValue * tiny10pow[j];
1434: if (t == 0.0) {
1435: /*
1436: * It did underflow.
1437: * Look more closely at the result.
1438: * If the exponent is just one too small,
1439: * then use the minimum finite as our estimate
1440: * value. Else call the result 0.0
1441: * and punt it.
1442: * ( I presume this could happen because
1443: * rounding forces the result here to be
1444: * an ULP or two less than
1445: * Double.MIN_VALUE ).
1446: */
1447: t = dValue * 2.0;
1448: t *= tiny10pow[j];
1449: if (t == 0.0) {
1450: return (isNegative) ? -0.0 : 0.0;
1451: }
1452: t = Double.MIN_VALUE;
1453: }
1454: dValue = t;
1455: }
1456: }
1457:
1458: /*
1459: * dValue is now approximately the result.
1460: * The hard part is adjusting it, by comparison
1461: * with FDBigInt arithmetic.
1462: * Formulate the EXACT big-number result as
1463: * bigD0 * 10^exp
1464: */
1465: FDBigInt bigD0 = new FDBigInt(lValue, digits, kDigits, nDigits);
1466: exp = decExponent - nDigits;
1467:
1468: correctionLoop: while (true) {
1469: /* AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES
1470: * bigIntExp and bigIntNBits
1471: */
1472: FDBigInt bigB = doubleToBigInt(dValue);
1473:
1474: /*
1475: * Scale bigD, bigB appropriately for
1476: * big-integer operations.
1477: * Naively, we multiply by powers of ten
1478: * and powers of two. What we actually do
1479: * is keep track of the powers of 5 and
1480: * powers of 2 we would use, then factor out
1481: * common divisors before doing the work.
1482: */
1483: int B2, B5; // powers of 2, 5 in bigB
1484: int D2, D5; // powers of 2, 5 in bigD
1485: int Ulp2; // powers of 2 in halfUlp.
1486: if (exp >= 0) {
1487: B2 = B5 = 0;
1488: D2 = D5 = exp;
1489: } else {
1490: B2 = B5 = -exp;
1491: D2 = D5 = 0;
1492: }
1493: if (bigIntExp >= 0) {
1494: B2 += bigIntExp;
1495: } else {
1496: D2 -= bigIntExp;
1497: }
1498: Ulp2 = B2;
1499: // shift bigB and bigD left by a number s. t.
1500: // halfUlp is still an integer.
1501: int hulpbias;
1502: if (bigIntExp + bigIntNBits <= -expBias + 1) {
1503: // This is going to be a denormalized number
1504: // (if not actually zero).
1505: // half an ULP is at 2^-(expBias+expShift+1)
1506: hulpbias = bigIntExp + expBias + expShift;
1507: } else {
1508: hulpbias = expShift + 2 - bigIntNBits;
1509: }
1510: B2 += hulpbias;
1511: D2 += hulpbias;
1512: // if there are common factors of 2, we might just as well
1513: // factor them out, as they add nothing useful.
1514: int common2 = Math.min(B2, Math.min(D2, Ulp2));
1515: B2 -= common2;
1516: D2 -= common2;
1517: Ulp2 -= common2;
1518: // do multiplications by powers of 5 and 2
1519: bigB = multPow52(bigB, B5, B2);
1520: FDBigInt bigD = multPow52(new FDBigInt(bigD0), D5, D2);
1521: //
1522: // to recap:
1523: // bigB is the scaled-big-int version of our floating-point
1524: // candidate.
1525: // bigD is the scaled-big-int version of the exact value
1526: // as we understand it.
1527: // halfUlp is 1/2 an ulp of bigB, except for special cases
1528: // of exact powers of 2
1529: //
1530: // the plan is to compare bigB with bigD, and if the difference
1531: // is less than halfUlp, then we're satisfied. Otherwise,
1532: // use the ratio of difference to halfUlp to calculate a fudge
1533: // factor to add to the floating value, then go 'round again.
1534: //
1535: FDBigInt diff;
1536: int cmpResult;
1537: boolean overvalue;
1538: if ((cmpResult = bigB.cmp(bigD)) > 0) {
1539: overvalue = true; // our candidate is too big.
1540: diff = bigB.sub(bigD);
1541: if ((bigIntNBits == 1) && (bigIntExp > -expBias)) {
1542: // candidate is a normalized exact power of 2 and
1543: // is too big. We will be subtracting.
1544: // For our purposes, ulp is the ulp of the
1545: // next smaller range.
1546: Ulp2 -= 1;
1547: if (Ulp2 < 0) {
1548: // rats. Cannot de-scale ulp this far.
1549: // must scale diff in other direction.
1550: Ulp2 = 0;
1551: diff.lshiftMe(1);
1552: }
1553: }
1554: } else if (cmpResult < 0) {
1555: overvalue = false; // our candidate is too small.
1556: diff = bigD.sub(bigB);
1557: } else {
1558: // the candidate is exactly right!
1559: // this happens with surprising frequency
1560: break correctionLoop;
1561: }
1562: FDBigInt halfUlp = constructPow52(B5, Ulp2);
1563: if ((cmpResult = diff.cmp(halfUlp)) < 0) {
1564: // difference is small.
1565: // this is close enough
1566: roundDir = overvalue ? -1 : 1;
1567: break correctionLoop;
1568: } else if (cmpResult == 0) {
1569: // difference is exactly half an ULP
1570: // round to some other value maybe, then finish
1571: dValue += 0.5 * ulp(dValue, overvalue);
1572: // should check for bigIntNBits == 1 here??
1573: roundDir = overvalue ? -1 : 1;
1574: break correctionLoop;
1575: } else {
1576: // difference is non-trivial.
1577: // could scale addend by ratio of difference to
1578: // halfUlp here, if we bothered to compute that difference.
1579: // Most of the time ( I hope ) it is about 1 anyway.
1580: dValue += ulp(dValue, overvalue);
1581: if (dValue == 0.0 || dValue == Double.POSITIVE_INFINITY)
1582: break correctionLoop; // oops. Fell off end of range.
1583: continue; // try again.
1584: }
1585:
1586: }
1587: return (isNegative) ? -dValue : dValue;
1588: }
1589:
1590: /*
1591: * Take a FloatingDecimal, which we presumably just scanned in,
1592: * and find out what its value is, as a float.
1593: * This is distinct from doubleValue() to avoid the
1594: * case of a double rounding error, wherein the conversion
1595: * to double has one rounding error, and the conversion of that double
1596: * to a float has another rounding error, IN THE WRONG DIRECTION,
1597: * ( because of the preference to a zero low-order bit ).
1598: */
1599:
1600: public float floatValue() {
1601: int kDigits = Math.min(nDigits, singleMaxDecimalDigits + 1);
1602: int iValue;
1603: float fValue;
1604:
1605: /*
1606: * convert the lead kDigits to an integer.
1607: */
1608: iValue = (int) digits[0] - (int) '0';
1609: for (int i = 1; i < kDigits; i++) {
1610: iValue = iValue * 10 + (int) digits[i] - (int) '0';
1611: }
1612: fValue = (float) iValue;
1613: int exp = decExponent - kDigits;
1614: /*
1615: * iValue now contains an integer with the value of
1616: * the first kDigits digits of the number.
1617: * fValue contains the (float) of the same.
1618: */
1619:
1620: if (nDigits <= singleMaxDecimalDigits) {
1621: /*
1622: * possibly an easy case.
1623: * We know that the digits can be represented
1624: * exactly. And if the exponent isn't too outrageous,
1625: * the whole thing can be done with one operation,
1626: * thus one rounding error.
1627: * Note that all our constructors trim all leading and
1628: * trailing zeros, so simple values (including zero)
1629: * will always end up here.
1630: */
1631: if (exp == 0 || fValue == 0.0f)
1632: return (isNegative) ? -fValue : fValue; // small floating integer
1633: else if (exp >= 0) {
1634: if (exp <= singleMaxSmallTen) {
1635: /*
1636: * Can get the answer with one operation,
1637: * thus one roundoff.
1638: */
1639: fValue *= singleSmall10pow[exp];
1640: return (isNegative) ? -fValue : fValue;
1641: }
1642: int slop = singleMaxDecimalDigits - kDigits;
1643: if (exp <= singleMaxSmallTen + slop) {
1644: /*
1645: * We can multiply dValue by 10^(slop)
1646: * and it is still "small" and exact.
1647: * Then we can multiply by 10^(exp-slop)
1648: * with one rounding.
1649: */
1650: fValue *= singleSmall10pow[slop];
1651: fValue *= singleSmall10pow[exp - slop];
1652: return (isNegative) ? -fValue : fValue;
1653: }
1654: /*
1655: * Else we have a hard case with a positive exp.
1656: */
1657: } else {
1658: if (exp >= -singleMaxSmallTen) {
1659: /*
1660: * Can get the answer in one division.
1661: */
1662: fValue /= singleSmall10pow[-exp];
1663: return (isNegative) ? -fValue : fValue;
1664: }
1665: /*
1666: * Else we have a hard case with a negative exp.
1667: */
1668: }
1669: } else if ((decExponent >= nDigits)
1670: && (nDigits + decExponent <= maxDecimalDigits)) {
1671: /*
1672: * In double-precision, this is an exact floating integer.
1673: * So we can compute to double, then shorten to float
1674: * with one round, and get the right answer.
1675: *
1676: * First, finish accumulating digits.
1677: * Then convert that integer to a double, multiply
1678: * by the appropriate power of ten, and convert to float.
1679: */
1680: long lValue = (long) iValue;
1681: for (int i = kDigits; i < nDigits; i++) {
1682: lValue = lValue * 10L
1683: + (long) ((int) digits[i] - (int) '0');
1684: }
1685: double dValue = (double) lValue;
1686: exp = decExponent - nDigits;
1687: dValue *= small10pow[exp];
1688: fValue = (float) dValue;
1689: return (isNegative) ? -fValue : fValue;
1690:
1691: }
1692: /*
1693: * Harder cases:
1694: * The sum of digits plus exponent is greater than
1695: * what we think we can do with one error.
1696: *
1697: * Start by weeding out obviously out-of-range
1698: * results, then convert to double and go to
1699: * common hard-case code.
1700: */
1701: if (decExponent > singleMaxDecimalExponent + 1) {
1702: /*
1703: * Lets face it. This is going to be
1704: * Infinity. Cut to the chase.
1705: */
1706: return (isNegative) ? Float.NEGATIVE_INFINITY
1707: : Float.POSITIVE_INFINITY;
1708: } else if (decExponent < singleMinDecimalExponent - 1) {
1709: /*
1710: * Lets face it. This is going to be
1711: * zero. Cut to the chase.
1712: */
1713: return (isNegative) ? -0.0f : 0.0f;
1714: }
1715:
1716: /*
1717: * Here, we do 'way too much work, but throwing away
1718: * our partial results, and going and doing the whole
1719: * thing as double, then throwing away half the bits that computes
1720: * when we convert back to float.
1721: *
1722: * The alternative is to reproduce the whole multiple-precision
1723: * algorithm for float precision, or to try to parameterize it
1724: * for common usage. The former will take about 400 lines of code,
1725: * and the latter I tried without success. Thus we use the following
1726: * solution here.
1727: */
1728: mustSetRoundDir = true;
1729: double dValue = doubleValue();
1730: return stickyRound(dValue);
1731: }
1732:
1733: /*
1734: * All the positive powers of 10 that can be
1735: * represented exactly in double/float.
1736: */
1737: private static final double small10pow[] = { 1.0e0, 1.0e1, 1.0e2,
1738: 1.0e3, 1.0e4, 1.0e5, 1.0e6, 1.0e7, 1.0e8, 1.0e9, 1.0e10,
1739: 1.0e11, 1.0e12, 1.0e13, 1.0e14, 1.0e15, 1.0e16, 1.0e17,
1740: 1.0e18, 1.0e19, 1.0e20, 1.0e21, 1.0e22 };
1741:
1742: private static final float singleSmall10pow[] = { 1.0e0f, 1.0e1f,
1743: 1.0e2f, 1.0e3f, 1.0e4f, 1.0e5f, 1.0e6f, 1.0e7f, 1.0e8f,
1744: 1.0e9f, 1.0e10f };
1745:
1746: private static final double big10pow[] = { 1e16, 1e32, 1e64, 1e128,
1747: 1e256 };
1748: private static final double tiny10pow[] = { 1e-16, 1e-32, 1e-64,
1749: 1e-128, 1e-256 };
1750:
1751: private static final int maxSmallTen = small10pow.length - 1;
1752: private static final int singleMaxSmallTen = singleSmall10pow.length - 1;
1753:
1754: private static final int small5pow[] = { 1, 5, 5 * 5, 5 * 5 * 5,
1755: 5 * 5 * 5 * 5, 5 * 5 * 5 * 5 * 5, 5 * 5 * 5 * 5 * 5 * 5,
1756: 5 * 5 * 5 * 5 * 5 * 5 * 5, 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1757: 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1758: 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1759: 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1760: 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1761: 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 };
1762:
1763: private static final long long5pow[] = {
1764: 1L,
1765: 5L,
1766: 5L * 5,
1767: 5L * 5 * 5,
1768: 5L * 5 * 5 * 5,
1769: 5L * 5 * 5 * 5 * 5,
1770: 5L * 5 * 5 * 5 * 5 * 5,
1771: 5L * 5 * 5 * 5 * 5 * 5 * 5,
1772: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1773: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1774: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1775: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1776: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1777: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1778: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1779: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1780: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5
1781: * 5,
1782: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5
1783: * 5 * 5,
1784: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5
1785: * 5 * 5 * 5,
1786: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5
1787: * 5 * 5 * 5 * 5,
1788: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5
1789: * 5 * 5 * 5 * 5 * 5,
1790: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5
1791: * 5 * 5 * 5 * 5 * 5 * 5,
1792: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5
1793: * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1794: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5
1795: * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1796: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5
1797: * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1798: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5
1799: * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
1800: 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5
1801: * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5, };
1802:
1803: // approximately ceil( log2( long5pow[i] ) )
1804: private static final int n5bits[] = { 0, 3, 5, 7, 10, 12, 14, 17,
1805: 19, 21, 24, 26, 28, 31, 33, 35, 38, 40, 42, 45, 47, 49, 52,
1806: 54, 56, 59, 61, };
1807:
1808: private static final char infinity[] = { 'I', 'n', 'f', 'i', 'n',
1809: 'i', 't', 'y' };
1810: private static final char notANumber[] = { 'N', 'a', 'N' };
1811: private static final char zero[] = { '0', '0', '0', '0', '0', '0',
1812: '0', '0' };
1813:
1814: }
1815:
1816: /*
1817: * A really, really simple bigint package
1818: * tailored to the needs of floating base conversion.
1819: */
1820: class FDBigInt {
1821: int nWords; // number of words used
1822: int data[]; // value: data[0] is least significant
1823:
1824: /* public */FDBigInt(int v) {
1825: nWords = 1;
1826: data = new int[1];
1827: data[0] = v;
1828: }
1829:
1830: /* public */FDBigInt(long v) {
1831: data = new int[2];
1832: data[0] = (int) v;
1833: data[1] = (int) (v >>> 32);
1834: nWords = (data[1] == 0) ? 1 : 2;
1835: }
1836:
1837: /* public */FDBigInt(FDBigInt other) {
1838: data = new int[nWords = other.nWords];
1839: JVM.unchecked_int_arraycopy(other.data, 0, data, 0, nWords);
1840: }
1841:
1842: private FDBigInt(int[] d, int n) {
1843: data = d;
1844: nWords = n;
1845: }
1846:
1847: /* public */FDBigInt(long seed, char digit[], int nd0, int nd) {
1848: int n = (nd + 8) / 9; // estimate size needed.
1849: if (n < 2)
1850: n = 2;
1851: data = new int[n]; // allocate enough space
1852: data[0] = (int) seed; // starting value
1853: data[1] = (int) (seed >>> 32);
1854: nWords = (data[1] == 0) ? 1 : 2;
1855: int i = nd0;
1856: int limit = nd - 5; // slurp digits 5 at a time.
1857: int v;
1858: while (i < limit) {
1859: int ilim = i + 5;
1860: v = (int) digit[i++] - (int) '0';
1861: while (i < ilim) {
1862: v = 10 * v + (int) digit[i++] - (int) '0';
1863: }
1864: multaddMe(100000, v); // ... where 100000 is 10^5.
1865: }
1866: int factor = 1;
1867: v = 0;
1868: while (i < nd) {
1869: v = 10 * v + (int) digit[i++] - (int) '0';
1870: factor *= 10;
1871: }
1872: if (factor != 1) {
1873: multaddMe(factor, v);
1874: }
1875: }
1876:
1877: /*
1878: * Left shift by c bits.
1879: * Shifts this in place.
1880: */
1881: public void lshiftMe(int c) throws IllegalArgumentException {
1882: if (c <= 0) {
1883: if (c == 0)
1884: return; // silly.
1885: else {
1886: throw new IllegalArgumentException(
1887: /* #ifdef VERBOSE_EXCEPTIONS */
1888: /// skipped "negative shift count"
1889: /* #endif */
1890: );
1891: }
1892: }
1893: int wordcount = c >> 5;
1894: int bitcount = c & 0x1f;
1895: int anticount = 32 - bitcount;
1896: int t[] = data;
1897: int s[] = data;
1898: if (nWords + wordcount + 1 > t.length) {
1899: // reallocate.
1900: t = new int[nWords + wordcount + 1];
1901: }
1902: int target = nWords + wordcount;
1903: int src = nWords - 1;
1904: if (bitcount == 0) {
1905: // special solution, since an anticount of 32 won't go!
1906: JVM.unchecked_int_arraycopy(s, 0, t, wordcount, nWords);
1907: target = wordcount - 1;
1908: } else {
1909: t[target--] = s[src] >>> anticount;
1910: while (src >= 1) {
1911: t[target--] = (s[src] << bitcount)
1912: | (s[--src] >>> anticount);
1913: }
1914: t[target--] = s[src] << bitcount;
1915: }
1916: while (target >= 0) {
1917: t[target--] = 0;
1918: }
1919: data = t;
1920: nWords += wordcount + 1;
1921: // may have constructed high-order word of 0.
1922: // if so, trim it
1923: while (nWords > 1 && data[nWords - 1] == 0)
1924: nWords--;
1925: }
1926:
1927: /*
1928: * normalize this number by shifting until
1929: * the MSB of the number is at 0x08000000.
1930: * This is in preparation for quoRemIteration, below.
1931: * The idea is that, to make division easier, we want the
1932: * divisor to be "normalized" -- usually this means shifting
1933: * the MSB into the high words sign bit. But because we know that
1934: * the quotient will be 0 < q < 10, we would like to arrange that
1935: * the dividend not span up into another word of precision.
1936: * (This needs to be explained more clearly!)
1937: */
1938: public int normalizeMe() throws IllegalArgumentException {
1939: int src;
1940: int wordcount = 0;
1941: int bitcount = 0;
1942: int v = 0;
1943: for (src = nWords - 1; src >= 0 && (v = data[src]) == 0; src--) {
1944: wordcount += 1;
1945: }
1946: if (src < 0) {
1947: // oops. Value is zero. Cannot normalize it!
1948: throw new IllegalArgumentException(
1949: /* #ifdef VERBOSE_EXCEPTIONS */
1950: /// skipped "zero value"
1951: /* #endif */
1952: );
1953: }
1954: /*
1955: * In most cases, we assume that wordcount is zero. This only
1956: * makes sense, as we try not to maintain any high-order
1957: * words full of zeros. In fact, if there are zeros, we will
1958: * simply SHORTEN our number at this point. Watch closely...
1959: */
1960: nWords -= wordcount;
1961: /*
1962: * Compute how far left we have to shift v s.t. its highest-
1963: * order bit is in the right place. Then call lshiftMe to
1964: * do the work.
1965: */
1966: if ((v & 0xf0000000) != 0) {
1967: // will have to shift up into the next word.
1968: // too bad.
1969: bitcount = 32;
1970: while ((v & 0xf0000000) != 0) {
1971: v >>>= 1;
1972: bitcount--;
1973: }
1974: } else {
1975: while (v <= 0x000fffff) {
1976: // solution: byte-at-a-time shifting
1977: v <<= 8;
1978: bitcount += 8;
1979: }
1980: while (v <= 0x07ffffff) {
1981: v <<= 1;
1982: bitcount += 1;
1983: }
1984: }
1985: if (bitcount != 0)
1986: lshiftMe(bitcount);
1987: return bitcount;
1988: }
1989:
1990: /*
1991: * Multiply a FDBigInt by an int.
1992: * Result is a new FDBigInt.
1993: */
1994: public FDBigInt mult(int iv) {
1995: long v = iv;
1996: int r[];
1997: long p;
1998:
1999: // guess adequate size of r.
2000: r = new int[(v * ((long) data[nWords - 1] & 0xffffffffL) > 0xfffffffL) ? nWords + 1
2001: : nWords];
2002: p = 0L;
2003: for (int i = 0; i < nWords; i++) {
2004: p += v * ((long) data[i] & 0xffffffffL);
2005: r[i] = (int) p;
2006: p >>>= 32;
2007: }
2008: if (p == 0L) {
2009: return new FDBigInt(r, nWords);
2010: } else {
2011: r[nWords] = (int) p;
2012: return new FDBigInt(r, nWords + 1);
2013: }
2014: }
2015:
2016: /*
2017: * Multiply a FDBigInt by an int and add another int.
2018: * Result is computed in place.
2019: * Hope it fits!
2020: */
2021: public void multaddMe(int iv, int addend) {
2022: long v = iv;
2023: long p;
2024:
2025: // unroll 0th iteration, doing addition.
2026: p = v * ((long) data[0] & 0xffffffffL)
2027: + ((long) addend & 0xffffffffL);
2028: data[0] = (int) p;
2029: p >>>= 32;
2030: for (int i = 1; i < nWords; i++) {
2031: p += v * ((long) data[i] & 0xffffffffL);
2032: data[i] = (int) p;
2033: p >>>= 32;
2034: }
2035: if (p != 0L) {
2036: data[nWords] = (int) p; // will fail noisily if illegal!
2037: nWords++;
2038: }
2039: }
2040:
2041: /*
2042: * Multiply a FDBigInt by another FDBigInt.
2043: * Result is a new FDBigInt.
2044: */
2045: public FDBigInt mult(FDBigInt other) {
2046: // crudely guess adequate size for r
2047: int r[] = new int[nWords + other.nWords];
2048: int i;
2049: // I think I am promised zeros...
2050:
2051: for (i = 0; i < this .nWords; i++) {
2052: long v = (long) this .data[i] & 0xffffffffL; // UNSIGNED CONVERSION
2053: long p = 0L;
2054: int j;
2055: for (j = 0; j < other.nWords; j++) {
2056: p += ((long) r[i + j] & 0xffffffffL) + v
2057: * ((long) other.data[j] & 0xffffffffL); // UNSIGNED CONVERSIONS ALL 'ROUND.
2058: r[i + j] = (int) p;
2059: p >>>= 32;
2060: }
2061: r[i + j] = (int) p;
2062: }
2063: // compute how much of r we actually needed for all that.
2064: for (i = r.length - 1; i > 0; i--)
2065: if (r[i] != 0)
2066: break;
2067: return new FDBigInt(r, i + 1);
2068: }
2069:
2070: /*
2071: * Add one FDBigInt to another. Return a FDBigInt
2072: */
2073: public FDBigInt add(FDBigInt other) {
2074: int i;
2075: int a[], b[];
2076: int n, m;
2077: long c = 0L;
2078: // arrange such that a.nWords >= b.nWords;
2079: // n = a.nWords, m = b.nWords
2080: if (this .nWords >= other.nWords) {
2081: a = this .data;
2082: n = this .nWords;
2083: b = other.data;
2084: m = other.nWords;
2085: } else {
2086: a = other.data;
2087: n = other.nWords;
2088: b = this .data;
2089: m = this .nWords;
2090: }
2091: int r[] = new int[n];
2092: for (i = 0; i < n; i++) {
2093: c += (long) a[i] & 0xffffffffL;
2094: if (i < m) {
2095: c += (long) b[i] & 0xffffffffL;
2096: }
2097: r[i] = (int) c;
2098: c >>= 32; // signed shift.
2099: }
2100: if (c != 0L) {
2101: // oops -- carry out -- need longer result.
2102: int s[] = new int[r.length + 1];
2103: JVM.unchecked_int_arraycopy(r, 0, s, 0, r.length);
2104: s[i++] = (int) c;
2105: return new FDBigInt(s, i);
2106: }
2107: return new FDBigInt(r, i);
2108: }
2109:
2110: /*
2111: * Subtract one FDBigInt from another. Return a FDBigInt
2112: * Assert that the result is positive.
2113: */
2114: public FDBigInt sub(FDBigInt other) {
2115: int r[] = new int[this .nWords];
2116: int i;
2117: int n = this .nWords;
2118: int m = other.nWords;
2119: int nzeros = 0;
2120: long c = 0L;
2121: for (i = 0; i < n; i++) {
2122: c += (long) this .data[i] & 0xffffffffL;
2123: if (i < m) {
2124: c -= (long) other.data[i] & 0xffffffffL;
2125: }
2126: if ((r[i] = (int) c) == 0)
2127: nzeros++;
2128: else
2129: nzeros = 0;
2130: c >>= 32; // signed shift.
2131: }
2132: if (c != 0L) {
2133: throw new RuntimeException(
2134: /* #ifdef VERBOSE_EXCEPTIONS */
2135: /// skipped "Assertion botch: borrow out of subtract"
2136: /* #endif */
2137: );
2138: }
2139: while (i < m)
2140: if (other.data[i++] != 0) {
2141: throw new RuntimeException(
2142: /* #ifdef VERBOSE_EXCEPTIONS */
2143: /// skipped "Assertion botch: negative result of subtract"
2144: /* #endif */
2145: );
2146: }
2147: return new FDBigInt(r, n - nzeros);
2148: }
2149:
2150: /*
2151: * Compare FDBigInt with another FDBigInt. Return an integer
2152: * >0: this > other
2153: * 0: this == other
2154: * <0: this < other
2155: */
2156: public int cmp(FDBigInt other) {
2157: int i;
2158: if (this .nWords > other.nWords) {
2159: // if any of my high-order words is non-zero,
2160: // then the answer is evident
2161: int j = other.nWords - 1;
2162: for (i = this .nWords - 1; i > j; i--)
2163: if (this .data[i] != 0)
2164: return 1;
2165: } else if (this .nWords < other.nWords) {
2166: // if any of other's high-order words is non-zero,
2167: // then the answer is evident
2168: int j = this .nWords - 1;
2169: for (i = other.nWords - 1; i > j; i--)
2170: if (other.data[i] != 0)
2171: return -1;
2172: } else {
2173: i = this .nWords - 1;
2174: }
2175: for (; i > 0; i--)
2176: if (this .data[i] != other.data[i])
2177: break;
2178: // careful! want unsigned compare!
2179: // use brute force here.
2180: int a = this .data[i];
2181: int b = other.data[i];
2182: if (a < 0) {
2183: // a is really big, unsigned
2184: if (b < 0) {
2185: return a - b; // both big, negative
2186: } else {
2187: return 1; // b not big, answer is obvious;
2188: }
2189: } else {
2190: // a is not really big
2191: if (b < 0) {
2192: // but b is really big
2193: return -1;
2194: } else {
2195: return a - b;
2196: }
2197: }
2198: }
2199:
2200: /*
2201: * Compute
2202: * q = (int)( this / S )
2203: * this = 10 * ( this mod S )
2204: * Return q.
2205: * This is the iteration step of digit development for output.
2206: * We assume that S has been normalized, as above, and that
2207: * "this" has been lshift'ed accordingly.
2208: * Also assume, of course, that the result, q, can be expressed
2209: * as an integer, 0 <= q < 10.
2210: */
2211: public int quoRemIteration(FDBigInt S)
2212: throws IllegalArgumentException {
2213: // ensure that this and S have the same number of
2214: // digits. If S is properly normalized and q < 10 then
2215: // this must be so.
2216: if (nWords != S.nWords) {
2217: throw new IllegalArgumentException(
2218: /* #ifdef VERBOSE_EXCEPTIONS */
2219: /// skipped "disparate values"
2220: /* #endif */
2221: );
2222: }
2223: // estimate q the obvious way. We will usually be
2224: // right. If not, then we're only off by a little and
2225: // will re-add.
2226: int n = nWords - 1;
2227: long q = ((long) data[n] & 0xffffffffL) / (long) S.data[n];
2228: long diff = 0L;
2229: for (int i = 0; i <= n; i++) {
2230: diff += ((long) data[i] & 0xffffffffL) - q
2231: * ((long) S.data[i] & 0xffffffffL);
2232: data[i] = (int) diff;
2233: diff >>= 32; // N.B. SIGNED shift.
2234: }
2235: if (diff != 0L) {
2236: // q is too big.
2237: // add S back in until this turns +. This should
2238: // not be very many times!
2239: long sum = 0L;
2240: while (sum == 0L) {
2241: sum = 0L;
2242: for (int i = 0; i <= n; i++) {
2243: sum += ((long) data[i] & 0xffffffffL)
2244: + ((long) S.data[i] & 0xffffffffL);
2245: data[i] = (int) sum;
2246: sum >>= 32; // Signed or unsigned, answer is 0 or 1
2247: }
2248: /*
2249: * Originally the following line read
2250: * "if ( sum !=0 && sum != -1 )"
2251: * but that would be wrong, because of the
2252: * treatment of the two values as entirely unsigned,
2253: * it would be impossible for a carry-out to be interpreted
2254: * as -1 -- it would have to be a single-bit carry-out, or
2255: * +1.
2256: */
2257: if (sum != 0 && sum != 1) {
2258: throw new RuntimeException(
2259: /* #ifdef VERBOSE_EXCEPTIONS */
2260: /// skipped "Assertion botch: "+sum+" carry out of division correction"
2261: /* #endif */
2262: );
2263: }
2264: q -= 1;
2265: }
2266: }
2267: // finally, we can multiply this by 10.
2268: // it cannot overflow, right, as the high-order word has
2269: // at least 4 high-order zeros!
2270: long p = 0L;
2271: for (int i = 0; i <= n; i++) {
2272: p += 10 * ((long) data[i] & 0xffffffffL);
2273: data[i] = (int) p;
2274: p >>= 32; // SIGNED shift.
2275: }
2276: if (p != 0L) {
2277: throw new RuntimeException(
2278: /* #ifdef VERBOSE_EXCEPTIONS */
2279: /// skipped "Assertion botch: carry out of *10"
2280: /* #endif */
2281: );
2282: }
2283:
2284: return (int) q;
2285: }
2286:
2287: public long longValue() {
2288: // if this can be represented as a long,
2289: // return the value
2290: int i;
2291: for (i = this .nWords - 1; i > 1; i--) {
2292: if (data[i] != 0) {
2293: throw new RuntimeException(
2294: /* #ifdef VERBOSE_EXCEPTIONS */
2295: /// skipped "Assertion botch: value too big"
2296: /* #endif */
2297: );
2298: }
2299: }
2300: switch (i) {
2301: case 1:
2302: if (data[1] < 0) {
2303: throw new RuntimeException(
2304: /* #ifdef VERBOSE_EXCEPTIONS */
2305: /// skipped "Assertion botch: value too big"
2306: /* #endif */
2307: );
2308: }
2309: return ((long) (data[1]) << 32)
2310: | ((long) data[0] & 0xffffffffL);
2311: case 0:
2312: return ((long) data[0] & 0xffffffffL);
2313: default:
2314: throw new RuntimeException(
2315: /* #ifdef VERBOSE_EXCEPTIONS */
2316: /// skipped "Assertion botch: longValue confused"
2317: /* #endif */
2318: );
2319: }
2320: }
2321:
2322: public String toString() {
2323: StringBuffer r = new StringBuffer(30);
2324: r.append('[');
2325: int i = Math.min(nWords - 1, data.length - 1);
2326: if (nWords > data.length) {
2327: r.append("(" + data.length + "<" + nWords + "!)");
2328: }
2329: for (; i > 0; i--) {
2330: r.append(Integer.toHexString(data[i]));
2331: r.append((char) ' ');
2332: }
2333: r.append(Integer.toHexString(data[0]));
2334: r.append((char) ']');
2335: return new String(r);
2336: }
2337: }
|