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

Save This Page
Home » openjdk-7 » sun » misc » [javadoc | source]