1 /*
   2  * Copyright (c) 2003, 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 
  26 package sun.misc;
  27 

  28 import sun.misc.DoubleConsts;
  29 import sun.misc.FloatConsts;
  30 import java.util.regex.*;
  31 
  32 public class FormattedFloatingDecimal{
  33     boolean     isExceptional;
  34     boolean     isNegative;
  35     int         decExponent;  // value set at construction, then immutable
  36     int         decExponentRounded;
  37     char        digits[];
  38     int         nDigits;
  39     int         bigIntExp;
  40     int         bigIntNBits;
  41     boolean     mustSetRoundDir = false;
  42     boolean     fromHex = false;
  43     int         roundDir = 0; // set by doubleValue
  44     int         precision;    // number of digits to the right of decimal
  45 
  46     public enum Form { SCIENTIFIC, COMPATIBLE, DECIMAL_FLOAT, GENERAL };
  47 
  48     private Form form;
  49 
  50     private     FormattedFloatingDecimal( boolean negSign, int decExponent, char []digits, int n, boolean e, int precision, Form form )
  51     {
  52         isNegative = negSign;
  53         isExceptional = e;
  54         this.decExponent = decExponent;
  55         this.digits = digits;
  56         this.nDigits = n;
  57         this.precision = precision;
  58         this.form = form;
  59     }
  60 
  61     /*
  62      * Constants of the implementation
  63      * Most are IEEE-754 related.
  64      * (There are more really boring constants at the end.)
  65      */
  66     static final long   signMask = 0x8000000000000000L;
  67     static final long   expMask  = 0x7ff0000000000000L;
  68     static final long   fractMask= ~(signMask|expMask);
  69     static final int    expShift = 52;
  70     static final int    expBias  = 1023;
  71     static final long   fractHOB = ( 1L<<expShift ); // assumed High-Order bit
  72     static final long   expOne   = ((long)expBias)<<expShift; // exponent of 1.0
  73     static final int    maxSmallBinExp = 62;
  74     static final int    minSmallBinExp = -( 63 / 3 );
  75     static final int    maxDecimalDigits = 15;
  76     static final int    maxDecimalExponent = 308;
  77     static final int    minDecimalExponent = -324;
  78     static final int    bigDecimalExponent = 324; // i.e. abs(minDecimalExponent)
  79 
  80     static final long   highbyte = 0xff00000000000000L;
  81     static final long   highbit  = 0x8000000000000000L;
  82     static final long   lowbytes = ~highbyte;
  83 
  84     static final int    singleSignMask =    0x80000000;
  85     static final int    singleExpMask  =    0x7f800000;
  86     static final int    singleFractMask =   ~(singleSignMask|singleExpMask);
  87     static final int    singleExpShift  =   23;
  88     static final int    singleFractHOB  =   1<<singleExpShift;
  89     static final int    singleExpBias   =   127;
  90     static final int    singleMaxDecimalDigits = 7;
  91     static final int    singleMaxDecimalExponent = 38;
  92     static final int    singleMinDecimalExponent = -45;
  93 
  94     static final int    intDecimalDigits = 9;
  95 
  96 
  97     /*
  98      * count number of bits from high-order 1 bit to low-order 1 bit,
  99      * inclusive.
 100      */
 101     private static int
 102     countBits( long v ){
 103         //
 104         // the strategy is to shift until we get a non-zero sign bit
 105         // then shift until we have no bits left, counting the difference.
 106         // we do byte shifting as a hack. Hope it helps.
 107         //
 108         if ( v == 0L ) return 0;
 109 
 110         while ( ( v & highbyte ) == 0L ){
 111             v <<= 8;
 112         }
 113         while ( v > 0L ) { // i.e. while ((v&highbit) == 0L )
 114             v <<= 1;
 115         }
 116 
 117         int n = 0;
 118         while (( v & lowbytes ) != 0L ){
 119             v <<= 8;
 120             n += 8;
 121         }
 122         while ( v != 0L ){
 123             v <<= 1;
 124             n += 1;
 125         }
 126         return n;
 127     }
 128 
 129     /*
 130      * Keep big powers of 5 handy for future reference.
 131      */
 132     private static FDBigInt b5p[];
 133 
 134     private static synchronized FDBigInt
 135     big5pow( int p ){
 136         assert p >= 0 : p; // negative power of 5
 137         if ( b5p == null ){
 138             b5p = new FDBigInt[ p+1 ];
 139         }else if (b5p.length <= p ){
 140             FDBigInt t[] = new FDBigInt[ p+1 ];
 141             System.arraycopy( b5p, 0, t, 0, b5p.length );
 142             b5p = t;
 143         }
 144         if ( b5p[p] != null )
 145             return b5p[p];
 146         else if ( p < small5pow.length )
 147             return b5p[p] = new FDBigInt( small5pow[p] );
 148         else if ( p < long5pow.length )
 149             return b5p[p] = new FDBigInt( long5pow[p] );
 150         else {
 151             // construct the value.
 152             // recursively.
 153             int q, r;
 154             // in order to compute 5^p,
 155             // compute its square root, 5^(p/2) and square.
 156             // or, let q = p / 2, r = p -q, then
 157             // 5^p = 5^(q+r) = 5^q * 5^r
 158             q = p >> 1;
 159             r = p - q;
 160             FDBigInt bigq =  b5p[q];
 161             if ( bigq == null )
 162                 bigq = big5pow ( q );
 163             if ( r < small5pow.length ){
 164                 return (b5p[p] = bigq.mult( small5pow[r] ) );
 165             }else{
 166                 FDBigInt bigr = b5p[ r ];
 167                 if ( bigr == null )
 168                     bigr = big5pow( r );
 169                 return (b5p[p] = bigq.mult( bigr ) );
 170             }
 171         }
 172     }
 173 
 174     //
 175     // a common operation
 176     //
 177     private static FDBigInt
 178     multPow52( FDBigInt v, int p5, int p2 ){
 179         if ( p5 != 0 ){
 180             if ( p5 < small5pow.length ){
 181                 v = v.mult( small5pow[p5] );
 182             } else {
 183                 v = v.mult( big5pow( p5 ) );
 184             }
 185         }
 186         if ( p2 != 0 ){
 187             v.lshiftMe( p2 );
 188         }
 189         return v;
 190     }
 191 
 192     //
 193     // another common operation
 194     //
 195     private static FDBigInt
 196     constructPow52( int p5, int p2 ){
 197         FDBigInt v = new FDBigInt( big5pow( p5 ) );
 198         if ( p2 != 0 ){
 199             v.lshiftMe( p2 );
 200         }
 201         return v;
 202     }
 203 
 204     /*
 205      * Make a floating double into a FDBigInt.
 206      * This could also be structured as a FDBigInt
 207      * constructor, but we'd have to build a lot of knowledge
 208      * about floating-point representation into it, and we don't want to.
 209      *
 210      * AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES
 211      * bigIntExp and bigIntNBits
 212      *
 213      */
 214     private FDBigInt
 215     doubleToBigInt( double dval ){
 216         long lbits = Double.doubleToLongBits( dval ) & ~signMask;
 217         int binexp = (int)(lbits >>> expShift);
 218         lbits &= fractMask;
 219         if ( binexp > 0 ){
 220             lbits |= fractHOB;
 221         } else {
 222             assert lbits != 0L : lbits; // doubleToBigInt(0.0)
 223             binexp +=1;
 224             while ( (lbits & fractHOB ) == 0L){
 225                 lbits <<= 1;
 226                 binexp -= 1;
 227             }
 228         }
 229         binexp -= expBias;
 230         int nbits = countBits( lbits );
 231         /*
 232          * We now know where the high-order 1 bit is,
 233          * and we know how many there are.
 234          */
 235         int lowOrderZeros = expShift+1-nbits;
 236         lbits >>>= lowOrderZeros;
 237 
 238         bigIntExp = binexp+1-nbits;
 239         bigIntNBits = nbits;
 240         return new FDBigInt( lbits );
 241     }
 242 
 243     /*
 244      * Compute a number that is the ULP of the given value,
 245      * for purposes of addition/subtraction. Generally easy.
 246      * More difficult if subtracting and the argument
 247      * is a normalized a power of 2, as the ULP changes at these points.
 248      */
 249     private static double ulp( double dval, boolean subtracting ){
 250         long lbits = Double.doubleToLongBits( dval ) & ~signMask;
 251         int binexp = (int)(lbits >>> expShift);
 252         double ulpval;
 253         if ( subtracting && ( binexp >= expShift ) && ((lbits&fractMask) == 0L) ){
 254             // for subtraction from normalized, powers of 2,
 255             // use next-smaller exponent
 256             binexp -= 1;
 257         }
 258         if ( binexp > expShift ){
 259             ulpval = Double.longBitsToDouble( ((long)(binexp-expShift))<<expShift );
 260         } else if ( binexp == 0 ){
 261             ulpval = Double.MIN_VALUE;
 262         } else {
 263             ulpval = Double.longBitsToDouble( 1L<<(binexp-1) );
 264         }
 265         if ( subtracting ) ulpval = - ulpval;
 266 
 267         return ulpval;
 268     }
 269 
 270     /*
 271      * Round a double to a float.
 272      * In addition to the fraction bits of the double,
 273      * look at the class instance variable roundDir,
 274      * which should help us avoid double-rounding error.
 275      * roundDir was set in hardValueOf if the estimate was
 276      * close enough, but not exact. It tells us which direction
 277      * of rounding is preferred.
 278      */
 279     float
 280     stickyRound( double dval ){
 281         long lbits = Double.doubleToLongBits( dval );
 282         long binexp = lbits & expMask;
 283         if ( binexp == 0L || binexp == expMask ){
 284             // what we have here is special.
 285             // don't worry, the right thing will happen.
 286             return (float) dval;
 287         }
 288         lbits += (long)roundDir; // hack-o-matic.
 289         return (float)Double.longBitsToDouble( lbits );
 290     }
 291 
 292 
 293     /*
 294      * This is the easy subcase --
 295      * all the significant bits, after scaling, are held in lvalue.
 296      * negSign and decExponent tell us what processing and scaling
 297      * has already been done. Exceptional cases have already been
 298      * stripped out.
 299      * In particular:
 300      * lvalue is a finite number (not Inf, nor NaN)
 301      * lvalue > 0L (not zero, nor negative).
 302      *
 303      * The only reason that we develop the digits here, rather than
 304      * calling on Long.toString() is that we can do it a little faster,
 305      * and besides want to treat trailing 0s specially. If Long.toString
 306      * changes, we should re-evaluate this strategy!
 307      */
 308     private void
 309     developLongDigits( int decExponent, long lvalue, long insignificant ){
 310         char digits[];
 311         int  ndigits;
 312         int  digitno;
 313         int  c;
 314         //
 315         // Discard non-significant low-order bits, while rounding,
 316         // up to insignificant value.
 317         int i;
 318         for ( i = 0; insignificant >= 10L; i++ )
 319             insignificant /= 10L;
 320         if ( i != 0 ){
 321             long pow10 = long5pow[i] << i; // 10^i == 5^i * 2^i;
 322             long residue = lvalue % pow10;
 323             lvalue /= pow10;
 324             decExponent += i;
 325             if ( residue >= (pow10>>1) ){
 326                 // round up based on the low-order bits we're discarding
 327                 lvalue++;
 328             }
 329         }
 330         if ( lvalue <= Integer.MAX_VALUE ){
 331             assert lvalue > 0L : lvalue; // lvalue <= 0
 332             // even easier subcase!
 333             // can do int arithmetic rather than long!
 334             int  ivalue = (int)lvalue;
 335             ndigits = 10;
 336             digits = (char[])(perThreadBuffer.get());
 337             digitno = ndigits-1;
 338             c = ivalue%10;
 339             ivalue /= 10;
 340             while ( c == 0 ){
 341                 decExponent++;
 342                 c = ivalue%10;
 343                 ivalue /= 10;
 344             }
 345             while ( ivalue != 0){
 346                 digits[digitno--] = (char)(c+'0');
 347                 decExponent++;
 348                 c = ivalue%10;
 349                 ivalue /= 10;
 350             }
 351             digits[digitno] = (char)(c+'0');
 352         } else {
 353             // same algorithm as above (same bugs, too )
 354             // but using long arithmetic.
 355             ndigits = 20;
 356             digits = (char[])(perThreadBuffer.get());
 357             digitno = ndigits-1;
 358             c = (int)(lvalue%10L);
 359             lvalue /= 10L;
 360             while ( c == 0 ){
 361                 decExponent++;
 362                 c = (int)(lvalue%10L);
 363                 lvalue /= 10L;
 364             }
 365             while ( lvalue != 0L ){
 366                 digits[digitno--] = (char)(c+'0');
 367                 decExponent++;
 368                 c = (int)(lvalue%10L);
 369                 lvalue /= 10;
 370             }
 371             digits[digitno] = (char)(c+'0');
 372         }
 373         char result [];
 374         ndigits -= digitno;
 375         result = new char[ ndigits ];
 376         System.arraycopy( digits, digitno, result, 0, ndigits );
 377         this.digits = result;
 378         this.decExponent = decExponent+1;
 379         this.nDigits = ndigits;
 380     }
 381 
 382     //
 383     // add one to the least significant digit.
 384     // in the unlikely event there is a carry out,
 385     // deal with it.
 386     // assert that this will only happen where there
 387     // is only one digit, e.g. (float)1e-44 seems to do it.
 388     //
 389     private void
 390     roundup(){
 391         int i;
 392         int q = digits[ i = (nDigits-1)];
 393         if ( q == '9' ){
 394             while ( q == '9' && i > 0 ){
 395                 digits[i] = '0';
 396                 q = digits[--i];
 397             }
 398             if ( q == '9' ){
 399                 // carryout! High-order 1, rest 0s, larger exp.
 400                 decExponent += 1;
 401                 digits[0] = '1';
 402                 return;
 403             }
 404             // else fall through.
 405         }
 406         digits[i] = (char)(q+1);
 407     }
 408 
 409     // Given the desired number of digits predict the result's exponent.
 410     private int checkExponent(int length) {
 411         if (length >= nDigits || length < 0)
 412             return decExponent;
 413 
 414         for (int i = 0; i < length; i++)
 415             if (digits[i] != '9')
 416                 // a '9' anywhere in digits will absorb the round
 417                 return decExponent;
 418         return decExponent + (digits[length] >= '5' ? 1 : 0);
 419     }
 420 
 421     // Unlike roundup(), this method does not modify digits.  It also
 422     // rounds at a particular precision.
 423     private char [] applyPrecision(int length) {
 424         char [] result = new char[nDigits];
 425         for (int i = 0; i < result.length; i++) result[i] = '0';
 426 
 427         if (length >= nDigits || length < 0) {
 428             // no rounding necessary
 429             System.arraycopy(digits, 0, result, 0, nDigits);
 430             return result;
 431         }
 432         if (length == 0) {
 433             // only one digit (0 or 1) is returned because the precision
 434             // excludes all significant digits
 435             if (digits[0] >= '5') {
 436                 result[0] = '1';
 437             }
 438             return result;
 439         }
 440 
 441         int i = length;
 442         int q = digits[i];
 443         if (q >= '5' && i > 0) {
 444             q = digits[--i];
 445             if ( q == '9' ) {
 446                 while ( q == '9' && i > 0 ){
 447                     q = digits[--i];
 448                 }
 449                 if ( q == '9' ){
 450                     // carryout! High-order 1, rest 0s, larger exp.
 451                     result[0] = '1';
 452                     return result;
 453                 }
 454             }
 455             result[i] = (char)(q + 1);
 456         }
 457         while (--i >= 0) {
 458             result[i] = digits[i];
 459         }
 460         return result;
 461     }
 462 
 463     /*
 464      * FIRST IMPORTANT CONSTRUCTOR: DOUBLE
 465      */
 466     public FormattedFloatingDecimal( double d )
 467     {
 468         this(d, Integer.MAX_VALUE, Form.COMPATIBLE);
 469     }
 470 
 471     public FormattedFloatingDecimal( double d, int precision, Form form )
 472     {
 473         long    dBits = Double.doubleToLongBits( d );
 474         long    fractBits;
 475         int     binExp;
 476         int     nSignificantBits;
 477 
 478         this.precision = precision;
 479         this.form      = form;
 480 
 481         // discover and delete sign
 482         if ( (dBits&signMask) != 0 ){
 483             isNegative = true;
 484             dBits ^= signMask;
 485         } else {
 486             isNegative = false;
 487         }
 488         // Begin to unpack
 489         // Discover obvious special cases of NaN and Infinity.
 490         binExp = (int)( (dBits&expMask) >> expShift );
 491         fractBits = dBits&fractMask;
 492         if ( binExp == (int)(expMask>>expShift) ) {
 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;
 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 == 0L ){
 510                 // not a denorm, just a 0!
 511                 decExponent = 0;
 512                 digits = zero;
 513                 nDigits = 1;
 514                 return;
 515             }
 516             while ( (fractBits&fractHOB) == 0L ){
 517                 fractBits <<= 1;
 518                 binExp -= 1;
 519             }
 520             nSignificantBits = expShift + binExp +1; // recall binExp is  - shift count.
 521             binExp += 1;
 522         } else {
 523             fractBits |= fractHOB;
 524             nSignificantBits = expShift+1;
 525         }
 526         binExp -= expBias;
 527         // call the routine that actually does all the hard work.
 528         dtoa( binExp, fractBits, nSignificantBits );
 529     }
 530 
 531     /*
 532      * SECOND IMPORTANT CONSTRUCTOR: SINGLE
 533      */
 534     public FormattedFloatingDecimal( float f )
 535     {
 536         this(f, Integer.MAX_VALUE, Form.COMPATIBLE);
 537     }
 538     public FormattedFloatingDecimal( float f, int precision, Form form )
 539     {
 540         int     fBits = Float.floatToIntBits( f );
 541         int     fractBits;
 542         int     binExp;
 543         int     nSignificantBits;
 544 
 545         this.precision = precision;
 546         this.form      = form;
 547 
 548         // discover and delete sign
 549         if ( (fBits&singleSignMask) != 0 ){
 550             isNegative = true;
 551             fBits ^= singleSignMask;
 552         } else {
 553             isNegative = false;
 554         }
 555         // Begin to unpack
 556         // Discover obvious special cases of NaN and Infinity.
 557         binExp = (int)( (fBits&singleExpMask) >> singleExpShift );
 558         fractBits = fBits&singleFractMask;
 559         if ( binExp == (int)(singleExpMask>>singleExpShift) ) {
 560             isExceptional = true;
 561             if ( fractBits == 0L ){
 562                 digits =  infinity;
 563             } else {
 564                 digits = notANumber;
 565                 isNegative = false; // NaN has no sign!
 566             }
 567             nDigits = digits.length;
 568             return;
 569         }
 570         isExceptional = false;
 571         // Finish unpacking
 572         // Normalize denormalized numbers.
 573         // Insert assumed high-order bit for normalized numbers.
 574         // Subtract exponent bias.
 575         if ( binExp == 0 ){
 576             if ( fractBits == 0 ){
 577                 // not a denorm, just a 0!
 578                 decExponent = 0;
 579                 digits = zero;
 580                 nDigits = 1;
 581                 return;
 582             }
 583             while ( (fractBits&singleFractHOB) == 0 ){
 584                 fractBits <<= 1;
 585                 binExp -= 1;
 586             }
 587             nSignificantBits = singleExpShift + binExp +1; // recall binExp is  - shift count.
 588             binExp += 1;
 589         } else {
 590             fractBits |= singleFractHOB;
 591             nSignificantBits = singleExpShift+1;
 592         }
 593         binExp -= singleExpBias;
 594         // call the routine that actually does all the hard work.
 595         dtoa( binExp, ((long)fractBits)<<(expShift-singleExpShift), nSignificantBits );
 596     }
 597 
 598     private void
 599     dtoa( int binExp, long fractBits, int nSignificantBits )
 600     {
 601         int     nFractBits; // number of significant bits of fractBits;
 602         int     nTinyBits;  // number of these to the right of the point.
 603         int     decExp;
 604 
 605         // Examine number. Determine if it is an easy case,
 606         // which we can do pretty trivially using float/long conversion,
 607         // or whether we must do real work.
 608         nFractBits = countBits( fractBits );
 609         nTinyBits = Math.max( 0, nFractBits - binExp - 1 );
 610         if ( binExp <= maxSmallBinExp && binExp >= minSmallBinExp ){
 611             // Look more closely at the number to decide if,
 612             // with scaling by 10^nTinyBits, the result will fit in
 613             // a long.
 614             if ( (nTinyBits < long5pow.length) && ((nFractBits + n5bits[nTinyBits]) < 64 ) ){
 615                 /*
 616                  * We can do this:
 617                  * take the fraction bits, which are normalized.
 618                  * (a) nTinyBits == 0: Shift left or right appropriately
 619                  *     to align the binary point at the extreme right, i.e.
 620                  *     where a long int point is expected to be. The integer
 621                  *     result is easily converted to a string.
 622                  * (b) nTinyBits > 0: Shift right by expShift-nFractBits,
 623                  *     which effectively converts to long and scales by
 624                  *     2^nTinyBits. Then multiply by 5^nTinyBits to
 625                  *     complete the scaling. We know this won't overflow
 626                  *     because we just counted the number of bits necessary
 627                  *     in the result. The integer you get from this can
 628                  *     then be converted to a string pretty easily.
 629                  */
 630                 long halfULP;
 631                 if ( nTinyBits == 0 ) {
 632                     if ( binExp > nSignificantBits ){
 633                         halfULP = 1L << ( binExp-nSignificantBits-1);
 634                     } else {
 635                         halfULP = 0L;
 636                     }
 637                     if ( binExp >= expShift ){
 638                         fractBits <<= (binExp-expShift);
 639                     } else {
 640                         fractBits >>>= (expShift-binExp) ;
 641                     }
 642                     developLongDigits( 0, fractBits, halfULP );
 643                     return;
 644                 }
 645                 /*
 646                  * The following causes excess digits to be printed
 647                  * out in the single-float case. Our manipulation of
 648                  * halfULP here is apparently not correct. If we
 649                  * better understand how this works, perhaps we can
 650                  * use this special case again. But for the time being,
 651                  * we do not.
 652                  * else {
 653                  *     fractBits >>>= expShift+1-nFractBits;
 654                  *     fractBits *= long5pow[ nTinyBits ];
 655                  *     halfULP = long5pow[ nTinyBits ] >> (1+nSignificantBits-nFractBits);
 656                  *     developLongDigits( -nTinyBits, fractBits, halfULP );
 657                  *     return;
 658                  * }
 659                  */
 660             }
 661         }
 662         /*
 663          * This is the hard case. We are going to compute large positive
 664          * integers B and S and integer decExp, s.t.
 665          *      d = ( B / S ) * 10^decExp
 666          *      1 <= B / S < 10
 667          * Obvious choices are:
 668          *      decExp = floor( log10(d) )
 669          *      B      = d * 2^nTinyBits * 10^max( 0, -decExp )
 670          *      S      = 10^max( 0, decExp) * 2^nTinyBits
 671          * (noting that nTinyBits has already been forced to non-negative)
 672          * I am also going to compute a large positive integer
 673          *      M      = (1/2^nSignificantBits) * 2^nTinyBits * 10^max( 0, -decExp )
 674          * i.e. M is (1/2) of the ULP of d, scaled like B.
 675          * When we iterate through dividing B/S and picking off the
 676          * quotient bits, we will know when to stop when the remainder
 677          * is <= M.
 678          *
 679          * We keep track of powers of 2 and powers of 5.
 680          */
 681 
 682         /*
 683          * Estimate decimal exponent. (If it is small-ish,
 684          * we could double-check.)
 685          *
 686          * First, scale the mantissa bits such that 1 <= d2 < 2.
 687          * We are then going to estimate
 688          *          log10(d2) ~=~  (d2-1.5)/1.5 + log(1.5)
 689          * and so we can estimate
 690          *      log10(d) ~=~ log10(d2) + binExp * log10(2)
 691          * take the floor and call it decExp.
 692          * FIXME -- use more precise constants here. It costs no more.
 693          */
 694         double d2 = Double.longBitsToDouble(
 695             expOne | ( fractBits &~ fractHOB ) );
 696         decExp = (int)Math.floor(
 697             (d2-1.5D)*0.289529654D + 0.176091259 + (double)binExp * 0.301029995663981 );
 698         int B2, B5; // powers of 2 and powers of 5, respectively, in B
 699         int S2, S5; // powers of 2 and powers of 5, respectively, in S
 700         int M2, M5; // powers of 2 and powers of 5, respectively, in M
 701         int Bbits; // binary digits needed to represent B, approx.
 702         int tenSbits; // binary digits needed to represent 10*S, approx.
 703         FDBigInt Sval, Bval, Mval;
 704 
 705         B5 = Math.max( 0, -decExp );
 706         B2 = B5 + nTinyBits + binExp;
 707 
 708         S5 = Math.max( 0, decExp );
 709         S2 = S5 + nTinyBits;
 710 
 711         M5 = B5;
 712         M2 = B2 - nSignificantBits;
 713 
 714         /*
 715          * the long integer fractBits contains the (nFractBits) interesting
 716          * bits from the mantissa of d ( hidden 1 added if necessary) followed
 717          * by (expShift+1-nFractBits) zeros. In the interest of compactness,
 718          * I will shift out those zeros before turning fractBits into a
 719          * FDBigInt. The resulting whole number will be
 720          *      d * 2^(nFractBits-1-binExp).
 721          */
 722         fractBits >>>= (expShift+1-nFractBits);
 723         B2 -= nFractBits-1;
 724         int common2factor = Math.min( B2, S2 );
 725         B2 -= common2factor;
 726         S2 -= common2factor;
 727         M2 -= common2factor;
 728 
 729         /*
 730          * HACK!! For exact powers of two, the next smallest number
 731          * is only half as far away as we think (because the meaning of
 732          * ULP changes at power-of-two bounds) for this reason, we
 733          * hack M2. Hope this works.
 734          */
 735         if ( nFractBits == 1 )
 736             M2 -= 1;
 737 
 738         if ( M2 < 0 ){
 739             // oops.
 740             // since we cannot scale M down far enough,
 741             // we must scale the other values up.
 742             B2 -= M2;
 743             S2 -= M2;
 744             M2 =  0;
 745         }
 746         /*
 747          * Construct, Scale, iterate.
 748          * Some day, we'll write a stopping test that takes
 749          * account of the assymetry of the spacing of floating-point
 750          * numbers below perfect powers of 2
 751          * 26 Sept 96 is not that day.
 752          * So we use a symmetric test.
 753          */
 754         char digits[] = this.digits = new char[18];
 755         int  ndigit = 0;
 756         boolean low, high;
 757         long lowDigitDifference;
 758         int  q;
 759 
 760         /*
 761          * Detect the special cases where all the numbers we are about
 762          * to compute will fit in int or long integers.
 763          * In these cases, we will avoid doing FDBigInt arithmetic.
 764          * We use the same algorithms, except that we "normalize"
 765          * our FDBigInts before iterating. This is to make division easier,
 766          * as it makes our fist guess (quotient of high-order words)
 767          * more accurate!
 768          *
 769          * Some day, we'll write a stopping test that takes
 770          * account of the assymetry of the spacing of floating-point
 771          * numbers below perfect powers of 2
 772          * 26 Sept 96 is not that day.
 773          * So we use a symmetric test.
 774          */
 775         Bbits = nFractBits + B2 + (( B5 < n5bits.length )? n5bits[B5] : ( B5*3 ));
 776         tenSbits = S2+1 + (( (S5+1) < n5bits.length )? n5bits[(S5+1)] : ( (S5+1)*3 ));
 777         if ( Bbits < 64 && tenSbits < 64){
 778             if ( Bbits < 32 && tenSbits < 32){
 779                 // wa-hoo! They're all ints!
 780                 int b = ((int)fractBits * small5pow[B5] ) << B2;
 781                 int s = small5pow[S5] << S2;
 782                 int m = small5pow[M5] << M2;
 783                 int tens = s * 10;
 784                 /*
 785                  * Unroll the first iteration. If our decExp estimate
 786                  * was too high, our first quotient will be zero. In this
 787                  * case, we discard it and decrement decExp.
 788                  */
 789                 ndigit = 0;
 790                 q = b / s;
 791                 b = 10 * ( b % s );
 792                 m *= 10;
 793                 low  = (b <  m );
 794                 high = (b+m > tens );
 795                 assert q < 10 : q; // excessively large digit
 796                 if ( (q == 0) && ! high ){
 797                     // oops. Usually ignore leading zero.
 798                     decExp--;
 799                 } else {
 800                     digits[ndigit++] = (char)('0' + q);
 801                 }
 802                 /*
 803                  * HACK! Java spec sez that we always have at least
 804                  * one digit after the . in either F- or E-form output.
 805                  * Thus we will need more than one digit if we're using
 806                  * E-form
 807                  */
 808                 if (! (form == Form.COMPATIBLE && -3 < decExp && decExp < 8)) {
 809                     high = low = false;
 810                 }
 811                 while( ! low && ! high ){
 812                     q = b / s;
 813                     b = 10 * ( b % s );
 814                     m *= 10;
 815                     assert q < 10 : q; // excessively large digit
 816                     if ( m > 0L ){
 817                         low  = (b <  m );
 818                         high = (b+m > tens );
 819                     } else {
 820                         // hack -- m might overflow!
 821                         // in this case, it is certainly > b,
 822                         // which won't
 823                         // and b+m > tens, too, since that has overflowed
 824                         // either!
 825                         low = true;
 826                         high = true;
 827                     }
 828                     digits[ndigit++] = (char)('0' + q);
 829                 }
 830                 lowDigitDifference = (b<<1) - tens;
 831             } else {
 832                 // still good! they're all longs!
 833                 long b = (fractBits * long5pow[B5] ) << B2;
 834                 long s = long5pow[S5] << S2;
 835                 long m = long5pow[M5] << M2;
 836                 long tens = s * 10L;
 837                 /*
 838                  * Unroll the first iteration. If our decExp estimate
 839                  * was too high, our first quotient will be zero. In this
 840                  * case, we discard it and decrement decExp.
 841                  */
 842                 ndigit = 0;
 843                 q = (int) ( b / s );
 844                 b = 10L * ( b % s );
 845                 m *= 10L;
 846                 low  = (b <  m );
 847                 high = (b+m > tens );
 848                 assert q < 10 : q; // excessively large digit
 849                 if ( (q == 0) && ! high ){
 850                     // oops. Usually ignore leading zero.
 851                     decExp--;
 852                 } else {
 853                     digits[ndigit++] = (char)('0' + q);
 854                 }
 855                 /*
 856                  * HACK! Java spec sez that we always have at least
 857                  * one digit after the . in either F- or E-form output.
 858                  * Thus we will need more than one digit if we're using
 859                  * E-form
 860                  */
 861                 if (! (form == Form.COMPATIBLE && -3 < decExp && decExp < 8)) {
 862                     high = low = false;
 863                 }
 864                 while( ! low && ! high ){
 865                     q = (int) ( b / s );
 866                     b = 10 * ( b % s );
 867                     m *= 10;
 868                     assert q < 10 : q;  // excessively large digit
 869                     if ( m > 0L ){
 870                         low  = (b <  m );
 871                         high = (b+m > tens );
 872                     } else {
 873                         // hack -- m might overflow!
 874                         // in this case, it is certainly > b,
 875                         // which won't
 876                         // and b+m > tens, too, since that has overflowed
 877                         // either!
 878                         low = true;
 879                         high = true;
 880                     }
 881                     digits[ndigit++] = (char)('0' + q);
 882                 }
 883                 lowDigitDifference = (b<<1) - tens;
 884             }
 885         } else {
 886             FDBigInt tenSval;
 887             int  shiftBias;
 888 
 889             /*
 890              * We really must do FDBigInt arithmetic.
 891              * Fist, construct our FDBigInt initial values.
 892              */
 893             Bval = multPow52( new FDBigInt( fractBits  ), B5, B2 );
 894             Sval = constructPow52( S5, S2 );
 895             Mval = constructPow52( M5, M2 );
 896 
 897 
 898             // normalize so that division works better
 899             Bval.lshiftMe( shiftBias = Sval.normalizeMe() );
 900             Mval.lshiftMe( shiftBias );
 901             tenSval = Sval.mult( 10 );
 902             /*
 903              * Unroll the first iteration. If our decExp estimate
 904              * was too high, our first quotient will be zero. In this
 905              * case, we discard it and decrement decExp.
 906              */
 907             ndigit = 0;
 908             q = Bval.quoRemIteration( Sval );
 909             Mval = Mval.mult( 10 );
 910             low  = (Bval.cmp( Mval ) < 0);
 911             high = (Bval.add( Mval ).cmp( tenSval ) > 0 );
 912             assert q < 10 : q; // excessively large digit
 913             if ( (q == 0) && ! high ){
 914                 // oops. Usually ignore leading zero.
 915                 decExp--;
 916             } else {
 917                 digits[ndigit++] = (char)('0' + q);
 918             }
 919             /*
 920              * HACK! Java spec sez that we always have at least
 921              * one digit after the . in either F- or E-form output.
 922              * Thus we will need more than one digit if we're using
 923              * E-form
 924              */
 925             if (! (form == Form.COMPATIBLE && -3 < decExp && decExp < 8)) {
 926                 high = low = false;
 927             }
 928             while( ! low && ! high ){
 929                 q = Bval.quoRemIteration( Sval );
 930                 Mval = Mval.mult( 10 );
 931                 assert q < 10 : q;  // excessively large digit
 932                 low  = (Bval.cmp( Mval ) < 0);
 933                 high = (Bval.add( Mval ).cmp( tenSval ) > 0 );
 934                 digits[ndigit++] = (char)('0' + q);
 935             }
 936             if ( high && low ){
 937                 Bval.lshiftMe(1);
 938                 lowDigitDifference = Bval.cmp(tenSval);
 939             } else
 940                 lowDigitDifference = 0L; // this here only for flow analysis!
 941         }
 942         this.decExponent = decExp+1;
 943         this.digits = digits;
 944         this.nDigits = ndigit;
 945         /*
 946          * Last digit gets rounded based on stopping condition.
 947          */
 948         if ( high ){
 949             if ( low ){
 950                 if ( lowDigitDifference == 0L ){
 951                     // it's a tie!
 952                     // choose based on which digits we like.
 953                     if ( (digits[nDigits-1]&1) != 0 ) roundup();
 954                 } else if ( lowDigitDifference > 0 ){
 955                     roundup();
 956                 }
 957             } else {
 958                 roundup();
 959             }
 960         }
 961     }
 962 
 963     public String
 964     toString(){
 965         // most brain-dead version
 966         StringBuffer result = new StringBuffer( nDigits+8 );
 967         if ( isNegative ){ result.append( '-' ); }
 968         if ( isExceptional ){
 969             result.append( digits, 0, nDigits );
 970         } else {
 971             result.append( "0.");
 972             result.append( digits, 0, nDigits );
 973             result.append('e');
 974             result.append( decExponent );
 975         }
 976         return new String(result);
 977     }
 978 
 979     // returns the exponent before rounding
 980     public int getExponent() {
 981         return decExponent - 1;
 982     }
 983 
 984     // returns the exponent after rounding has been done by applyPrecision
 985     public int getExponentRounded() {
 986         return decExponentRounded - 1;
 987     }
 988 
 989     public int getChars(char[] result) {
 990         assert nDigits <= 19 : nDigits; // generous bound on size of nDigits
 991         int i = 0;
 992         if (isNegative) { result[0] = '-'; i = 1; }
 993         if (isExceptional) {
 994             System.arraycopy(digits, 0, result, i, nDigits);
 995             i += nDigits;
 996         } else {
 997             char digits [] = this.digits;
 998             int exp = decExponent;
 999             switch (form) {
1000             case COMPATIBLE:
1001                 break;
1002             case DECIMAL_FLOAT:
1003                 exp = checkExponent(decExponent + precision);
1004                 digits = applyPrecision(decExponent + precision);
1005                 break;
1006             case SCIENTIFIC:
1007                 exp = checkExponent(precision + 1);
1008                 digits = applyPrecision(precision + 1);
1009                 break;
1010             case GENERAL:
1011                 exp = checkExponent(precision);
1012                 digits = applyPrecision(precision);
1013                 // adjust precision to be the number of digits to right of decimal
1014                 // the real exponent to be output is actually exp - 1, not exp
1015                 if (exp - 1 < -4 || exp - 1 >= precision) {
1016                     form = Form.SCIENTIFIC;
1017                     precision--;
1018                 } else {
1019                     form = Form.DECIMAL_FLOAT;
1020                     precision = precision - exp;
1021                 }
1022                 break;
1023             default:
1024                 assert false;
1025             }
1026             decExponentRounded = exp;
1027 
1028             if (exp > 0
1029                 && ((form == Form.COMPATIBLE && (exp < 8))
1030                     || (form == Form.DECIMAL_FLOAT)))
1031             {
1032                 // print digits.digits.
1033                 int charLength = Math.min(nDigits, exp);
1034                 System.arraycopy(digits, 0, result, i, charLength);
1035                 i += charLength;
1036                 if (charLength < exp) {
1037                     charLength = exp-charLength;
1038                     for (int nz = 0; nz < charLength; nz++)
1039                         result[i++] = '0';
1040                     // Do not append ".0" for formatted floats since the user
1041                     // may request that it be omitted. It is added as necessary
1042                     // by the Formatter.
1043                     if (form == Form.COMPATIBLE) {
1044                         result[i++] = '.';
1045                         result[i++] = '0';
1046                     }
1047                 } else {
1048                     // Do not append ".0" for formatted floats since the user
1049                     // may request that it be omitted. It is added as necessary
1050                     // by the Formatter.
1051                     if (form == Form.COMPATIBLE) {
1052                         result[i++] = '.';
1053                         if (charLength < nDigits) {
1054                             int t = Math.min(nDigits - charLength, precision);
1055                             System.arraycopy(digits, charLength, result, i, t);
1056                             i += t;
1057                         } else {
1058                             result[i++] = '0';
1059                         }
1060                     } else {
1061                         int t = Math.min(nDigits - charLength, precision);
1062                         if (t > 0) {
1063                             result[i++] = '.';
1064                             System.arraycopy(digits, charLength, result, i, t);
1065                             i += t;
1066                         }
1067                     }
1068                 }
1069             } else if (exp <= 0
1070                        && ((form == Form.COMPATIBLE && exp > -3)
1071                        || (form == Form.DECIMAL_FLOAT)))
1072             {
1073                 // print 0.0* digits
1074                 result[i++] = '0';
1075                 if (exp != 0) {
1076                     // write '0' s before the significant digits
1077                     int t = Math.min(-exp, precision);
1078                     if (t > 0) {
1079                         result[i++] = '.';
1080                         for (int nz = 0; nz < t; nz++)
1081                             result[i++] = '0';
1082                     }
1083                 }
1084                 int t = Math.min(digits.length, precision + exp);
1085                 if (t > 0) {
1086                     if (i == 1)
1087                         result[i++] = '.';
1088                     // copy only when significant digits are within the precision
1089                     System.arraycopy(digits, 0, result, i, t);
1090                     i += t;
1091                 }
1092             } else {
1093                 result[i++] = digits[0];
1094                 if (form == Form.COMPATIBLE) {
1095                     result[i++] = '.';
1096                     if (nDigits > 1) {
1097                         System.arraycopy(digits, 1, result, i, nDigits-1);
1098                         i += nDigits-1;
1099                     } else {
1100                         result[i++] = '0';
1101                     }
1102                     result[i++] = 'E';
1103                 } else {
1104                     if (nDigits > 1) {
1105                         int t = Math.min(nDigits -1, precision);
1106                         if (t > 0) {
1107                             result[i++] = '.';
1108                             System.arraycopy(digits, 1, result, i, t);
1109                             i += t;
1110                         }
1111                     }
1112                     result[i++] = 'e';
1113                 }
1114                 int e;
1115                 if (exp <= 0) {
1116                     result[i++] = '-';
1117                     e = -exp+1;
1118                 } else {
1119                     if (form != Form.COMPATIBLE)
1120                         result[i++] = '+';
1121                     e = exp-1;
1122                 }
1123                 // decExponent has 1, 2, or 3, digits
1124                 if (e <= 9) {
1125                     if (form != Form.COMPATIBLE)
1126                         result[i++] = '0';
1127                     result[i++] = (char)(e+'0');
1128                 } else if (e <= 99) {
1129                     result[i++] = (char)(e/10 +'0');
1130                     result[i++] = (char)(e%10 + '0');
1131                 } else {
1132                     result[i++] = (char)(e/100+'0');
1133                     e %= 100;
1134                     result[i++] = (char)(e/10+'0');
1135                     result[i++] = (char)(e%10 + '0');
1136                 }
1137             }
1138         }
1139         return i;
1140     }
1141 
1142     // Per-thread buffer for string/stringbuffer conversion
1143     private static ThreadLocal perThreadBuffer = new ThreadLocal() {
1144             protected synchronized Object initialValue() {
1145                 return new char[26];
1146             }
1147         };
1148 
1149     /*
1150      * Take a FormattedFloatingDecimal, which we presumably just scanned in,
1151      * and find out what its value is, as a double.
1152      *
1153      * AS A SIDE EFFECT, SET roundDir TO INDICATE PREFERRED
1154      * ROUNDING DIRECTION in case the result is really destined
1155      * for a single-precision float.
1156      */
1157 
1158     public strictfp double doubleValue(){
1159         int     kDigits = Math.min( nDigits, maxDecimalDigits+1 );
1160         long    lValue;
1161         double  dValue;
1162         double  rValue, tValue;
1163 
1164         // First, check for NaN and Infinity values
1165         if(digits == infinity || digits == notANumber) {
1166             if(digits == notANumber)
1167                 return Double.NaN;
1168             else
1169                 return (isNegative?Double.NEGATIVE_INFINITY:Double.POSITIVE_INFINITY);
1170         }
1171         else {
1172             if (mustSetRoundDir) {
1173                 roundDir = 0;
1174             }
1175             /*
1176              * convert the lead kDigits to a long integer.
1177              */
1178             // (special performance hack: start to do it using int)
1179             int iValue = (int)digits[0]-(int)'0';
1180             int iDigits = Math.min( kDigits, intDecimalDigits );
1181             for ( int i=1; i < iDigits; i++ ){
1182                 iValue = iValue*10 + (int)digits[i]-(int)'0';
1183             }
1184             lValue = (long)iValue;
1185             for ( int i=iDigits; i < kDigits; i++ ){
1186                 lValue = lValue*10L + (long)((int)digits[i]-(int)'0');
1187             }
1188             dValue = (double)lValue;
1189             int exp = decExponent-kDigits;
1190             /*
1191              * lValue now contains a long integer with the value of
1192              * the first kDigits digits of the number.
1193              * dValue contains the (double) of the same.
1194              */
1195 
1196             if ( nDigits <= maxDecimalDigits ){
1197                 /*
1198                  * possibly an easy case.
1199                  * We know that the digits can be represented
1200                  * exactly. And if the exponent isn't too outrageous,
1201                  * the whole thing can be done with one operation,
1202                  * thus one rounding error.
1203                  * Note that all our constructors trim all leading and
1204                  * trailing zeros, so simple values (including zero)
1205                  * will always end up here
1206                  */
1207                 if (exp == 0 || dValue == 0.0)
1208                     return (isNegative)? -dValue : dValue; // small floating integer
1209                 else if ( exp >= 0 ){
1210                     if ( exp <= maxSmallTen ){
1211                         /*
1212                          * Can get the answer with one operation,
1213                          * thus one roundoff.
1214                          */
1215                         rValue = dValue * small10pow[exp];
1216                         if ( mustSetRoundDir ){
1217                             tValue = rValue / small10pow[exp];
1218                             roundDir = ( tValue ==  dValue ) ? 0
1219                                 :( tValue < dValue ) ? 1
1220                                 : -1;
1221                         }
1222                         return (isNegative)? -rValue : rValue;
1223                     }
1224                     int slop = maxDecimalDigits - kDigits;
1225                     if ( exp <= maxSmallTen+slop ){
1226                         /*
1227                          * We can multiply dValue by 10^(slop)
1228                          * and it is still "small" and exact.
1229                          * Then we can multiply by 10^(exp-slop)
1230                          * with one rounding.
1231                          */
1232                         dValue *= small10pow[slop];
1233                         rValue = dValue * small10pow[exp-slop];
1234 
1235                         if ( mustSetRoundDir ){
1236                             tValue = rValue / small10pow[exp-slop];
1237                             roundDir = ( tValue ==  dValue ) ? 0
1238                                 :( tValue < dValue ) ? 1
1239                                 : -1;
1240                         }
1241                         return (isNegative)? -rValue : rValue;
1242                     }
1243                     /*
1244                      * Else we have a hard case with a positive exp.
1245                      */
1246                 } else {
1247                     if ( exp >= -maxSmallTen ){
1248                         /*
1249                          * Can get the answer in one division.
1250                          */
1251                         rValue = dValue / small10pow[-exp];
1252                         tValue = rValue * small10pow[-exp];
1253                         if ( mustSetRoundDir ){
1254                             roundDir = ( tValue ==  dValue ) ? 0
1255                                 :( tValue < dValue ) ? 1
1256                                 : -1;
1257                         }
1258                         return (isNegative)? -rValue : rValue;
1259                     }
1260                     /*
1261                      * Else we have a hard case with a negative exp.
1262                      */
1263                 }
1264             }
1265 
1266             /*
1267              * Harder cases:
1268              * The sum of digits plus exponent is greater than
1269              * what we think we can do with one error.
1270              *
1271              * Start by approximating the right answer by,
1272              * naively, scaling by powers of 10.
1273              */
1274             if ( exp > 0 ){
1275                 if ( decExponent > maxDecimalExponent+1 ){
1276                     /*
1277                      * Lets face it. This is going to be
1278                      * Infinity. Cut to the chase.
1279                      */
1280                     return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
1281                 }
1282                 if ( (exp&15) != 0 ){
1283                     dValue *= small10pow[exp&15];
1284                 }
1285                 if ( (exp>>=4) != 0 ){
1286                     int j;
1287                     for( j = 0; exp > 1; j++, exp>>=1 ){
1288                         if ( (exp&1)!=0)
1289                             dValue *= big10pow[j];
1290                     }
1291                     /*
1292                      * The reason for the weird exp > 1 condition
1293                      * in the above loop was so that the last multiply
1294                      * would get unrolled. We handle it here.
1295                      * It could overflow.
1296                      */
1297                     double t = dValue * big10pow[j];
1298                     if ( Double.isInfinite( t ) ){
1299                         /*
1300                          * It did overflow.
1301                          * Look more closely at the result.
1302                          * If the exponent is just one too large,
1303                          * then use the maximum finite as our estimate
1304                          * value. Else call the result infinity
1305                          * and punt it.
1306                          * ( I presume this could happen because
1307                          * rounding forces the result here to be
1308                          * an ULP or two larger than
1309                          * Double.MAX_VALUE ).
1310                          */
1311                         t = dValue / 2.0;
1312                         t *= big10pow[j];
1313                         if ( Double.isInfinite( t ) ){
1314                             return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
1315                         }
1316                         t = Double.MAX_VALUE;
1317                     }
1318                     dValue = t;
1319                 }
1320             } else if ( exp < 0 ){
1321                 exp = -exp;
1322                 if ( decExponent < minDecimalExponent-1 ){
1323                     /*
1324                      * Lets face it. This is going to be
1325                      * zero. Cut to the chase.
1326                      */
1327                     return (isNegative)? -0.0 : 0.0;
1328                 }
1329                 if ( (exp&15) != 0 ){
1330                     dValue /= small10pow[exp&15];
1331                 }
1332                 if ( (exp>>=4) != 0 ){
1333                     int j;
1334                     for( j = 0; exp > 1; j++, exp>>=1 ){
1335                         if ( (exp&1)!=0)
1336                             dValue *= tiny10pow[j];
1337                     }
1338                     /*
1339                      * The reason for the weird exp > 1 condition
1340                      * in the above loop was so that the last multiply
1341                      * would get unrolled. We handle it here.
1342                      * It could underflow.
1343                      */
1344                     double t = dValue * tiny10pow[j];
1345                     if ( t == 0.0 ){
1346                         /*
1347                          * It did underflow.
1348                          * Look more closely at the result.
1349                          * If the exponent is just one too small,
1350                          * then use the minimum finite as our estimate
1351                          * value. Else call the result 0.0
1352                          * and punt it.
1353                          * ( I presume this could happen because
1354                          * rounding forces the result here to be
1355                          * an ULP or two less than
1356                          * Double.MIN_VALUE ).
1357                          */
1358                         t = dValue * 2.0;
1359                         t *= tiny10pow[j];
1360                         if ( t == 0.0 ){
1361                             return (isNegative)? -0.0 : 0.0;
1362                         }
1363                         t = Double.MIN_VALUE;
1364                     }
1365                     dValue = t;
1366                 }
1367             }
1368 
1369             /*
1370              * dValue is now approximately the result.
1371              * The hard part is adjusting it, by comparison
1372              * with FDBigInt arithmetic.
1373              * Formulate the EXACT big-number result as
1374              * bigD0 * 10^exp
1375              */
1376             FDBigInt bigD0 = new FDBigInt( lValue, digits, kDigits, nDigits );
1377             exp   = decExponent - nDigits;
1378 
1379             correctionLoop:
1380             while(true){
1381                 /* AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES
1382                  * bigIntExp and bigIntNBits
1383                  */
1384                 FDBigInt bigB = doubleToBigInt( dValue );
1385 
1386                 /*
1387                  * Scale bigD, bigB appropriately for
1388                  * big-integer operations.
1389                  * Naively, we multipy by powers of ten
1390                  * and powers of two. What we actually do
1391                  * is keep track of the powers of 5 and
1392                  * powers of 2 we would use, then factor out
1393                  * common divisors before doing the work.
1394                  */
1395                 int B2, B5; // powers of 2, 5 in bigB
1396                 int     D2, D5; // powers of 2, 5 in bigD
1397                 int Ulp2;   // powers of 2 in halfUlp.
1398                 if ( exp >= 0 ){
1399                     B2 = B5 = 0;
1400                     D2 = D5 = exp;
1401                 } else {
1402                     B2 = B5 = -exp;
1403                     D2 = D5 = 0;
1404                 }
1405                 if ( bigIntExp >= 0 ){
1406                     B2 += bigIntExp;
1407                 } else {
1408                     D2 -= bigIntExp;
1409                 }
1410                 Ulp2 = B2;
1411                 // shift bigB and bigD left by a number s. t.
1412                 // halfUlp is still an integer.
1413                 int hulpbias;
1414                 if ( bigIntExp+bigIntNBits <= -expBias+1 ){
1415                     // This is going to be a denormalized number
1416                     // (if not actually zero).
1417                     // half an ULP is at 2^-(expBias+expShift+1)
1418                     hulpbias = bigIntExp+ expBias + expShift;
1419                 } else {
1420                     hulpbias = expShift + 2 - bigIntNBits;
1421                 }
1422                 B2 += hulpbias;
1423                 D2 += hulpbias;
1424                 // if there are common factors of 2, we might just as well
1425                 // factor them out, as they add nothing useful.
1426                 int common2 = Math.min( B2, Math.min( D2, Ulp2 ) );
1427                 B2 -= common2;
1428                 D2 -= common2;
1429                 Ulp2 -= common2;
1430                 // do multiplications by powers of 5 and 2
1431                 bigB = multPow52( bigB, B5, B2 );
1432                 FDBigInt bigD = multPow52( new FDBigInt( bigD0 ), D5, D2 );
1433                 //
1434                 // to recap:
1435                 // bigB is the scaled-big-int version of our floating-point
1436                 // candidate.
1437                 // bigD is the scaled-big-int version of the exact value
1438                 // as we understand it.
1439                 // halfUlp is 1/2 an ulp of bigB, except for special cases
1440                 // of exact powers of 2
1441                 //
1442                 // the plan is to compare bigB with bigD, and if the difference
1443                 // is less than halfUlp, then we're satisfied. Otherwise,
1444                 // use the ratio of difference to halfUlp to calculate a fudge
1445                 // factor to add to the floating value, then go 'round again.
1446                 //
1447                 FDBigInt diff;
1448                 int cmpResult;
1449                 boolean overvalue;
1450                 if ( (cmpResult = bigB.cmp( bigD ) ) > 0 ){
1451                     overvalue = true; // our candidate is too big.
1452                     diff = bigB.sub( bigD );
1453                     if ( (bigIntNBits == 1) && (bigIntExp > -expBias) ){
1454                         // candidate is a normalized exact power of 2 and
1455                         // is too big. We will be subtracting.
1456                         // For our purposes, ulp is the ulp of the
1457                         // next smaller range.
1458                         Ulp2 -= 1;
1459                         if ( Ulp2 < 0 ){
1460                             // rats. Cannot de-scale ulp this far.
1461                             // must scale diff in other direction.
1462                             Ulp2 = 0;
1463                             diff.lshiftMe( 1 );
1464                         }
1465                     }
1466                 } else if ( cmpResult < 0 ){
1467                     overvalue = false; // our candidate is too small.
1468                     diff = bigD.sub( bigB );
1469                 } else {
1470                     // the candidate is exactly right!
1471                     // this happens with surprising fequency
1472                     break correctionLoop;
1473                 }
1474                 FDBigInt halfUlp = constructPow52( B5, Ulp2 );
1475                 if ( (cmpResult = diff.cmp( halfUlp ) ) < 0 ){
1476                     // difference is small.
1477                     // this is close enough
1478                     if (mustSetRoundDir) {
1479                         roundDir = overvalue ? -1 : 1;
1480                     }
1481                     break correctionLoop;
1482                 } else if ( cmpResult == 0 ){
1483                     // difference is exactly half an ULP
1484                     // round to some other value maybe, then finish
1485                     dValue += 0.5*ulp( dValue, overvalue );
1486                     // should check for bigIntNBits == 1 here??
1487                     if (mustSetRoundDir) {
1488                         roundDir = overvalue ? -1 : 1;
1489                     }
1490                     break correctionLoop;
1491                 } else {
1492                     // difference is non-trivial.
1493                     // could scale addend by ratio of difference to
1494                     // halfUlp here, if we bothered to compute that difference.
1495                     // Most of the time ( I hope ) it is about 1 anyway.
1496                     dValue += ulp( dValue, overvalue );
1497                     if ( dValue == 0.0 || dValue == Double.POSITIVE_INFINITY )
1498                         break correctionLoop; // oops. Fell off end of range.
1499                     continue; // try again.
1500                 }
1501 
1502             }
1503             return (isNegative)? -dValue : dValue;
1504         }
1505     }
1506 
1507     /*
1508      * Take a FormattedFloatingDecimal, which we presumably just scanned in,
1509      * and find out what its value is, as a float.
1510      * This is distinct from doubleValue() to avoid the extremely
1511      * unlikely case of a double rounding error, wherein the converstion
1512      * to double has one rounding error, and the conversion of that double
1513      * to a float has another rounding error, IN THE WRONG DIRECTION,
1514      * ( because of the preference to a zero low-order bit ).
1515      */
1516 
1517     public strictfp float floatValue(){
1518         int     kDigits = Math.min( nDigits, singleMaxDecimalDigits+1 );
1519         int     iValue;
1520         float   fValue;
1521 
1522         // First, check for NaN and Infinity values
1523         if(digits == infinity || digits == notANumber) {
1524             if(digits == notANumber)
1525                 return Float.NaN;
1526             else
1527                 return (isNegative?Float.NEGATIVE_INFINITY:Float.POSITIVE_INFINITY);
1528         }
1529         else {
1530             /*
1531              * convert the lead kDigits to an integer.
1532              */
1533             iValue = (int)digits[0]-(int)'0';
1534             for ( int i=1; i < kDigits; i++ ){
1535                 iValue = iValue*10 + (int)digits[i]-(int)'0';
1536             }
1537             fValue = (float)iValue;
1538             int exp = decExponent-kDigits;
1539             /*
1540              * iValue now contains an integer with the value of
1541              * the first kDigits digits of the number.
1542              * fValue contains the (float) of the same.
1543              */
1544 
1545             if ( nDigits <= singleMaxDecimalDigits ){
1546                 /*
1547                  * possibly an easy case.
1548                  * We know that the digits can be represented
1549                  * exactly. And if the exponent isn't too outrageous,
1550                  * the whole thing can be done with one operation,
1551                  * thus one rounding error.
1552                  * Note that all our constructors trim all leading and
1553                  * trailing zeros, so simple values (including zero)
1554                  * will always end up here.
1555                  */
1556                 if (exp == 0 || fValue == 0.0f)
1557                     return (isNegative)? -fValue : fValue; // small floating integer
1558                 else if ( exp >= 0 ){
1559                     if ( exp <= singleMaxSmallTen ){
1560                         /*
1561                          * Can get the answer with one operation,
1562                          * thus one roundoff.
1563                          */
1564                         fValue *= singleSmall10pow[exp];
1565                         return (isNegative)? -fValue : fValue;
1566                     }
1567                     int slop = singleMaxDecimalDigits - kDigits;
1568                     if ( exp <= singleMaxSmallTen+slop ){
1569                         /*
1570                          * We can multiply dValue by 10^(slop)
1571                          * and it is still "small" and exact.
1572                          * Then we can multiply by 10^(exp-slop)
1573                          * with one rounding.
1574                          */
1575                         fValue *= singleSmall10pow[slop];
1576                         fValue *= singleSmall10pow[exp-slop];
1577                         return (isNegative)? -fValue : fValue;
1578                     }
1579                     /*
1580                      * Else we have a hard case with a positive exp.
1581                      */
1582                 } else {
1583                     if ( exp >= -singleMaxSmallTen ){
1584                         /*
1585                          * Can get the answer in one division.
1586                          */
1587                         fValue /= singleSmall10pow[-exp];
1588                         return (isNegative)? -fValue : fValue;
1589                     }
1590                     /*
1591                      * Else we have a hard case with a negative exp.
1592                      */
1593                 }
1594             } else if ( (decExponent >= nDigits) && (nDigits+decExponent <= maxDecimalDigits) ){
1595                 /*
1596                  * In double-precision, this is an exact floating integer.
1597                  * So we can compute to double, then shorten to float
1598                  * with one round, and get the right answer.
1599                  *
1600                  * First, finish accumulating digits.
1601                  * Then convert that integer to a double, multiply
1602                  * by the appropriate power of ten, and convert to float.
1603                  */
1604                 long lValue = (long)iValue;
1605                 for ( int i=kDigits; i < nDigits; i++ ){
1606                     lValue = lValue*10L + (long)((int)digits[i]-(int)'0');
1607                 }
1608                 double dValue = (double)lValue;
1609                 exp = decExponent-nDigits;
1610                 dValue *= small10pow[exp];
1611                 fValue = (float)dValue;
1612                 return (isNegative)? -fValue : fValue;
1613 
1614             }
1615             /*
1616              * Harder cases:
1617              * The sum of digits plus exponent is greater than
1618              * what we think we can do with one error.
1619              *
1620              * Start by weeding out obviously out-of-range
1621              * results, then convert to double and go to
1622              * common hard-case code.
1623              */
1624             if ( decExponent > singleMaxDecimalExponent+1 ){
1625                 /*
1626                  * Lets face it. This is going to be
1627                  * Infinity. Cut to the chase.
1628                  */
1629                 return (isNegative)? Float.NEGATIVE_INFINITY : Float.POSITIVE_INFINITY;
1630             } else if ( decExponent < singleMinDecimalExponent-1 ){
1631                 /*
1632                  * Lets face it. This is going to be
1633                  * zero. Cut to the chase.
1634                  */
1635                 return (isNegative)? -0.0f : 0.0f;
1636             }
1637 
1638             /*
1639              * Here, we do 'way too much work, but throwing away
1640              * our partial results, and going and doing the whole
1641              * thing as double, then throwing away half the bits that computes
1642              * when we convert back to float.
1643              *
1644              * The alternative is to reproduce the whole multiple-precision
1645              * algorythm for float precision, or to try to parameterize it
1646              * for common usage. The former will take about 400 lines of code,
1647              * and the latter I tried without success. Thus the semi-hack
1648              * answer here.
1649              */
1650             mustSetRoundDir = !fromHex;
1651             double dValue = doubleValue();
1652             return stickyRound( dValue );
1653         }
1654     }
1655 
1656 
1657     /*
1658      * All the positive powers of 10 that can be
1659      * represented exactly in double/float.
1660      */
1661     private static final double small10pow[] = {
1662         1.0e0,
1663         1.0e1, 1.0e2, 1.0e3, 1.0e4, 1.0e5,
1664         1.0e6, 1.0e7, 1.0e8, 1.0e9, 1.0e10,
1665         1.0e11, 1.0e12, 1.0e13, 1.0e14, 1.0e15,
1666         1.0e16, 1.0e17, 1.0e18, 1.0e19, 1.0e20,
1667         1.0e21, 1.0e22
1668     };
1669 
1670     private static final float singleSmall10pow[] = {
1671         1.0e0f,
1672         1.0e1f, 1.0e2f, 1.0e3f, 1.0e4f, 1.0e5f,
1673         1.0e6f, 1.0e7f, 1.0e8f, 1.0e9f, 1.0e10f
1674     };
1675 
1676     private static final double big10pow[] = {
1677         1e16, 1e32, 1e64, 1e128, 1e256 };
1678     private static final double tiny10pow[] = {
1679         1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1680 
1681     private static final int maxSmallTen = small10pow.length-1;
1682     private static final int singleMaxSmallTen = singleSmall10pow.length-1;
1683 
1684     private static final int small5pow[] = {
1685         1,
1686         5,
1687         5*5,
1688         5*5*5,
1689         5*5*5*5,
1690         5*5*5*5*5,
1691         5*5*5*5*5*5,
1692         5*5*5*5*5*5*5,
1693         5*5*5*5*5*5*5*5,
1694         5*5*5*5*5*5*5*5*5,
1695         5*5*5*5*5*5*5*5*5*5,
1696         5*5*5*5*5*5*5*5*5*5*5,
1697         5*5*5*5*5*5*5*5*5*5*5*5,
1698         5*5*5*5*5*5*5*5*5*5*5*5*5
1699     };
1700 
1701 
1702     private static final long long5pow[] = {
1703         1L,
1704         5L,
1705         5L*5,
1706         5L*5*5,
1707         5L*5*5*5,
1708         5L*5*5*5*5,
1709         5L*5*5*5*5*5,
1710         5L*5*5*5*5*5*5,
1711         5L*5*5*5*5*5*5*5,
1712         5L*5*5*5*5*5*5*5*5,
1713         5L*5*5*5*5*5*5*5*5*5,
1714         5L*5*5*5*5*5*5*5*5*5*5,
1715         5L*5*5*5*5*5*5*5*5*5*5*5,
1716         5L*5*5*5*5*5*5*5*5*5*5*5*5,
1717         5L*5*5*5*5*5*5*5*5*5*5*5*5*5,
1718         5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1719         5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1720         5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1721         5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1722         5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1723         5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1724         5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1725         5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1726         5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1727         5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1728         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,
1729         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,
1730     };
1731 
1732     // approximately ceil( log2( long5pow[i] ) )
1733     private static final int n5bits[] = {
1734         0,
1735         3,
1736         5,
1737         7,
1738         10,
1739         12,
1740         14,
1741         17,
1742         19,
1743         21,
1744         24,
1745         26,
1746         28,
1747         31,
1748         33,
1749         35,
1750         38,
1751         40,
1752         42,
1753         45,
1754         47,
1755         49,
1756         52,
1757         54,
1758         56,
1759         59,
1760         61,
1761     };
1762 
1763     private static final char infinity[] = { 'I', 'n', 'f', 'i', 'n', 'i', 't', 'y' };
1764     private static final char notANumber[] = { 'N', 'a', 'N' };
1765     private static final char zero[] = { '0', '0', '0', '0', '0', '0', '0', '0' };
1766 }
--- EOF ---