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 }