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

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