1   /*
2 * Copyright (c) 1996, 2011, Oracle and/or its affiliates. All rights reserved.
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * This code is free software; you can redistribute it and/or modify it
6 * under the terms of the GNU General Public License version 2 only, as
7 * published by the Free Software Foundation.  Oracle designates this
8 * particular file as subject to the "Classpath" exception as provided
9 * by Oracle in the LICENSE file that accompanied this code.
10 *
11 * This code is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
14 * version 2 for more details (a copy is included in the LICENSE file that
15 * accompanied this code).
16 *
17 * You should have received a copy of the GNU General Public License version
18 * 2 along with this work; if not, write to the Free Software Foundation,
19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20 *
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22 * or visit www.oracle.com if you need additional information or have any
23 * questions.
24 */
25
26package java.lang;
27
28import sun.misc.FpUtils;
29import sun.misc.FDBigInt;
30import sun.misc.DoubleConsts;
31import sun.misc.FloatConsts;
32import java.util.regex.*;
33
34public class FloatingDecimal {
35    boolean     isExceptional;
36    boolean     isNegative;
37    int         decExponent;
38    char        digits[];
39    int         nDigits;
40    int         bigIntExp;
41    int         bigIntNBits;
42    boolean     mustSetRoundDir = false;
43    boolean     fromHex = false;
44    int         roundDir = 0; // set by doubleValue
45
46    /*
47     * Constants of the implementation
48     * Most are IEEE-754 related.
49     * (There are more really boring constants at the end.)
50     */
51    static final long   signMask = 0x8000000000000000L;
52    static final long   expMask  = 0x7ff0000000000000L;
53    static final long   fractMask= ~(signMask|expMask);
54    static final int    expShift = 52;
55    static final int    expBias  = 1023;
56    static final long   fractHOB = ( 1L<<expShift ); // assumed High-Order bit
57    static final long   expOne   = ((long)expBias)<<expShift; // exponent of 1.0
58    static final int    maxSmallBinExp = 62;
59    static final int    minSmallBinExp = -( 63 / 3 );
60    static final int    maxDecimalDigits = 15;
61    static final int    maxDecimalExponent = 308;
62    static final int    minDecimalExponent = -324;
63    static final int    bigDecimalExponent = 324; // i.e. abs(minDecimalExponent)
64
65    static final long   highbyte = 0xff00000000000000L;
66    static final long   highbit  = 0x8000000000000000L;
67    static final long   lowbytes = ~highbyte;
68
69    static final int    singleSignMask =    0x80000000;
70    static final int    singleExpMask  =    0x7f800000;
71    static final int    singleFractMask =   ~(singleSignMask|singleExpMask);
72    static final int    singleExpShift  =   23;
73    static final int    singleFractHOB  =   1<<singleExpShift;
74    static final int    singleExpBias   =   127;
75    static final int    singleMaxDecimalDigits = 7;
76    static final int    singleMaxDecimalExponent = 38;
77    static final int    singleMinDecimalExponent = -45;
78
79    static final int    intDecimalDigits = 9;
80
81    /*
82     * count number of bits from high-order 1 bit to low-order 1 bit,
83     * inclusive.
84     */
85    private static int
86    countBits( long v ){
87        //
88        // the strategy is to shift until we get a non-zero sign bit
89        // then shift until we have no bits left, counting the difference.
90        // we do byte shifting as a hack. Hope it helps.
91        //
92        if ( v == 0L ) return 0;
93
94        while ( ( v & highbyte ) == 0L ){
95            v <<= 8;
96        }
97        while ( v > 0L ) { // i.e. while ((v&highbit) == 0L )
98            v <<= 1;
99        }
100
101        int n = 0;
102        while (( v & lowbytes ) != 0L ){
103            v <<= 8;
104            n += 8;
105        }
106        while ( v != 0L ){
107            v <<= 1;
108            n += 1;
109        }
110        return n;
111    }
112
113    /*
114     * Keep big powers of 5 handy for future reference.
115     */
116    private static FDBigInt b5p[];
117
118    private static synchronized FDBigInt
119    big5pow( int p ){
120        assert p >= 0 : p; // negative power of 5
121        if ( b5p == null ){
122            b5p = new FDBigInt[ p+1 ];
123        }else if (b5p.length <= p ){
124            FDBigInt t[] = new FDBigInt[ p+1 ];
125            System.arraycopy( b5p, 0, t, 0, b5p.length );
126            b5p = t;
127        }
128        if ( b5p[p] != null )
129            return b5p[p];
130        else if ( p < small5pow.length )
131            return b5p[p] = new FDBigInt( small5pow[p] );
132        else if ( p < long5pow.length )
133            return b5p[p] = new FDBigInt( long5pow[p] );
134        else {
135            // construct the value.
136            // recursively.
137            int q, r;
138            // in order to compute 5^p,
139            // compute its square root, 5^(p/2) and square.
140            // or, let q = p / 2, r = p -q, then
141            // 5^p = 5^(q+r) = 5^q * 5^r
142            q = p >> 1;
143            r = p - q;
144            FDBigInt bigq =  b5p[q];
145            if ( bigq == null )
146                bigq = big5pow ( q );
147            if ( r < small5pow.length ){
148                return (b5p[p] = bigq.mult( small5pow[r] ) );
149            }else{
150                FDBigInt bigr = b5p[ r ];
151                if ( bigr == null )
152                    bigr = big5pow( r );
153                return (b5p[p] = bigq.mult( bigr ) );
154            }
155        }
156    }
157
158    //
159    // a common operation
160    //
161    private static FDBigInt
162    multPow52( FDBigInt v, int p5, int p2 ){
163        if ( p5 != 0 ){
164            if ( p5 < small5pow.length ){
165                v = v.mult( small5pow[p5] );
166            } else {
167                v = v.mult( big5pow( p5 ) );
168            }
169        }
170        if ( p2 != 0 ){
171            v.lshiftMe( p2 );
172        }
173        return v;
174    }
175
176    //
177    // another common operation
178    //
179    private static FDBigInt
180    constructPow52( int p5, int p2 ){
181        FDBigInt v = new FDBigInt( big5pow( p5 ) );
182        if ( p2 != 0 ){
183            v.lshiftMe( p2 );
184        }
185        return v;
186    }
187
188    /*
189     * Make a floating double into a FDBigInt.
190     * This could also be structured as a FDBigInt
191     * constructor, but we'd have to build a lot of knowledge
192     * about floating-point representation into it, and we don't want to.
193     *
194     * AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES
195     * bigIntExp and bigIntNBits
196     *
197     */
198    private FDBigInt
199    doubleToBigInt( double dval ){
200        long lbits = Double.doubleToLongBits( dval ) & ~signMask;
201        int binexp = (int)(lbits >>> expShift);
202        lbits &= fractMask;
203        if ( binexp > 0 ){
204            lbits |= fractHOB;
205        } else {
206            assert lbits != 0L : lbits; // doubleToBigInt(0.0)
207            binexp +=1;
208            while ( (lbits & fractHOB ) == 0L){
209                lbits <<= 1;
210                binexp -= 1;
211            }
212        }
213        binexp -= expBias;
214        int nbits = countBits( lbits );
215        /*
216         * We now know where the high-order 1 bit is,
217         * and we know how many there are.
218         */
219        int lowOrderZeros = expShift+1-nbits;
220        lbits >>>= lowOrderZeros;
221
222        bigIntExp = binexp+1-nbits;
223        bigIntNBits = nbits;
224        return new FDBigInt( lbits );
225    }
226
227    /*
228     * Compute a number that is the ULP of the given value,
229     * for purposes of addition/subtraction. Generally easy.
230     * More difficult if subtracting and the argument
231     * is a normalized a power of 2, as the ULP changes at these points.
232     */
233    private static double ulp( double dval, boolean subtracting ){
234        long lbits = Double.doubleToLongBits( dval ) & ~signMask;
235        int binexp = (int)(lbits >>> expShift);
236        double ulpval;
237        if ( subtracting && ( binexp >= expShift ) && ((lbits&fractMask) == 0L) ){
238            // for subtraction from normalized, powers of 2,
239            // use next-smaller exponent
240            binexp -= 1;
241        }
242        if ( binexp > expShift ){
243            ulpval = Double.longBitsToDouble( ((long)(binexp-expShift))<<expShift );
244        } else if ( binexp == 0 ){
245            ulpval = Double.MIN_VALUE;
246        } else {
247            ulpval = Double.longBitsToDouble( 1L<<(binexp-1) );
248        }
249        if ( subtracting ) ulpval = - ulpval;
250
251        return ulpval;
252    }
253
254    /*
255     * Round a double to a float.
256     * In addition to the fraction bits of the double,
257     * look at the class instance variable roundDir,
258     * which should help us avoid double-rounding error.
259     * roundDir was set in hardValueOf if the estimate was
260     * close enough, but not exact. It tells us which direction
261     * of rounding is preferred.
262     */
263    float
264    stickyRound( double dval ){
265        long lbits = Double.doubleToLongBits( dval );
266        long binexp = lbits & expMask;
267        if ( binexp == 0L || binexp == expMask ){
268            // what we have here is special.
269            // don't worry, the right thing will happen.
270            return (float) dval;
271        }
272        lbits += (long)roundDir; // hack-o-matic.
273        return (float)Double.longBitsToDouble( lbits );
274    }
275
276
277    /*
278     * This is the easy subcase --
279     * all the significant bits, after scaling, are held in lvalue.
280     * negSign and decExponent tell us what processing and scaling
281     * has already been done. Exceptional cases have already been
282     * stripped out.
283     * In particular:
284     * lvalue is a finite number (not Inf, nor NaN)
285     * lvalue > 0L (not zero, nor negative).
286     *
287     * The only reason that we develop the digits here, rather than
288     * calling on Long.toString() is that we can do it a little faster,
289     * and besides want to treat trailing 0s specially. If Long.toString
290     * changes, we should re-evaluate this strategy!
291     */
292    private void
293    developLongDigits( int decExponent, long lvalue, long insignificant ){
294        char digits[];
295        int  ndigits;
296        int  digitno;
297        int  c;
298        //
299        // Discard non-significant low-order bits, while rounding,
300        // up to insignificant value.
301        int i;
302        for ( i = 0; insignificant >= 10L; i++ )
303            insignificant /= 10L;
304        if ( i != 0 ){
305            long pow10 = long5pow[i] << i; // 10^i == 5^i * 2^i;
306            long residue = lvalue % pow10;
307            lvalue /= pow10;
308            decExponent += i;
309            if ( residue >= (pow10>>1) ){
310                // round up based on the low-order bits we're discarding
311                lvalue++;
312            }
313        }
314        if ( lvalue <= Integer.MAX_VALUE ){
315            assert lvalue > 0L : lvalue; // lvalue <= 0
316            // even easier subcase!
317            // can do int arithmetic rather than long!
318            int  ivalue = (int)lvalue;
319            ndigits = 10;
320            digits = (char[])(perThreadBuffer.get());
321            digitno = ndigits-1;
322            c = ivalue%10;
323            ivalue /= 10;
324            while ( c == 0 ){
325                decExponent++;
326                c = ivalue%10;
327                ivalue /= 10;
328            }
329            while ( ivalue != 0){
330                digits[digitno--] = (char)(c+'0');
331                decExponent++;
332                c = ivalue%10;
333                ivalue /= 10;
334            }
335            digits[digitno] = (char)(c+'0');
336        } else {
337            // same algorithm as above (same bugs, too )
338            // but using long arithmetic.
339            ndigits = 20;
340            digits = (char[])(perThreadBuffer.get());
341            digitno = ndigits-1;
342            c = (int)(lvalue%10L);
343            lvalue /= 10L;
344            while ( c == 0 ){
345                decExponent++;
346                c = (int)(lvalue%10L);
347                lvalue /= 10L;
348            }
349            while ( lvalue != 0L ){
350                digits[digitno--] = (char)(c+'0');
351                decExponent++;
352                c = (int)(lvalue%10L);
353                lvalue /= 10;
354            }
355            digits[digitno] = (char)(c+'0');
356        }
357        char result [];
358        ndigits -= digitno;
359        result = new char[ ndigits ];
360        System.arraycopy( digits, digitno, result, 0, ndigits );
361        this.digits = result;
362        this.decExponent = decExponent+1;
363        this.nDigits = ndigits;
364    }
365
366    //
367    // add one to the least significant digit.
368    // in the unlikely event there is a carry out,
369    // deal with it.
370    // assert that this will only happen where there
371    // is only one digit, e.g. (float)1e-44 seems to do it.
372    //
373    private void
374    roundup(){
375        int i;
376        int q = digits[ i = (nDigits-1)];
377        if ( q == '9' ){
378            while ( q == '9' && i > 0 ){
379                digits[i] = '0';
380                q = digits[--i];
381            }
382            if ( q == '9' ){
383                // carryout! High-order 1, rest 0s, larger exp.
384                decExponent += 1;
385                digits[0] = '1';
386                return;
387            }
388            // else fall through.
389        }
390        digits[i] = (char)(q+1);
391    }
392
393    private FloatingDecimal() {}
394
395    private static final ThreadLocal<FloatingDecimal> TL_INSTANCE = new ThreadLocal<FloatingDecimal>() {
396        @Override protected FloatingDecimal initialValue() {
397            return new FloatingDecimal();
398        }
399    };
400    public static FloatingDecimal getThreadLocalInstance() {
401        return TL_INSTANCE.get();
402    }
403
404    /*
405     * FIRST IMPORTANT LOAD: DOUBLE
406     */
407    public FloatingDecimal loadDouble(double d) {
408        long    dBits = Double.doubleToLongBits( d );
409        long    fractBits;
410        int     binExp;
411        int     nSignificantBits;
412
413        mustSetRoundDir = false;
414        fromHex = false;
415        roundDir = 0;
416
417        // discover and delete sign
418        if ( (dBits&signMask) != 0 ){
419            isNegative = true;
420            dBits ^= signMask;
421        } else {
422            isNegative = false;
423        }
424        // Begin to unpack
425        // Discover obvious special cases of NaN and Infinity.
426        binExp = (int)( (dBits&expMask) >> expShift );
427        fractBits = dBits&fractMask;
428        if ( binExp == (int)(expMask>>expShift) ) {
429            isExceptional = true;
430            if ( fractBits == 0L ){
431                digits =  infinity;
432            } else {
433                digits = notANumber;
434                isNegative = false; // NaN has no sign!
435            }
436            nDigits = digits.length;
437            return this;
438        }
439        isExceptional = false;
440        // Finish unpacking
441        // Normalize denormalized numbers.
442        // Insert assumed high-order bit for normalized numbers.
443        // Subtract exponent bias.
444        if ( binExp == 0 ){
445            if ( fractBits == 0L ){
446                // not a denorm, just a 0!
447                decExponent = 0;
448                digits = zero;
449                nDigits = 1;
450                return this;
451            }
452            while ( (fractBits&fractHOB) == 0L ){
453                fractBits <<= 1;
454                binExp -= 1;
455            }
456            nSignificantBits = expShift + binExp +1; // recall binExp is  - shift count.
457            binExp += 1;
458        } else {
459            fractBits |= fractHOB;
460            nSignificantBits = expShift+1;
461        }
462        binExp -= expBias;
463        // call the routine that actually does all the hard work.
464        dtoa( binExp, fractBits, nSignificantBits );
465        return this;
466    }
467
468    /*
469     * SECOND IMPORTANT LOAD: SINGLE
470     */
471    public FloatingDecimal loadFloat(float f) {
472        int     fBits = Float.floatToIntBits( f );
473        int     fractBits;
474        int     binExp;
475        int     nSignificantBits;
476
477        mustSetRoundDir = false;
478        fromHex = false;
479        roundDir = 0;
480
481        // discover and delete sign
482        if ( (fBits&singleSignMask) != 0 ){
483            isNegative = true;
484            fBits ^= singleSignMask;
485        } else {
486            isNegative = false;
487        }
488        // Begin to unpack
489        // Discover obvious special cases of NaN and Infinity.
490        binExp = (int)( (fBits&singleExpMask) >> singleExpShift );
491        fractBits = fBits&singleFractMask;
492        if ( binExp == (int)(singleExpMask>>singleExpShift) ) {
493            isExceptional = true;
494            if ( fractBits == 0L ){
495                digits =  infinity;
496            } else {
497                digits = notANumber;
498                isNegative = false; // NaN has no sign!
499            }
500            nDigits = digits.length;
501            return this;
502        }
503        isExceptional = false;
504        // Finish unpacking
505        // Normalize denormalized numbers.
506        // Insert assumed high-order bit for normalized numbers.
507        // Subtract exponent bias.
508        if ( binExp == 0 ){
509            if ( fractBits == 0 ){
510                // not a denorm, just a 0!
511                decExponent = 0;
512                digits = zero;
513                nDigits = 1;
514                return this;
515            }
516            while ( (fractBits&singleFractHOB) == 0 ){
517                fractBits <<= 1;
518                binExp -= 1;
519            }
520            nSignificantBits = singleExpShift + binExp +1; // recall binExp is  - shift count.
521            binExp += 1;
522        } else {
523            fractBits |= singleFractHOB;
524            nSignificantBits = singleExpShift+1;
525        }
526        binExp -= singleExpBias;
527        // call the routine that actually does all the hard work.
528        dtoa( binExp, ((long)fractBits)<<(expShift-singleExpShift), nSignificantBits );
529        return this;
530    }
531
532    private void
533    dtoa( int binExp, long fractBits, int nSignificantBits )
534    {
535        int     nFractBits; // number of significant bits of fractBits;
536        int     nTinyBits;  // number of these to the right of the point.
537        int     decExp;
538
539        // Examine number. Determine if it is an easy case,
540        // which we can do pretty trivially using float/long conversion,
541        // or whether we must do real work.
542        nFractBits = countBits( fractBits );
543        nTinyBits = Math.max( 0, nFractBits - binExp - 1 );
544        if ( binExp <= maxSmallBinExp && binExp >= minSmallBinExp ){
545            // Look more closely at the number to decide if,
546            // with scaling by 10^nTinyBits, the result will fit in
547            // a long.
548            if ( (nTinyBits < long5pow.length) && ((nFractBits + n5bits[nTinyBits]) < 64 ) ){
549                /*
550                 * We can do this:
551                 * take the fraction bits, which are normalized.
552                 * (a) nTinyBits == 0: Shift left or right appropriately
553                 *     to align the binary point at the extreme right, i.e.
554                 *     where a long int point is expected to be. The integer
555                 *     result is easily converted to a string.
556                 * (b) nTinyBits > 0: Shift right by expShift-nFractBits,
557                 *     which effectively converts to long and scales by
558                 *     2^nTinyBits. Then multiply by 5^nTinyBits to
559                 *     complete the scaling. We know this won't overflow
560                 *     because we just counted the number of bits necessary
561                 *     in the result. The integer you get from this can
562                 *     then be converted to a string pretty easily.
563                 */
564                long halfULP;
565                if ( nTinyBits == 0 ) {
566                    if ( binExp > nSignificantBits ){
567                        halfULP = 1L << ( binExp-nSignificantBits-1);
568                    } else {
569                        halfULP = 0L;
570                    }
571                    if ( binExp >= expShift ){
572                        fractBits <<= (binExp-expShift);
573                    } else {
574                        fractBits >>>= (expShift-binExp) ;
575                    }
576                    developLongDigits( 0, fractBits, halfULP );
577                    return;
578                }
579                /*
580                 * The following causes excess digits to be printed
581                 * out in the single-float case. Our manipulation of
582                 * halfULP here is apparently not correct. If we
583                 * better understand how this works, perhaps we can
584                 * use this special case again. But for the time being,
585                 * we do not.
586                 * else {
587                 *     fractBits >>>= expShift+1-nFractBits;
588                 *     fractBits *= long5pow[ nTinyBits ];
589                 *     halfULP = long5pow[ nTinyBits ] >> (1+nSignificantBits-nFractBits);
590                 *     developLongDigits( -nTinyBits, fractBits, halfULP );
591                 *     return;
592                 * }
593                 */
594            }
595        }
596        /*
597         * This is the hard case. We are going to compute large positive
598         * integers B and S and integer decExp, s.t.
599         *      d = ( B / S ) * 10^decExp
600         *      1 <= B / S < 10
601         * Obvious choices are:
602         *      decExp = floor( log10(d) )
603         *      B      = d * 2^nTinyBits * 10^max( 0, -decExp )
604         *      S      = 10^max( 0, decExp) * 2^nTinyBits
605         * (noting that nTinyBits has already been forced to non-negative)
606         * I am also going to compute a large positive integer
607         *      M      = (1/2^nSignificantBits) * 2^nTinyBits * 10^max( 0, -decExp )
608         * i.e. M is (1/2) of the ULP of d, scaled like B.
609         * When we iterate through dividing B/S and picking off the
610         * quotient bits, we will know when to stop when the remainder
611         * is <= M.
612         *
613         * We keep track of powers of 2 and powers of 5.
614         */
615
616        /*
617         * Estimate decimal exponent. (If it is small-ish,
618         * we could double-check.)
619         *
620         * First, scale the mantissa bits such that 1 <= d2 < 2.
621         * We are then going to estimate
622         *          log10(d2) ~=~  (d2-1.5)/1.5 + log(1.5)
623         * and so we can estimate
624         *      log10(d) ~=~ log10(d2) + binExp * log10(2)
625         * take the floor and call it decExp.
626         * FIXME -- use more precise constants here. It costs no more.
627         */
628        double d2 = Double.longBitsToDouble(
629            expOne | ( fractBits &~ fractHOB ) );
630        decExp = (int)Math.floor(
631            (d2-1.5D)*0.289529654D + 0.176091259 + (double)binExp * 0.301029995663981 );
632        int B2, B5; // powers of 2 and powers of 5, respectively, in B
633        int S2, S5; // powers of 2 and powers of 5, respectively, in S
634        int M2, M5; // powers of 2 and powers of 5, respectively, in M
635        int Bbits; // binary digits needed to represent B, approx.
636        int tenSbits; // binary digits needed to represent 10*S, approx.
637        FDBigInt Sval, Bval, Mval;
638
639        B5 = Math.max( 0, -decExp );
640        B2 = B5 + nTinyBits + binExp;
641
642        S5 = Math.max( 0, decExp );
643        S2 = S5 + nTinyBits;
644
645        M5 = B5;
646        M2 = B2 - nSignificantBits;
647
648        /*
649         * the long integer fractBits contains the (nFractBits) interesting
650         * bits from the mantissa of d ( hidden 1 added if necessary) followed
651         * by (expShift+1-nFractBits) zeros. In the interest of compactness,
652         * I will shift out those zeros before turning fractBits into a
653         * FDBigInt. The resulting whole number will be
654         *      d * 2^(nFractBits-1-binExp).
655         */
656        fractBits >>>= (expShift+1-nFractBits);
657        B2 -= nFractBits-1;
658        int common2factor = Math.min( B2, S2 );
659        B2 -= common2factor;
660        S2 -= common2factor;
661        M2 -= common2factor;
662
663        /*
664         * HACK!! For exact powers of two, the next smallest number
665         * is only half as far away as we think (because the meaning of
666         * ULP changes at power-of-two bounds) for this reason, we
667         * hack M2. Hope this works.
668         */
669        if ( nFractBits == 1 )
670            M2 -= 1;
671
672        if ( M2 < 0 ){
673            // oops.
674            // since we cannot scale M down far enough,
675            // we must scale the other values up.
676            B2 -= M2;
677            S2 -= M2;
678            M2 =  0;
679        }
680        /*
681         * Construct, Scale, iterate.
682         * Some day, we'll write a stopping test that takes
683         * account of the asymmetry of the spacing of floating-point
684         * numbers below perfect powers of 2
685         * 26 Sept 96 is not that day.
686         * So we use a symmetric test.
687         */
688        char digits[] = this.digits = new char[18];
689        int  ndigit = 0;
690        boolean low, high;
691        long lowDigitDifference;
692        int  q;
693
694        /*
695         * Detect the special cases where all the numbers we are about
696         * to compute will fit in int or long integers.
697         * In these cases, we will avoid doing FDBigInt arithmetic.
698         * We use the same algorithms, except that we "normalize"
699         * our FDBigInts before iterating. This is to make division easier,
700         * as it makes our fist guess (quotient of high-order words)
701         * more accurate!
702         *
703         * Some day, we'll write a stopping test that takes
704         * account of the asymmetry of the spacing of floating-point
705         * numbers below perfect powers of 2
706         * 26 Sept 96 is not that day.
707         * So we use a symmetric test.
708         */
709        Bbits = nFractBits + B2 + (( B5 < n5bits.length )? n5bits[B5] : ( B5*3 ));
710        tenSbits = S2+1 + (( (S5+1) < n5bits.length )? n5bits[(S5+1)] : ( (S5+1)*3 ));
711        if ( Bbits < 64 && tenSbits < 64){
712            if ( Bbits < 32 && tenSbits < 32){
713                // wa-hoo! They're all ints!
714                int b = ((int)fractBits * small5pow[B5] ) << B2;
715                int s = small5pow[S5] << S2;
716                int m = small5pow[M5] << M2;
717                int tens = s * 10;
718                /*
719                 * Unroll the first iteration. If our decExp estimate
720                 * was too high, our first quotient will be zero. In this
721                 * case, we discard it and decrement decExp.
722                 */
723                ndigit = 0;
724                q = b / s;
725                b = 10 * ( b % s );
726                m *= 10;
727                low  = (b <  m );
728                high = (b+m > tens );
729                assert q < 10 : q; // excessively large digit
730                if ( (q == 0) && ! high ){
731                    // oops. Usually ignore leading zero.
732                    decExp--;
733                } else {
734                    digits[ndigit++] = (char)('0' + q);
735                }
736                /*
737                 * HACK! Java spec sez that we always have at least
738                 * one digit after the . in either F- or E-form output.
739                 * Thus we will need more than one digit if we're using
740                 * E-form
741                 */
742                if ( decExp < -3 || decExp >= 8 ){
743                    high = low = false;
744                }
745                while( ! low && ! high ){
746                    q = b / s;
747                    b = 10 * ( b % s );
748                    m *= 10;
749                    assert q < 10 : q; // excessively large digit
750                    if ( m > 0L ){
751                        low  = (b <  m );
752                        high = (b+m > tens );
753                    } else {
754                        // hack -- m might overflow!
755                        // in this case, it is certainly > b,
756                        // which won't
757                        // and b+m > tens, too, since that has overflowed
758                        // either!
759                        low = true;
760                        high = true;
761                    }
762                    digits[ndigit++] = (char)('0' + q);
763                }
764                lowDigitDifference = (b<<1) - tens;
765            } else {
766                // still good! they're all longs!
767                long b = (fractBits * long5pow[B5] ) << B2;
768                long s = long5pow[S5] << S2;
769                long m = long5pow[M5] << M2;
770                long tens = s * 10L;
771                /*
772                 * Unroll the first iteration. If our decExp estimate
773                 * was too high, our first quotient will be zero. In this
774                 * case, we discard it and decrement decExp.
775                 */
776                ndigit = 0;
777                q = (int) ( b / s );
778                b = 10L * ( b % s );
779                m *= 10L;
780                low  = (b <  m );
781                high = (b+m > tens );
782                assert q < 10 : q; // excessively large digit
783                if ( (q == 0) && ! high ){
784                    // oops. Usually ignore leading zero.
785                    decExp--;
786                } else {
787                    digits[ndigit++] = (char)('0' + q);
788                }
789                /*
790                 * HACK! Java spec sez that we always have at least
791                 * one digit after the . in either F- or E-form output.
792                 * Thus we will need more than one digit if we're using
793                 * E-form
794                 */
795                if ( decExp < -3 || decExp >= 8 ){
796                    high = low = false;
797                }
798                while( ! low && ! high ){
799                    q = (int) ( b / s );
800                    b = 10 * ( b % s );
801                    m *= 10;
802                    assert q < 10 : q;  // excessively large digit
803                    if ( m > 0L ){
804                        low  = (b <  m );
805                        high = (b+m > tens );
806                    } else {
807                        // hack -- m might overflow!
808                        // in this case, it is certainly > b,
809                        // which won't
810                        // and b+m > tens, too, since that has overflowed
811                        // either!
812                        low = true;
813                        high = true;
814                    }
815                    digits[ndigit++] = (char)('0' + q);
816                }
817                lowDigitDifference = (b<<1) - tens;
818            }
819        } else {
820            FDBigInt tenSval;
821            int  shiftBias;
822
823            /*
824             * We really must do FDBigInt arithmetic.
825             * Fist, construct our FDBigInt initial values.
826             */
827            Bval = multPow52( new FDBigInt( fractBits  ), B5, B2 );
828            Sval = constructPow52( S5, S2 );
829            Mval = constructPow52( M5, M2 );
830
831
832            // normalize so that division works better
833            Bval.lshiftMe( shiftBias = Sval.normalizeMe() );
834            Mval.lshiftMe( shiftBias );
835            tenSval = Sval.mult( 10 );
836            /*
837             * Unroll the first iteration. If our decExp estimate
838             * was too high, our first quotient will be zero. In this
839             * case, we discard it and decrement decExp.
840             */
841            ndigit = 0;
842            q = Bval.quoRemIteration( Sval );
843            Mval = Mval.mult( 10 );
844            low  = (Bval.cmp( Mval ) < 0);
845            high = (Bval.add( Mval ).cmp( tenSval ) > 0 );
846            assert q < 10 : q; // excessively large digit
847            if ( (q == 0) && ! high ){
848                // oops. Usually ignore leading zero.
849                decExp--;
850            } else {
851                digits[ndigit++] = (char)('0' + q);
852            }
853            /*
854             * HACK! Java spec sez that we always have at least
855             * one digit after the . in either F- or E-form output.
856             * Thus we will need more than one digit if we're using
857             * E-form
858             */
859            if ( decExp < -3 || decExp >= 8 ){
860                high = low = false;
861            }
862            while( ! low && ! high ){
863                q = Bval.quoRemIteration( Sval );
864                Mval = Mval.mult( 10 );
865                assert q < 10 : q;  // excessively large digit
866                low  = (Bval.cmp( Mval ) < 0);
867                high = (Bval.add( Mval ).cmp( tenSval ) > 0 );
868                digits[ndigit++] = (char)('0' + q);
869            }
870            if ( high && low ){
871                Bval.lshiftMe(1);
872                lowDigitDifference = Bval.cmp(tenSval);
873            } else
874                lowDigitDifference = 0L; // this here only for flow analysis!
875        }
876        this.decExponent = decExp+1;
877        this.digits = digits;
878        this.nDigits = ndigit;
879        /*
880         * Last digit gets rounded based on stopping condition.
881         */
882        if ( high ){
883            if ( low ){
884                if ( lowDigitDifference == 0L ){
885                    // it's a tie!
886                    // choose based on which digits we like.
887                    if ( (digits[nDigits-1]&1) != 0 ) roundup();
888                } else if ( lowDigitDifference > 0 ){
889                    roundup();
890                }
891            } else {
892                roundup();
893            }
894        }
895    }
896
897    public String
898    toString(){
899        // most brain-dead version
900        StringBuffer result = new StringBuffer( nDigits+8 );
901        if ( isNegative ){ result.append( '-' ); }
902        if ( isExceptional ){
903            result.append( digits, 0, nDigits );
904        } else {
905            result.append( "0.");
906            result.append( digits, 0, nDigits );
907            result.append('e');
908            result.append( decExponent );
909        }
910        return new String(result);
911    }
912
913    public String toJavaFormatString() {
914        char result[] = (char[])(perThreadBuffer.get());
915        int i = getChars(result);
916        return new String(result, 0, i);
917    }
918
919    private int getChars(char[] result) {
920        assert nDigits <= 19 : nDigits; // generous bound on size of nDigits
921        int i = 0;
922        if (isNegative) { result[0] = '-'; i = 1; }
923        if (isExceptional) {
924            System.arraycopy(digits, 0, result, i, nDigits);
925            i += nDigits;
926        } else {
927            if (decExponent > 0 && decExponent < 8) {
928                // case with digits.digits
929                int charLength = Math.min(nDigits, decExponent);
930                System.arraycopy(digits, 0, result, i, charLength);
931                i += charLength;
932                if (charLength < decExponent) {
933                    charLength = decExponent-charLength;
934                    System.arraycopy(zero, 0, result, i, charLength);
935                    i += charLength;
936                    result[i++] = '.';
937                    result[i++] = '0';
938                } else {
939                    result[i++] = '.';
940                    if (charLength < nDigits) {
941                        int t = nDigits - charLength;
942                        System.arraycopy(digits, charLength, result, i, t);
943                        i += t;
944                    } else {
945                        result[i++] = '0';
946                    }
947                }
948            } else if (decExponent <=0 && decExponent > -3) {
949                // case with 0.digits
950                result[i++] = '0';
951                result[i++] = '.';
952                if (decExponent != 0) {
953                    System.arraycopy(zero, 0, result, i, -decExponent);
954                    i -= decExponent;
955                }
956                System.arraycopy(digits, 0, result, i, nDigits);
957                i += nDigits;
958            } else {
959                // case with digit.digitsEexponent
960                result[i++] = digits[0];
961                result[i++] = '.';
962                if (nDigits > 1) {
963                    System.arraycopy(digits, 1, result, i, nDigits-1);
964                    i += nDigits-1;
965                } else {
966                    result[i++] = '0';
967                }
968                result[i++] = 'E';
969                int e;
970                if (decExponent <= 0) {
971                    result[i++] = '-';
972                    e = -decExponent+1;
973                } else {
974                    e = decExponent-1;
975                }
976                // decExponent has 1, 2, or 3, digits
977                if (e <= 9) {
978                    result[i++] = (char)(e+'0');
979                } else if (e <= 99) {
980                    result[i++] = (char)(e/10 +'0');
981                    result[i++] = (char)(e%10 + '0');
982                } else {
983                    result[i++] = (char)(e/100+'0');
984                    e %= 100;
985                    result[i++] = (char)(e/10+'0');
986                    result[i++] = (char)(e%10 + '0');
987                }
988            }
989        }
990        return i;
991    }
992
993    // Per-thread buffer for string/stringbuffer conversion
994    private static ThreadLocal perThreadBuffer = new ThreadLocal() {
995            protected synchronized Object initialValue() {
996                return new char[26];
997            }
998        };
999
1000    public void appendTo(AbstractStringBuilder buf) {
1001        if (isNegative) { buf.append('-'); }
1002        if (isExceptional) {
1003            buf.append(digits, 0 , nDigits);
1004            return;
1005        }
1006        if (decExponent > 0 && decExponent < 8) {
1007            // print digits.digits.
1008            int charLength = Math.min(nDigits, decExponent);
1009            buf.append(digits, 0 , charLength);
1010            if (charLength < decExponent) {
1011                charLength = decExponent-charLength;
1012                buf.append(zero, 0 , charLength);
1013                buf.append(".0");
1014            } else {
1015                buf.append('.');
1016                if (charLength < nDigits) {
1017                    buf.append(digits, charLength, nDigits - charLength);
1018                } else {
1019                    buf.append('0');
1020                }
1021            }
1022        } else if (decExponent <=0 && decExponent > -3) {
1023            buf.append("0.");
1024            if (decExponent != 0) {
1025                buf.append(zero, 0, -decExponent);
1026            }
1027            buf.append(digits, 0, nDigits);
1028        } else {
1029            buf.append(digits[0]);
1030            buf.append('.');
1031            if (nDigits > 1) {
1032                buf.append(digits, 1, nDigits-1);
1033            } else {
1034                buf.append('0');
1035            }
1036            buf.append('E');
1037            int e;
1038            if (decExponent <= 0) {
1039                buf.append('-');
1040                e = -decExponent + 1;
1041            } else {
1042                e = decExponent - 1;
1043            }
1044            // decExponent has 1, 2, or 3, digits
1045            if (e <= 9) {
1046                buf.append((char)(e + '0'));
1047            } else if (e <= 99) {
1048                buf.append((char)(e/10 + '0'));
1049                buf.append((char)(e%10 + '0'));
1050            } else {
1051                buf.append((char)(e/100 + '0'));
1052                e %= 100;
1053                buf.append((char)(e/10 + '0'));
1054                buf.append((char)(e%10 + '0'));
1055            }
1056        }
1057    }
1058
1059    public FloatingDecimal
1060    readJavaFormatString( String in ) throws NumberFormatException {
1061        boolean isNegative = false;
1062        boolean signSeen   = false;
1063        int     decExp;
1064        char    c;
1065
1066    parseNumber:
1067        try{
1068            in = in.trim(); // don't fool around with white space.
1069                            // throws NullPointerException if null
1070            int l = in.length();
1071            if ( l == 0 ) throw new NumberFormatException("empty String");
1072            int i = 0;
1073            switch ( c = in.charAt( i ) ){
1074            case '-':
1075                isNegative = true;
1076                //FALLTHROUGH
1077            case '+':
1078                i++;
1079                signSeen = true;
1080            }
1081
1082            // Check for NaN and Infinity strings
1083            c = in.charAt(i);
1084            if(c == 'N' || c == 'I') { // possible NaN or infinity
1085                boolean potentialNaN = false;
1086                char targetChars[] = null;  // char array of "NaN" or "Infinity"
1087
1088                if(c == 'N') {
1089                    targetChars = notANumber;
1090                    potentialNaN = true;
1091                } else {
1092                    targetChars = infinity;
1093                }
1094
1095                // compare Input string to "NaN" or "Infinity"
1096                int j = 0;
1097                while(i < l && j < targetChars.length) {
1098                    if(in.charAt(i) == targetChars[j]) {
1099                        i++; j++;
1100                    }
1101                    else // something is amiss, throw exception
1102                        break parseNumber;
1103                }
1104
1105                // For the candidate string to be a NaN or infinity,
1106                // all characters in input string and target char[]
1107                // must be matched ==> j must equal targetChars.length
1108                // and i must equal l
1109                if( (j == targetChars.length) && (i == l) ) { // return NaN or infinity
1110                    return (potentialNaN ? loadDouble(Double.NaN) // NaN has no sign
1111                            : loadDouble(isNegative?
1112                                                  Double.NEGATIVE_INFINITY:
1113                                                  Double.POSITIVE_INFINITY)) ;
1114                }
1115                else { // something went wrong, throw exception
1116                    break parseNumber;
1117                }
1118
1119            } else if (c == '0')  { // check for hexadecimal floating-point number
1120                if (l > i+1 ) {
1121                    char ch = in.charAt(i+1);
1122                    if (ch == 'x' || ch == 'X' ) // possible hex string
1123                        return parseHexString(in);
1124                }
1125            }  // look for and process decimal floating-point string
1126
1127            char[] digits = new char[ l ];
1128            int    nDigits= 0;
1129            boolean decSeen = false;
1130            int decPt = 0;
1131            int nLeadZero = 0;
1132            int nTrailZero= 0;
1133        digitLoop:
1134            while ( i < l ){
1135                switch ( c = in.charAt( i ) ){
1136                case '0':
1137                    if ( nDigits > 0 ){
1138                        nTrailZero += 1;
1139                    } else {
1140                        nLeadZero += 1;
1141                    }
1142                    break; // out of switch.
1143                case '1':
1144                case '2':
1145                case '3':
1146                case '4':
1147                case '5':
1148                case '6':
1149                case '7':
1150                case '8':
1151                case '9':
1152                    while ( nTrailZero > 0 ){
1153                        digits[nDigits++] = '0';
1154                        nTrailZero -= 1;
1155                    }
1156                    digits[nDigits++] = c;
1157                    break; // out of switch.
1158                case '.':
1159                    if ( decSeen ){
1160                        // already saw one ., this is the 2nd.
1161                        throw new NumberFormatException("multiple points");
1162                    }
1163                    decPt = i;
1164                    if ( signSeen ){
1165                        decPt -= 1;
1166                    }
1167                    decSeen = true;
1168                    break; // out of switch.
1169                default:
1170                    break digitLoop;
1171                }
1172                i++;
1173            }
1174            /*
1175             * At this point, we've scanned all the digits and decimal
1176             * point we're going to see. Trim off leading and trailing
1177             * zeros, which will just confuse us later, and adjust
1178             * our initial decimal exponent accordingly.
1179             * To review:
1180             * we have seen i total characters.
1181             * nLeadZero of them were zeros before any other digits.
1182             * nTrailZero of them were zeros after any other digits.
1183             * if ( decSeen ), then a . was seen after decPt characters
1184             * ( including leading zeros which have been discarded )
1185             * nDigits characters were neither lead nor trailing
1186             * zeros, nor point
1187             */
1188            /*
1189             * special hack: if we saw no non-zero digits, then the
1190             * answer is zero!
1191             * Unfortunately, we feel honor-bound to keep parsing!
1192             */
1193            if ( nDigits == 0 ){
1194                digits = zero;
1195                nDigits = 1;
1196                if ( nLeadZero == 0 ){
1197                    // we saw NO DIGITS AT ALL,
1198                    // not even a crummy 0!
1199                    // this is not allowed.
1200                    break parseNumber; // go throw exception
1201                }
1202
1203            }
1204
1205            /* Our initial exponent is decPt, adjusted by the number of
1206             * discarded zeros. Or, if there was no decPt,
1207             * then its just nDigits adjusted by discarded trailing zeros.
1208             */
1209            if ( decSeen ){
1210                decExp = decPt - nLeadZero;
1211            } else {
1212                decExp = nDigits+nTrailZero;
1213            }
1214
1215            /*
1216             * Look for 'e' or 'E' and an optionally signed integer.
1217             */
1218            if ( (i < l) &&  (((c = in.charAt(i) )=='e') || (c == 'E') ) ){
1219                int expSign = 1;
1220                int expVal  = 0;
1221                int reallyBig = Integer.MAX_VALUE / 10;
1222                boolean expOverflow = false;
1223                switch( in.charAt(++i) ){
1224                case '-':
1225                    expSign = -1;
1226                    //FALLTHROUGH
1227                case '+':
1228                    i++;
1229                }
1230                int expAt = i;
1231            expLoop:
1232                while ( i < l  ){
1233                    if ( expVal >= reallyBig ){
1234                        // the next character will cause integer
1235                        // overflow.
1236                        expOverflow = true;
1237                    }
1238                    switch ( c = in.charAt(i++) ){
1239                    case '0':
1240                    case '1':
1241                    case '2':
1242                    case '3':
1243                    case '4':
1244                    case '5':
1245                    case '6':
1246                    case '7':
1247                    case '8':
1248                    case '9':
1249                        expVal = expVal*10 + ( (int)c - (int)'0' );
1250                        continue;
1251                    default:
1252                        i--;           // back up.
1253                        break expLoop; // stop parsing exponent.
1254                    }
1255                }
1256                int expLimit = bigDecimalExponent+nDigits+nTrailZero;
1257                if ( expOverflow || ( expVal > expLimit ) ){
1258                    //
1259                    // The intent here is to end up with
1260                    // infinity or zero, as appropriate.
1261                    // The reason for yielding such a small decExponent,
1262                    // rather than something intuitive such as
1263                    // expSign*Integer.MAX_VALUE, is that this value
1264                    // is subject to further manipulation in
1265                    // doubleValue() and floatValue(), and I don't want
1266                    // it to be able to cause overflow there!
1267                    // (The only way we can get into trouble here is for
1268                    // really outrageous nDigits+nTrailZero, such as 2 billion. )
1269                    //
1270                    decExp = expSign*expLimit;
1271                } else {
1272                    // this should not overflow, since we tested
1273                    // for expVal > (MAX+N), where N >= abs(decExp)
1274                    decExp = decExp + expSign*expVal;
1275                }
1276
1277                // if we saw something not a digit ( or end of string )
1278                // after the [Ee][+-], without seeing any digits at all
1279                // this is certainly an error. If we saw some digits,
1280                // but then some trailing garbage, that might be ok.
1281                // so we just fall through in that case.
1282                // HUMBUG
1283                if ( i == expAt )
1284                    break parseNumber; // certainly bad
1285            }
1286            /*
1287             * We parsed everything we could.
1288             * If there are leftovers, then this is not good input!
1289             */
1290            if ( i < l &&
1291                ((i != l - 1) ||
1292                (in.charAt(i) != 'f' &&
1293                 in.charAt(i) != 'F' &&
1294                 in.charAt(i) != 'd' &&
1295                 in.charAt(i) != 'D'))) {
1296                break parseNumber; // go throw exception
1297            }
1298
1299            this.isNegative = isNegative;
1300            this.decExponent = decExp;
1301            this.digits = digits;
1302            this.nDigits = nDigits;
1303            this.isExceptional = false;
1304            return this;
1305        } catch ( StringIndexOutOfBoundsException e ){ }
1306        throw new NumberFormatException("For input string: \"" + in + "\"");
1307    }
1308
1309    /*
1310     * Take a FloatingDecimal, which we presumably just scanned in,
1311     * and find out what its value is, as a double.
1312     *
1313     * AS A SIDE EFFECT, SET roundDir TO INDICATE PREFERRED
1314     * ROUNDING DIRECTION in case the result is really destined
1315     * for a single-precision float.
1316     */
1317
1318    public strictfp double doubleValue(){
1319        int     kDigits = Math.min( nDigits, maxDecimalDigits+1 );
1320        long    lValue;
1321        double  dValue;
1322        double  rValue, tValue;
1323
1324        // First, check for NaN and Infinity values
1325        if(digits == infinity || digits == notANumber) {
1326            if(digits == notANumber)
1327                return Double.NaN;
1328            else
1329                return (isNegative?Double.NEGATIVE_INFINITY:Double.POSITIVE_INFINITY);
1330        }
1331        else {
1332            if (mustSetRoundDir) {
1333                roundDir = 0;
1334            }
1335            /*
1336             * convert the lead kDigits to a long integer.
1337             */
1338            // (special performance hack: start to do it using int)
1339            int iValue = (int)digits[0]-(int)'0';
1340            int iDigits = Math.min( kDigits, intDecimalDigits );
1341            for ( int i=1; i < iDigits; i++ ){
1342                iValue = iValue*10 + (int)digits[i]-(int)'0';
1343            }
1344            lValue = (long)iValue;
1345            for ( int i=iDigits; i < kDigits; i++ ){
1346                lValue = lValue*10L + (long)((int)digits[i]-(int)'0');
1347            }
1348            dValue = (double)lValue;
1349            int exp = decExponent-kDigits;
1350            /*
1351             * lValue now contains a long integer with the value of
1352             * the first kDigits digits of the number.
1353             * dValue contains the (double) of the same.
1354             */
1355
1356            if ( nDigits <= maxDecimalDigits ){
1357                /*
1358                 * possibly an easy case.
1359                 * We know that the digits can be represented
1360                 * exactly. And if the exponent isn't too outrageous,
1361                 * the whole thing can be done with one operation,
1362                 * thus one rounding error.
1363                 * Note that all our constructors trim all leading and
1364                 * trailing zeros, so simple values (including zero)
1365                 * will always end up here
1366                 */
1367                if (exp == 0 || dValue == 0.0)
1368                    return (isNegative)? -dValue : dValue; // small floating integer
1369                else if ( exp >= 0 ){
1370                    if ( exp <= maxSmallTen ){
1371                        /*
1372                         * Can get the answer with one operation,
1373                         * thus one roundoff.
1374                         */
1375                        rValue = dValue * small10pow[exp];
1376                        if ( mustSetRoundDir ){
1377                            tValue = rValue / small10pow[exp];
1378                            roundDir = ( tValue ==  dValue ) ? 0
1379                                :( tValue < dValue ) ? 1
1380                                : -1;
1381                        }
1382                        return (isNegative)? -rValue : rValue;
1383                    }
1384                    int slop = maxDecimalDigits - kDigits;
1385                    if ( exp <= maxSmallTen+slop ){
1386                        /*
1387                         * We can multiply dValue by 10^(slop)
1388                         * and it is still "small" and exact.
1389                         * Then we can multiply by 10^(exp-slop)
1390                         * with one rounding.
1391                         */
1392                        dValue *= small10pow[slop];
1393                        rValue = dValue * small10pow[exp-slop];
1394
1395                        if ( mustSetRoundDir ){
1396                            tValue = rValue / small10pow[exp-slop];
1397                            roundDir = ( tValue ==  dValue ) ? 0
1398                                :( tValue < dValue ) ? 1
1399                                : -1;
1400                        }
1401                        return (isNegative)? -rValue : rValue;
1402                    }
1403                    /*
1404                     * Else we have a hard case with a positive exp.
1405                     */
1406                } else {
1407                    if ( exp >= -maxSmallTen ){
1408                        /*
1409                         * Can get the answer in one division.
1410                         */
1411                        rValue = dValue / small10pow[-exp];
1412                        tValue = rValue * small10pow[-exp];
1413                        if ( mustSetRoundDir ){
1414                            roundDir = ( tValue ==  dValue ) ? 0
1415                                :( tValue < dValue ) ? 1
1416                                : -1;
1417                        }
1418                        return (isNegative)? -rValue : rValue;
1419                    }
1420                    /*
1421                     * Else we have a hard case with a negative exp.
1422                     */
1423                }
1424            }
1425
1426            /*
1427             * Harder cases:
1428             * The sum of digits plus exponent is greater than
1429             * what we think we can do with one error.
1430             *
1431             * Start by approximating the right answer by,
1432             * naively, scaling by powers of 10.
1433             */
1434            if ( exp > 0 ){
1435                if ( decExponent > maxDecimalExponent+1 ){
1436                    /*
1437                     * Lets face it. This is going to be
1438                     * Infinity. Cut to the chase.
1439                     */
1440                    return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
1441                }
1442                if ( (exp&15) != 0 ){
1443                    dValue *= small10pow[exp&15];
1444                }
1445                if ( (exp>>=4) != 0 ){
1446                    int j;
1447                    for( j = 0; exp > 1; j++, exp>>=1 ){
1448                        if ( (exp&1)!=0)
1449                            dValue *= big10pow[j];
1450                    }
1451                    /*
1452                     * The reason for the weird exp > 1 condition
1453                     * in the above loop was so that the last multiply
1454                     * would get unrolled. We handle it here.
1455                     * It could overflow.
1456                     */
1457                    double t = dValue * big10pow[j];
1458                    if ( Double.isInfinite( t ) ){
1459                        /*
1460                         * It did overflow.
1461                         * Look more closely at the result.
1462                         * If the exponent is just one too large,
1463                         * then use the maximum finite as our estimate
1464                         * value. Else call the result infinity
1465                         * and punt it.
1466                         * ( I presume this could happen because
1467                         * rounding forces the result here to be
1468                         * an ULP or two larger than
1469                         * Double.MAX_VALUE ).
1470                         */
1471                        t = dValue / 2.0;
1472                        t *= big10pow[j];
1473                        if ( Double.isInfinite( t ) ){
1474                            return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
1475                        }
1476                        t = Double.MAX_VALUE;
1477                    }
1478                    dValue = t;
1479                }
1480            } else if ( exp < 0 ){
1481                exp = -exp;
1482                if ( decExponent < minDecimalExponent-1 ){
1483                    /*
1484                     * Lets face it. This is going to be
1485                     * zero. Cut to the chase.
1486                     */
1487                    return (isNegative)? -0.0 : 0.0;
1488                }
1489                if ( (exp&15) != 0 ){
1490                    dValue /= small10pow[exp&15];
1491                }
1492                if ( (exp>>=4) != 0 ){
1493                    int j;
1494                    for( j = 0; exp > 1; j++, exp>>=1 ){
1495                        if ( (exp&1)!=0)
1496                            dValue *= tiny10pow[j];
1497                    }
1498                    /*
1499                     * The reason for the weird exp > 1 condition
1500                     * in the above loop was so that the last multiply
1501                     * would get unrolled. We handle it here.
1502                     * It could underflow.
1503                     */
1504                    double t = dValue * tiny10pow[j];
1505                    if ( t == 0.0 ){
1506                        /*
1507                         * It did underflow.
1508                         * Look more closely at the result.
1509                         * If the exponent is just one too small,
1510                         * then use the minimum finite as our estimate
1511                         * value. Else call the result 0.0
1512                         * and punt it.
1513                         * ( I presume this could happen because
1514                         * rounding forces the result here to be
1515                         * an ULP or two less than
1516                         * Double.MIN_VALUE ).
1517                         */
1518                        t = dValue * 2.0;
1519                        t *= tiny10pow[j];
1520                        if ( t == 0.0 ){
1521                            return (isNegative)? -0.0 : 0.0;
1522                        }
1523                        t = Double.MIN_VALUE;
1524                    }
1525                    dValue = t;
1526                }
1527            }
1528
1529            /*
1530             * dValue is now approximately the result.
1531             * The hard part is adjusting it, by comparison
1532             * with FDBigInt arithmetic.
1533             * Formulate the EXACT big-number result as
1534             * bigD0 * 10^exp
1535             */
1536            FDBigInt bigD0 = new FDBigInt( lValue, digits, kDigits, nDigits );
1537            exp   = decExponent - nDigits;
1538
1539            correctionLoop:
1540            while(true){
1541                /* AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES
1542                 * bigIntExp and bigIntNBits
1543                 */
1544                FDBigInt bigB = doubleToBigInt( dValue );
1545
1546                /*
1547                 * Scale bigD, bigB appropriately for
1548                 * big-integer operations.
1549                 * Naively, we multiply by powers of ten
1550                 * and powers of two. What we actually do
1551                 * is keep track of the powers of 5 and
1552                 * powers of 2 we would use, then factor out
1553                 * common divisors before doing the work.
1554                 */
1555                int B2, B5; // powers of 2, 5 in bigB
1556                int     D2, D5; // powers of 2, 5 in bigD
1557                int Ulp2;   // powers of 2 in halfUlp.
1558                if ( exp >= 0 ){
1559                    B2 = B5 = 0;
1560                    D2 = D5 = exp;
1561                } else {
1562                    B2 = B5 = -exp;
1563                    D2 = D5 = 0;
1564                }
1565                if ( bigIntExp >= 0 ){
1566                    B2 += bigIntExp;
1567                } else {
1568                    D2 -= bigIntExp;
1569                }
1570                Ulp2 = B2;
1571                // shift bigB and bigD left by a number s. t.
1572                // halfUlp is still an integer.
1573                int hulpbias;
1574                if ( bigIntExp+bigIntNBits <= -expBias+1 ){
1575                    // This is going to be a denormalized number
1576                    // (if not actually zero).
1577                    // half an ULP is at 2^-(expBias+expShift+1)
1578                    hulpbias = bigIntExp+ expBias + expShift;
1579                } else {
1580                    hulpbias = expShift + 2 - bigIntNBits;
1581                }
1582                B2 += hulpbias;
1583                D2 += hulpbias;
1584                // if there are common factors of 2, we might just as well
1585                // factor them out, as they add nothing useful.
1586                int common2 = Math.min( B2, Math.min( D2, Ulp2 ) );
1587                B2 -= common2;
1588                D2 -= common2;
1589                Ulp2 -= common2;
1590                // do multiplications by powers of 5 and 2
1591                bigB = multPow52( bigB, B5, B2 );
1592                FDBigInt bigD = multPow52( new FDBigInt( bigD0 ), D5, D2 );
1593                //
1594                // to recap:
1595                // bigB is the scaled-big-int version of our floating-point
1596                // candidate.
1597                // bigD is the scaled-big-int version of the exact value
1598                // as we understand it.
1599                // halfUlp is 1/2 an ulp of bigB, except for special cases
1600                // of exact powers of 2
1601                //
1602                // the plan is to compare bigB with bigD, and if the difference
1603                // is less than halfUlp, then we're satisfied. Otherwise,
1604                // use the ratio of difference to halfUlp to calculate a fudge
1605                // factor to add to the floating value, then go 'round again.
1606                //
1607                FDBigInt diff;
1608                int cmpResult;
1609                boolean overvalue;
1610                if ( (cmpResult = bigB.cmp( bigD ) ) > 0 ){
1611                    overvalue = true; // our candidate is too big.
1612                    diff = bigB.sub( bigD );
1613                    if ( (bigIntNBits == 1) && (bigIntExp > -expBias+1) ){
1614                        // candidate is a normalized exact power of 2 and
1615                        // is too big. We will be subtracting.
1616                        // For our purposes, ulp is the ulp of the
1617                        // next smaller range.
1618                        Ulp2 -= 1;
1619                        if ( Ulp2 < 0 ){
1620                            // rats. Cannot de-scale ulp this far.
1621                            // must scale diff in other direction.
1622                            Ulp2 = 0;
1623                            diff.lshiftMe( 1 );
1624                        }
1625                    }
1626                } else if ( cmpResult < 0 ){
1627                    overvalue = false; // our candidate is too small.
1628                    diff = bigD.sub( bigB );
1629                } else {
1630                    // the candidate is exactly right!
1631                    // this happens with surprising frequency
1632                    break correctionLoop;
1633                }
1634                FDBigInt halfUlp = constructPow52( B5, Ulp2 );
1635                if ( (cmpResult = diff.cmp( halfUlp ) ) < 0 ){
1636                    // difference is small.
1637                    // this is close enough
1638                    if (mustSetRoundDir) {
1639                        roundDir = overvalue ? -1 : 1;
1640                    }
1641                    break correctionLoop;
1642                } else if ( cmpResult == 0 ){
1643                    // difference is exactly half an ULP
1644                    // round to some other value maybe, then finish
1645                    dValue += 0.5*ulp( dValue, overvalue );
1646                    // should check for bigIntNBits == 1 here??
1647                    if (mustSetRoundDir) {
1648                        roundDir = overvalue ? -1 : 1;
1649                    }
1650                    break correctionLoop;
1651                } else {
1652                    // difference is non-trivial.
1653                    // could scale addend by ratio of difference to
1654                    // halfUlp here, if we bothered to compute that difference.
1655                    // Most of the time ( I hope ) it is about 1 anyway.
1656                    dValue += ulp( dValue, overvalue );
1657                    if ( dValue == 0.0 || dValue == Double.POSITIVE_INFINITY )
1658                        break correctionLoop; // oops. Fell off end of range.
1659                    continue; // try again.
1660                }
1661
1662            }
1663            return (isNegative)? -dValue : dValue;
1664        }
1665    }
1666
1667    /*
1668     * Take a FloatingDecimal, which we presumably just scanned in,
1669     * and find out what its value is, as a float.
1670     * This is distinct from doubleValue() to avoid the extremely
1671     * unlikely case of a double rounding error, wherein the conversion
1672     * to double has one rounding error, and the conversion of that double
1673     * to a float has another rounding error, IN THE WRONG DIRECTION,
1674     * ( because of the preference to a zero low-order bit ).
1675     */
1676
1677    public strictfp float floatValue(){
1678        int     kDigits = Math.min( nDigits, singleMaxDecimalDigits+1 );
1679        int     iValue;
1680        float   fValue;
1681
1682        // First, check for NaN and Infinity values
1683        if(digits == infinity || digits == notANumber) {
1684            if(digits == notANumber)
1685                return Float.NaN;
1686            else
1687                return (isNegative?Float.NEGATIVE_INFINITY:Float.POSITIVE_INFINITY);
1688        }
1689        else {
1690            /*
1691             * convert the lead kDigits to an integer.
1692             */
1693            iValue = (int)digits[0]-(int)'0';
1694            for ( int i=1; i < kDigits; i++ ){
1695                iValue = iValue*10 + (int)digits[i]-(int)'0';
1696            }
1697            fValue = (float)iValue;
1698            int exp = decExponent-kDigits;
1699            /*
1700             * iValue now contains an integer with the value of
1701             * the first kDigits digits of the number.
1702             * fValue contains the (float) of the same.
1703             */
1704
1705            if ( nDigits <= singleMaxDecimalDigits ){
1706                /*
1707                 * possibly an easy case.
1708                 * We know that the digits can be represented
1709                 * exactly. And if the exponent isn't too outrageous,
1710                 * the whole thing can be done with one operation,
1711                 * thus one rounding error.
1712                 * Note that all our constructors trim all leading and
1713                 * trailing zeros, so simple values (including zero)
1714                 * will always end up here.
1715                 */
1716                if (exp == 0 || fValue == 0.0f)
1717                    return (isNegative)? -fValue : fValue; // small floating integer
1718                else if ( exp >= 0 ){
1719                    if ( exp <= singleMaxSmallTen ){
1720                        /*
1721                         * Can get the answer with one operation,
1722                         * thus one roundoff.
1723                         */
1724                        fValue *= singleSmall10pow[exp];
1725                        return (isNegative)? -fValue : fValue;
1726                    }
1727                    int slop = singleMaxDecimalDigits - kDigits;
1728                    if ( exp <= singleMaxSmallTen+slop ){
1729                        /*
1730                         * We can multiply dValue by 10^(slop)
1731                         * and it is still "small" and exact.
1732                         * Then we can multiply by 10^(exp-slop)
1733                         * with one rounding.
1734                         */
1735                        fValue *= singleSmall10pow[slop];
1736                        fValue *= singleSmall10pow[exp-slop];
1737                        return (isNegative)? -fValue : fValue;
1738                    }
1739                    /*
1740                     * Else we have a hard case with a positive exp.
1741                     */
1742                } else {
1743                    if ( exp >= -singleMaxSmallTen ){
1744                        /*
1745                         * Can get the answer in one division.
1746                         */
1747                        fValue /= singleSmall10pow[-exp];
1748                        return (isNegative)? -fValue : fValue;
1749                    }
1750                    /*
1751                     * Else we have a hard case with a negative exp.
1752                     */
1753                }
1754            } else if ( (decExponent >= nDigits) && (nDigits+decExponent <= maxDecimalDigits) ){
1755                /*
1756                 * In double-precision, this is an exact floating integer.
1757                 * So we can compute to double, then shorten to float
1758                 * with one round, and get the right answer.
1759                 *
1760                 * First, finish accumulating digits.
1761                 * Then convert that integer to a double, multiply
1762                 * by the appropriate power of ten, and convert to float.
1763                 */
1764                long lValue = (long)iValue;
1765                for ( int i=kDigits; i < nDigits; i++ ){
1766                    lValue = lValue*10L + (long)((int)digits[i]-(int)'0');
1767                }
1768                double dValue = (double)lValue;
1769                exp = decExponent-nDigits;
1770                dValue *= small10pow[exp];
1771                fValue = (float)dValue;
1772                return (isNegative)? -fValue : fValue;
1773
1774            }
1775            /*
1776             * Harder cases:
1777             * The sum of digits plus exponent is greater than
1778             * what we think we can do with one error.
1779             *
1780             * Start by weeding out obviously out-of-range
1781             * results, then convert to double and go to
1782             * common hard-case code.
1783             */
1784            if ( decExponent > singleMaxDecimalExponent+1 ){
1785                /*
1786                 * Lets face it. This is going to be
1787                 * Infinity. Cut to the chase.
1788                 */
1789                return (isNegative)? Float.NEGATIVE_INFINITY : Float.POSITIVE_INFINITY;
1790            } else if ( decExponent < singleMinDecimalExponent-1 ){
1791                /*
1792                 * Lets face it. This is going to be
1793                 * zero. Cut to the chase.
1794                 */
1795                return (isNegative)? -0.0f : 0.0f;
1796            }
1797
1798            /*
1799             * Here, we do 'way too much work, but throwing away
1800             * our partial results, and going and doing the whole
1801             * thing as double, then throwing away half the bits that computes
1802             * when we convert back to float.
1803             *
1804             * The alternative is to reproduce the whole multiple-precision
1805             * algorithm for float precision, or to try to parameterize it
1806             * for common usage. The former will take about 400 lines of code,
1807             * and the latter I tried without success. Thus the semi-hack
1808             * answer here.
1809             */
1810            mustSetRoundDir = !fromHex;
1811            double dValue = doubleValue();
1812            return stickyRound( dValue );
1813        }
1814    }
1815
1816
1817    /*
1818     * All the positive powers of 10 that can be
1819     * represented exactly in double/float.
1820     */
1821    private static final double small10pow[] = {
1822        1.0e0,
1823        1.0e1, 1.0e2, 1.0e3, 1.0e4, 1.0e5,
1824        1.0e6, 1.0e7, 1.0e8, 1.0e9, 1.0e10,
1825        1.0e11, 1.0e12, 1.0e13, 1.0e14, 1.0e15,
1826        1.0e16, 1.0e17, 1.0e18, 1.0e19, 1.0e20,
1827        1.0e21, 1.0e22
1828    };
1829
1830    private static final float singleSmall10pow[] = {
1831        1.0e0f,
1832        1.0e1f, 1.0e2f, 1.0e3f, 1.0e4f, 1.0e5f,
1833        1.0e6f, 1.0e7f, 1.0e8f, 1.0e9f, 1.0e10f
1834    };
1835
1836    private static final double big10pow[] = {
1837        1e16, 1e32, 1e64, 1e128, 1e256 };
1838    private static final double tiny10pow[] = {
1839        1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1840
1841    private static final int maxSmallTen = small10pow.length-1;
1842    private static final int singleMaxSmallTen = singleSmall10pow.length-1;
1843
1844    private static final int small5pow[] = {
1845        1,
1846        5,
1847        5*5,
1848        5*5*5,
1849        5*5*5*5,
1850        5*5*5*5*5,
1851        5*5*5*5*5*5,
1852        5*5*5*5*5*5*5,
1853        5*5*5*5*5*5*5*5,
1854        5*5*5*5*5*5*5*5*5,
1855        5*5*5*5*5*5*5*5*5*5,
1856        5*5*5*5*5*5*5*5*5*5*5,
1857        5*5*5*5*5*5*5*5*5*5*5*5,
1858        5*5*5*5*5*5*5*5*5*5*5*5*5
1859    };
1860
1861
1862    private static final long long5pow[] = {
1863        1L,
1864        5L,
1865        5L*5,
1866        5L*5*5,
1867        5L*5*5*5,
1868        5L*5*5*5*5,
1869        5L*5*5*5*5*5,
1870        5L*5*5*5*5*5*5,
1871        5L*5*5*5*5*5*5*5,
1872        5L*5*5*5*5*5*5*5*5,
1873        5L*5*5*5*5*5*5*5*5*5,
1874        5L*5*5*5*5*5*5*5*5*5*5,
1875        5L*5*5*5*5*5*5*5*5*5*5*5,
1876        5L*5*5*5*5*5*5*5*5*5*5*5*5,
1877        5L*5*5*5*5*5*5*5*5*5*5*5*5*5,
1878        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1879        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1880        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1881        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1882        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1883        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1884        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1885        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1886        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1887        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1888        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1889        5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1890    };
1891
1892    // approximately ceil( log2( long5pow[i] ) )
1893    private static final int n5bits[] = {
1894        0,
1895        3,
1896        5,
1897        7,
1898        10,
1899        12,
1900        14,
1901        17,
1902        19,
1903        21,
1904        24,
1905        26,
1906        28,
1907        31,
1908        33,
1909        35,
1910        38,
1911        40,
1912        42,
1913        45,
1914        47,
1915        49,
1916        52,
1917        54,
1918        56,
1919        59,
1920        61,
1921    };
1922
1923    private static final char infinity[] = { 'I', 'n', 'f', 'i', 'n', 'i', 't', 'y' };
1924    private static final char notANumber[] = { 'N', 'a', 'N' };
1925    private static final char zero[] = { '0', '0', '0', '0', '0', '0', '0', '0' };
1926
1927
1928    /*
1929     * Grammar is compatible with hexadecimal floating-point constants
1930     * described in section 6.4.4.2 of the C99 specification.
1931     */
1932    private static Pattern hexFloatPattern = null;
1933    private static synchronized Pattern getHexFloatPattern() {
1934        if (hexFloatPattern == null) {
1935           hexFloatPattern = Pattern.compile(
1936                   //1           234                   56                7                   8      9
1937                    "([-+])?0[xX](((\\p{XDigit}+)\\.?)|((\\p{XDigit}*)\\.(\\p{XDigit}+)))[pP]([-+])?(\\p{Digit}+)[fFdD]?"
1938                    );
1939        }
1940        return hexFloatPattern;
1941    }
1942
1943    /*
1944     * Convert string s to a suitable floating decimal; uses the
1945     * double constructor and set the roundDir variable appropriately
1946     * in case the value is later converted to a float.
1947     */
1948    FloatingDecimal parseHexString(String s) {
1949        // Verify string is a member of the hexadecimal floating-point
1950        // string language.
1951        Matcher m = getHexFloatPattern().matcher(s);
1952        boolean validInput = m.matches();
1953
1954        if (!validInput) {
1955            // Input does not match pattern
1956            throw new NumberFormatException("For input string: \"" + s + "\"");
1957        } else { // validInput
1958            /*
1959             * We must isolate the sign, significand, and exponent
1960             * fields.  The sign value is straightforward.  Since
1961             * floating-point numbers are stored with a normalized
1962             * representation, the significand and exponent are
1963             * interrelated.
1964             *
1965             * After extracting the sign, we normalized the
1966             * significand as a hexadecimal value, calculating an
1967             * exponent adjust for any shifts made during
1968             * normalization.  If the significand is zero, the
1969             * exponent doesn't need to be examined since the output
1970             * will be zero.
1971             *
1972             * Next the exponent in the input string is extracted.
1973             * Afterwards, the significand is normalized as a *binary*
1974             * value and the input value's normalized exponent can be
1975             * computed.  The significand bits are copied into a
1976             * double significand; if the string has more logical bits
1977             * than can fit in a double, the extra bits affect the
1978             * round and sticky bits which are used to round the final
1979             * value.
1980             */
1981
1982            //  Extract significand sign
1983            String group1 = m.group(1);
1984            double sign = (( group1 == null ) || group1.equals("+"))? 1.0 : -1.0;
1985
1986
1987            //  Extract Significand magnitude
1988            /*
1989             * Based on the form of the significand, calculate how the
1990             * binary exponent needs to be adjusted to create a
1991             * normalized *hexadecimal* floating-point number; that
1992             * is, a number where there is one nonzero hex digit to
1993             * the left of the (hexa)decimal point.  Since we are
1994             * adjusting a binary, not hexadecimal exponent, the
1995             * exponent is adjusted by a multiple of 4.
1996             *
1997             * There are a number of significand scenarios to consider;
1998             * letters are used in indicate nonzero digits:
1999             *
2000             * 1. 000xxxx       =>      x.xxx   normalized
2001             *    increase exponent by (number of x's - 1)*4
2002             *
2003             * 2. 000xxx.yyyy =>        x.xxyyyy        normalized
2004             *    increase exponent by (number of x's - 1)*4
2005             *
2006             * 3. .000yyy  =>   y.yy    normalized
2007             *    decrease exponent by (number of zeros + 1)*4
2008             *
2009             * 4. 000.00000yyy => y.yy normalized
2010             *    decrease exponent by (number of zeros to right of point + 1)*4
2011             *
2012             * If the significand is exactly zero, return a properly
2013             * signed zero.
2014             */
2015
2016            String significandString =null;
2017            int signifLength = 0;
2018            int exponentAdjust = 0;
2019            {
2020                int leftDigits  = 0; // number of meaningful digits to
2021                                     // left of "decimal" point
2022                                     // (leading zeros stripped)
2023                int rightDigits = 0; // number of digits to right of
2024                                     // "decimal" point; leading zeros
2025                                     // must always be accounted for
2026                /*
2027                 * The significand is made up of either
2028                 *
2029                 * 1. group 4 entirely (integer portion only)
2030                 *
2031                 * OR
2032                 *
2033                 * 2. the fractional portion from group 7 plus any
2034                 * (optional) integer portions from group 6.
2035                 */
2036                String group4;
2037                if( (group4 = m.group(4)) != null) {  // Integer-only significand
2038                    // Leading zeros never matter on the integer portion
2039                    significandString = stripLeadingZeros(group4);
2040                    leftDigits = significandString.length();
2041                }
2042                else {
2043                    // Group 6 is the optional integer; leading zeros
2044                    // never matter on the integer portion
2045                    String group6 = stripLeadingZeros(m.group(6));
2046                    leftDigits = group6.length();
2047
2048                    // fraction
2049                    String group7 = m.group(7);
2050                    rightDigits = group7.length();
2051
2052                    // Turn "integer.fraction" into "integer"+"fraction"
2053                    significandString =
2054                        ((group6 == null)?"":group6) + // is the null
2055                        // check necessary?
2056                        group7;
2057                }
2058
2059                significandString = stripLeadingZeros(significandString);
2060                signifLength  = significandString.length();
2061
2062                /*
2063                 * Adjust exponent as described above
2064                 */
2065                if (leftDigits >= 1) {  // Cases 1 and 2
2066                    exponentAdjust = 4*(leftDigits - 1);
2067                } else {                // Cases 3 and 4
2068                    exponentAdjust = -4*( rightDigits - signifLength + 1);
2069                }
2070
2071                // If the significand is zero, the exponent doesn't
2072                // matter; return a properly signed zero.
2073
2074                if (signifLength == 0) { // Only zeros in input
2075                    return loadDouble(sign * 0.0);
2076                }
2077            }
2078
2079            //  Extract Exponent
2080            /*
2081             * Use an int to read in the exponent value; this should
2082             * provide more than sufficient range for non-contrived
2083             * inputs.  If reading the exponent in as an int does
2084             * overflow, examine the sign of the exponent and
2085             * significand to determine what to do.
2086             */
2087            String group8 = m.group(8);
2088            boolean positiveExponent = ( group8 == null ) || group8.equals("+");
2089            long unsignedRawExponent;
2090            try {
2091                unsignedRawExponent = Integer.parseInt(m.group(9));
2092            }
2093            catch (NumberFormatException e) {
2094                // At this point, we know the exponent is
2095                // syntactically well-formed as a sequence of
2096                // digits.  Therefore, if an NumberFormatException
2097                // is thrown, it must be due to overflowing int's
2098                // range.  Also, at this point, we have already
2099                // checked for a zero significand.  Thus the signs
2100                // of the exponent and significand determine the
2101                // final result:
2102                //
2103                //                      significand
2104                //                      +               -
2105                // exponent     +       +infinity       -infinity
2106                //              -       +0.0            -0.0
2107                return loadDouble(sign * (positiveExponent ?
2108                                          Double.POSITIVE_INFINITY : 0.0));
2109            }
2110
2111            long rawExponent =
2112                (positiveExponent ? 1L : -1L) * // exponent sign
2113                unsignedRawExponent;            // exponent magnitude
2114
2115            // Calculate partially adjusted exponent
2116            long exponent = rawExponent + exponentAdjust ;
2117
2118            // Starting copying non-zero bits into proper position in
2119            // a long; copy explicit bit too; this will be masked
2120            // later for normal values.
2121
2122            boolean round = false;
2123            boolean sticky = false;
2124            int bitsCopied=0;
2125            int nextShift=0;
2126            long significand=0L;
2127            // First iteration is different, since we only copy
2128            // from the leading significand bit; one more exponent
2129            // adjust will be needed...
2130
2131            // IMPORTANT: make leadingDigit a long to avoid
2132            // surprising shift semantics!
2133            long leadingDigit = getHexDigit(significandString, 0);
2134
2135            /*
2136             * Left shift the leading digit (53 - (bit position of
2137             * leading 1 in digit)); this sets the top bit of the
2138             * significand to 1.  The nextShift value is adjusted
2139             * to take into account the number of bit positions of
2140             * the leadingDigit actually used.  Finally, the
2141             * exponent is adjusted to normalize the significand
2142             * as a binary value, not just a hex value.
2143             */
2144            if (leadingDigit == 1) {
2145                significand |= leadingDigit << 52;
2146                nextShift = 52 - 4;
2147                /* exponent += 0 */     }
2148            else if (leadingDigit <= 3) { // [2, 3]
2149                significand |= leadingDigit << 51;
2150                nextShift = 52 - 5;
2151                exponent += 1;
2152            }
2153            else if (leadingDigit <= 7) { // [4, 7]
2154                significand |= leadingDigit << 50;
2155                nextShift = 52 - 6;
2156                exponent += 2;
2157            }
2158            else if (leadingDigit <= 15) { // [8, f]
2159                significand |= leadingDigit << 49;
2160                nextShift = 52 - 7;
2161                exponent += 3;
2162            } else {
2163                throw new AssertionError("Result from digit conversion too large!");
2164            }
2165            // The preceding if-else could be replaced by a single
2166            // code block based on the high-order bit set in
2167            // leadingDigit.  Given leadingOnePosition,
2168
2169            // significand |= leadingDigit << (SIGNIFICAND_WIDTH - leadingOnePosition);
2170            // nextShift = 52 - (3 + leadingOnePosition);
2171            // exponent += (leadingOnePosition-1);
2172
2173
2174            /*
2175             * Now the exponent variable is equal to the normalized
2176             * binary exponent.  Code below will make representation
2177             * adjustments if the exponent is incremented after
2178             * rounding (includes overflows to infinity) or if the
2179             * result is subnormal.
2180             */
2181
2182            // Copy digit into significand until the significand can't
2183            // hold another full hex digit or there are no more input
2184            // hex digits.
2185            int i = 0;
2186            for(i = 1;
2187                i < signifLength && nextShift >= 0;
2188                i++) {
2189                long currentDigit = getHexDigit(significandString, i);
2190                significand |= (currentDigit << nextShift);
2191                nextShift-=4;
2192            }
2193
2194            // After the above loop, the bulk of the string is copied.
2195            // Now, we must copy any partial hex digits into the
2196            // significand AND compute the round bit and start computing
2197            // sticky bit.
2198
2199            if ( i < signifLength ) { // at least one hex input digit exists
2200                long currentDigit = getHexDigit(significandString, i);
2201
2202                // from nextShift, figure out how many bits need
2203                // to be copied, if any
2204                switch(nextShift) { // must be negative
2205                case -1:
2206                    // three bits need to be copied in; can
2207                    // set round bit
2208                    significand |= ((currentDigit & 0xEL) >> 1);
2209                    round = (currentDigit & 0x1L)  != 0L;
2210                    break;
2211
2212                case -2:
2213                    // two bits need to be copied in; can
2214                    // set round and start sticky
2215                    significand |= ((currentDigit & 0xCL) >> 2);
2216                    round = (currentDigit &0x2L)  != 0L;
2217                    sticky = (currentDigit & 0x1L) != 0;
2218                    break;
2219
2220                case -3:
2221                    // one bit needs to be copied in
2222                    significand |= ((currentDigit & 0x8L)>>3);
2223                    // Now set round and start sticky, if possible
2224                    round = (currentDigit &0x4L)  != 0L;
2225                    sticky = (currentDigit & 0x3L) != 0;
2226                    break;
2227
2228                case -4:
2229                    // all bits copied into significand; set
2230                    // round and start sticky
2231                    round = ((currentDigit & 0x8L) != 0);  // is top bit set?
2232                    // nonzeros in three low order bits?
2233                    sticky = (currentDigit & 0x7L) != 0;
2234                    break;
2235
2236                default:
2237                    throw new AssertionError("Unexpected shift distance remainder.");
2238                    // break;
2239                }
2240
2241                // Round is set; sticky might be set.
2242
2243                // For the sticky bit, it suffices to check the
2244                // current digit and test for any nonzero digits in
2245                // the remaining unprocessed input.
2246                i++;
2247                while(i < signifLength && !sticky) {
2248                    currentDigit =  getHexDigit(significandString,i);
2249                    sticky = sticky || (currentDigit != 0);
2250                    i++;
2251                }
2252
2253            }
2254            // else all of string was seen, round and sticky are
2255            // correct as false.
2256
2257
2258            // Check for overflow and update exponent accordingly.
2259
2260            if (exponent > DoubleConsts.MAX_EXPONENT) {         // Infinite result
2261                // overflow to properly signed infinity
2262                return loadDouble(sign * Double.POSITIVE_INFINITY);
2263            } else {  // Finite return value
2264                if (exponent <= DoubleConsts.MAX_EXPONENT && // (Usually) normal result
2265                    exponent >= DoubleConsts.MIN_EXPONENT) {
2266
2267                    // The result returned in this block cannot be a
2268                    // zero or subnormal; however after the
2269                    // significand is adjusted from rounding, we could
2270                    // still overflow in infinity.
2271
2272                    // AND exponent bits into significand; if the
2273                    // significand is incremented and overflows from
2274                    // rounding, this combination will update the
2275                    // exponent correctly, even in the case of
2276                    // Double.MAX_VALUE overflowing to infinity.
2277
2278                    significand = (( ((long)exponent +
2279                                     (long)DoubleConsts.EXP_BIAS) <<
2280                                     (DoubleConsts.SIGNIFICAND_WIDTH-1))
2281                                   & DoubleConsts.EXP_BIT_MASK) |
2282                        (DoubleConsts.SIGNIF_BIT_MASK & significand);
2283
2284                }  else  {  // Subnormal or zero
2285                    // (exponent < DoubleConsts.MIN_EXPONENT)
2286
2287                    if (exponent < (DoubleConsts.MIN_SUB_EXPONENT -1 )) {
2288                        // No way to round back to nonzero value
2289                        // regardless of significand if the exponent is
2290                        // less than -1075.
2291                        return loadDouble(sign * 0.0);
2292                    } else { //  -1075 <= exponent <= MIN_EXPONENT -1 = -1023
2293                        /*
2294                         * Find bit position to round to; recompute
2295                         * round and sticky bits, and shift
2296                         * significand right appropriately.
2297                         */
2298
2299                        sticky = sticky || round;
2300                        round = false;
2301
2302                        // Number of bits of significand to preserve is
2303                        // exponent - abs_min_exp +1
2304                        // check:
2305                        // -1075 +1074 + 1 = 0
2306                        // -1023 +1074 + 1 = 52
2307
2308                        int bitsDiscarded = 53 -
2309                            ((int)exponent - DoubleConsts.MIN_SUB_EXPONENT + 1);
2310                        assert bitsDiscarded >= 1 && bitsDiscarded <= 53;
2311
2312                        // What to do here:
2313                        // First, isolate the new round bit
2314                        round = (significand & (1L << (bitsDiscarded -1))) != 0L;
2315                        if (bitsDiscarded > 1) {
2316                            // create mask to update sticky bits; low
2317                            // order bitsDiscarded bits should be 1
2318                            long mask = ~((~0L) << (bitsDiscarded -1));
2319                            sticky = sticky || ((significand & mask) != 0L ) ;
2320                        }
2321
2322                        // Now, discard the bits
2323                        significand = significand >> bitsDiscarded;
2324
2325                        significand = (( ((long)(DoubleConsts.MIN_EXPONENT -1) + // subnorm exp.
2326                                          (long)DoubleConsts.EXP_BIAS) <<
2327                                         (DoubleConsts.SIGNIFICAND_WIDTH-1))
2328                                       & DoubleConsts.EXP_BIT_MASK) |
2329                            (DoubleConsts.SIGNIF_BIT_MASK & significand);
2330                    }
2331                }
2332
2333                // The significand variable now contains the currently
2334                // appropriate exponent bits too.
2335
2336                /*
2337                 * Determine if significand should be incremented;
2338                 * making this determination depends on the least
2339                 * significant bit and the round and sticky bits.
2340                 *
2341                 * Round to nearest even rounding table, adapted from
2342                 * table 4.7 in "Computer Arithmetic" by IsraelKoren.
2343                 * The digit to the left of the "decimal" point is the
2344                 * least significant bit, the digits to the right of
2345                 * the point are the round and sticky bits
2346                 *
2347                 * Number       Round(x)
2348                 * x0.00        x0.
2349                 * x0.01        x0.
2350                 * x0.10        x0.
2351                 * x0.11        x1. = x0. +1
2352                 * x1.00        x1.
2353                 * x1.01        x1.
2354                 * x1.10        x1. + 1
2355                 * x1.11        x1. + 1
2356                 */
2357                boolean incremented = false;
2358                boolean leastZero  = ((significand & 1L) == 0L);
2359                if( (  leastZero  && round && sticky ) ||
2360                    ((!leastZero) && round )) {
2361                    incremented = true;
2362                    significand++;
2363                }
2364
2365                loadDouble(FpUtils.rawCopySign(
2366                                               Double.longBitsToDouble(significand),
2367                                               sign));
2368
2369                /*
2370                 * Set roundingDir variable field of fd properly so
2371                 * that the input string can be properly rounded to a
2372                 * float value.  There are two cases to consider:
2373                 *
2374                 * 1. rounding to double discards sticky bit
2375                 * information that would change the result of a float
2376                 * rounding (near halfway case between two floats)
2377                 *
2378                 * 2. rounding to double rounds up when rounding up
2379                 * would not occur when rounding to float.
2380                 *
2381                 * For former case only needs to be considered when
2382                 * the bits rounded away when casting to float are all
2383                 * zero; otherwise, float round bit is properly set
2384                 * and sticky will already be true.
2385                 *
2386                 * The lower exponent bound for the code below is the
2387                 * minimum (normalized) subnormal exponent - 1 since a
2388                 * value with that exponent can round up to the
2389                 * minimum subnormal value and the sticky bit
2390                 * information must be preserved (i.e. case 1).
2391                 */
2392                if ((exponent >= FloatConsts.MIN_SUB_EXPONENT-1) &&
2393                    (exponent <= FloatConsts.MAX_EXPONENT ) ){
2394                    // Outside above exponent range, the float value
2395                    // will be zero or infinity.
2396
2397                    /*
2398                     * If the low-order 28 bits of a rounded double
2399                     * significand are 0, the double could be a
2400                     * half-way case for a rounding to float.  If the
2401                     * double value is a half-way case, the double
2402                     * significand may have to be modified to round
2403                     * the the right float value (see the stickyRound
2404                     * method).  If the rounding to double has lost
2405                     * what would be float sticky bit information, the
2406                     * double significand must be incremented.  If the
2407                     * double value's significand was itself
2408                     * incremented, the float value may end up too
2409                     * large so the increment should be undone.
2410                     */
2411                    if ((significand & 0xfffffffL) ==  0x0L) {
2412                        // For negative values, the sign of the
2413                        // roundDir is the same as for positive values
2414                        // since adding 1 increasing the significand's
2415                        // magnitude and subtracting 1 decreases the
2416                        // significand's magnitude.  If neither round
2417                        // nor sticky is true, the double value is
2418                        // exact and no adjustment is required for a
2419                        // proper float rounding.
2420                        if( round || sticky) {
2421                            if (leastZero) { // prerounding lsb is 0
2422                                // If round and sticky were both true,
2423                                // and the least significant
2424                                // significand bit were 0, the rounded
2425                                // significand would not have its
2426                                // low-order bits be zero.  Therefore,
2427                                // we only need to adjust the
2428                                // significand if round XOR sticky is
2429                                // true.
2430                                if (round ^ sticky) {
2431                                    this.roundDir =  1;
2432                                }
2433                            }
2434                            else { // prerounding lsb is 1
2435                                // If the prerounding lsb is 1 and the
2436                                // resulting significand has its
2437                                // low-order bits zero, the significand
2438                                // was incremented.  Here, we undo the
2439                                // increment, which will ensure the
2440                                // right guard and sticky bits for the
2441                                // float rounding.
2442                                if (round)
2443                                    this.roundDir =  -1;
2444                            }
2445                        }
2446                    }
2447                }
2448
2449                this.fromHex = true;
2450                return this;
2451            }
2452        }
2453    }
2454
2455    /**
2456     * Return <code>s</code> with any leading zeros removed.
2457     */
2458    static String stripLeadingZeros(String s) {
2459        return  s.replaceFirst("^0+", "");
2460    }
2461
2462    /**
2463     * Extract a hexadecimal digit from position <code>position</code>
2464     * of string <code>s</code>.
2465     */
2466    static int getHexDigit(String s, int position) {
2467        int value = Character.digit(s.charAt(position), 16);
2468        if (value <= -1 || value >= 16) {
2469            throw new AssertionError("Unexpected failure of digit conversion of " +
2470                                     s.charAt(position));
2471        }
2472        return value;
2473    }
2474
2475
2476}
2477