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