home *** CD-ROM | disk | FTP | other *** search
/ Apple Developer Connection Student Program / ADC Tools Sampler CD Disk 3 1999.iso / Metrowerks CodeWarrior / Java Support / Java_Source / Java2 / src / java / lang / FloatingDecimal.java < prev    next >
Encoding:
Java Source  |  1999-05-28  |  61.5 KB  |  2,248 lines  |  [TEXT/CWIE]
  1. /*
  2.  * @(#)FloatingDecimal.java    1.14 98/07/07
  3.  *
  4.  * Copyright 1996-1998 by Sun Microsystems, Inc.,
  5.  * 901 San Antonio Road, Palo Alto, California, 94303, U.S.A.
  6.  * All rights reserved.
  7.  *
  8.  * This software is the confidential and proprietary information
  9.  * of Sun Microsystems, Inc. ("Confidential Information").  You
  10.  * shall not disclose such Confidential Information and shall use
  11.  * it only in accordance with the terms of the license agreement
  12.  * you entered into with Sun.
  13.  */
  14.  
  15. package java.lang;
  16.  
  17. class FloatingDecimal{
  18.     boolean    isExceptional;
  19.     boolean    isNegative;
  20.     int        decExponent;
  21.     char    digits[];
  22.     int        nDigits;
  23.     int        bigIntExp;
  24.     int        bigIntNBits;
  25.     boolean    mustSetRoundDir = false;
  26.     int        roundDir; // set by doubleValue
  27.  
  28.     private    FloatingDecimal( boolean negSign, int decExponent, char []digits, int n,  boolean e )
  29.     {
  30.     isNegative = negSign;
  31.     isExceptional = e;
  32.     this.decExponent = decExponent;
  33.     this.digits = digits;
  34.     this.nDigits = n;
  35.     }
  36.  
  37.     /*
  38.      * Constants of the implementation
  39.      * Most are IEEE-754 related.
  40.      * (There are more really boring constants at the end.)
  41.      */
  42.     static final long    signMask = 0x8000000000000000L;
  43.     static final long    expMask  = 0x7ff0000000000000L;
  44.     static final long    fractMask= ~(signMask|expMask);
  45.     static final int    expShift = 52;
  46.     static final int    expBias  = 1023;
  47.     static final long    fractHOB = ( 1L<<expShift ); // assumed High-Order bit
  48.     static final long    expOne     = ((long)expBias)<<expShift; // exponent of 1.0
  49.     static final int    maxSmallBinExp = 62;
  50.     static final int    minSmallBinExp = -( 63 / 3 );
  51.     static final int    maxDecimalDigits = 15;
  52.     static final int    maxDecimalExponent = 308;
  53.     static final int    minDecimalExponent = -324;
  54.     static final int    bigDecimalExponent = 324; // i.e. abs(minDecimalExponent)
  55.  
  56.     static final long    highbyte = 0xff00000000000000L;
  57.     static final long    highbit  = 0x8000000000000000L;
  58.     static final long    lowbytes = ~highbyte;
  59.  
  60.     static final int    singleSignMask =    0x80000000;
  61.     static final int    singleExpMask  =    0x7f800000;
  62.     static final int    singleFractMask =   ~(singleSignMask|singleExpMask);
  63.     static final int    singleExpShift    =   23;
  64.     static final int    singleFractHOB    =   1<<singleExpShift;
  65.     static final int    singleExpBias    =   127;
  66.     static final int    singleMaxDecimalDigits = 7;
  67.     static final int    singleMaxDecimalExponent = 38;
  68.     static final int    singleMinDecimalExponent = -45;
  69.  
  70.     static final int    intDecimalDigits = 9;
  71.  
  72.  
  73.     /*
  74.      * count number of bits from high-order 1 bit to low-order 1 bit,
  75.      * inclusive.
  76.      */
  77.     private static int
  78.     countBits( long v ){
  79.     //
  80.     // the strategy is to shift until we get a non-zero sign bit
  81.     // then shift until we have no bits left, counting the difference.
  82.     // we do byte shifting as a hack. Hope it helps.
  83.     //
  84.     if ( v == 0L ) return 0;
  85.  
  86.     while ( ( v & highbyte ) == 0L ){
  87.         v <<= 8;
  88.     }
  89.     while ( v > 0L ) { // i.e. while ((v&highbit) == 0L )
  90.         v <<= 1;
  91.     }
  92.  
  93.     int n = 0;
  94.     while (( v & lowbytes ) != 0L ){
  95.         v <<= 8;
  96.         n += 8;
  97.     }
  98.     while ( v != 0L ){
  99.         v <<= 1;
  100.         n += 1;
  101.     }
  102.     return n;
  103.     }
  104.  
  105.     /*
  106.      * Keep big powers of 5 handy for future reference.
  107.      */
  108.     private static FDBigInt b5p[];
  109.  
  110.     private static synchronized FDBigInt
  111.     big5pow( int p ){
  112.     if ( p < 0 )
  113.         throw new RuntimeException( "Assertion botch: negative power of 5");
  114.     if ( b5p == null ){
  115.         b5p = new FDBigInt[ p+1 ];
  116.     }else if (b5p.length <= p ){
  117.         FDBigInt t[] = new FDBigInt[ p+1 ];
  118.         System.arraycopy( b5p, 0, t, 0, b5p.length );
  119.         b5p = t;
  120.     }
  121.     if ( b5p[p] != null )
  122.         return b5p[p];
  123.     else if ( p < small5pow.length )
  124.         return b5p[p] = new FDBigInt( small5pow[p] );
  125.     else if ( p < long5pow.length )
  126.         return b5p[p] = new FDBigInt( long5pow[p] );
  127.     else {
  128.         // construct the value.
  129.         // recursively.
  130.         int q, r;
  131.         // in order to compute 5^p,
  132.         // compute its square root, 5^(p/2) and square.
  133.         // or, let q = p / 2, r = p -q, then
  134.         // 5^p = 5^(q+r) = 5^q * 5^r
  135.         q = p >> 1;
  136.         r = p - q;
  137.         FDBigInt bigq =  b5p[q];
  138.         if ( bigq == null )
  139.         bigq = big5pow ( q );
  140.         if ( r < small5pow.length ){
  141.         return (b5p[p] = bigq.mult( small5pow[r] ) );
  142.         }else{
  143.         FDBigInt bigr = b5p[ r ];
  144.         if ( bigr == null )
  145.             bigr = big5pow( r );
  146.         return (b5p[p] = bigq.mult( bigr ) );
  147.         }
  148.     }
  149.     }
  150.  
  151.     //
  152.     // a common operation
  153.     //
  154.     private static FDBigInt
  155.     multPow52( FDBigInt v, int p5, int p2 ){
  156.     if ( p5 != 0 ){
  157.         if ( p5 < small5pow.length ){
  158.         v = v.mult( small5pow[p5] );
  159.         } else {
  160.         v = v.mult( big5pow( p5 ) );
  161.         }
  162.     }
  163.     if ( p2 != 0 ){
  164.         v.lshiftMe( p2 );
  165.     }
  166.     return v;
  167.     }
  168.  
  169.     //
  170.     // another common operation
  171.     //
  172.     private static FDBigInt
  173.     constructPow52( int p5, int p2 ){
  174.     FDBigInt v = new FDBigInt( big5pow( p5 ) );
  175.     if ( p2 != 0 ){
  176.         v.lshiftMe( p2 );
  177.     }
  178.     return v;
  179.     }
  180.  
  181.     /*
  182.      * Make a floating double into a FDBigInt.
  183.      * This could also be structured as a FDBigInt
  184.      * constructor, but we'd have to build a lot of knowledge
  185.      * about floating-point representation into it, and we don't want to.
  186.      *
  187.      * AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES
  188.      * bigIntExp and bigIntNBits
  189.      *
  190.      */
  191.     private FDBigInt
  192.     doubleToBigInt( double dval ){
  193.     long lbits = Double.doubleToLongBits( dval ) & ~signMask;
  194.     int binexp = (int)(lbits >>> expShift);
  195.     lbits &= fractMask;
  196.     if ( binexp > 0 ){
  197.         lbits |= fractHOB;
  198.     } else {
  199.         if ( lbits == 0L )
  200.         throw new RuntimeException("Assertion botch: doubleToBigInt(0.0)");
  201.         binexp +=1;
  202.         while ( (lbits & fractHOB ) == 0L){
  203.         lbits <<= 1;
  204.         binexp -= 1;
  205.         }
  206.     }
  207.     binexp -= expBias;
  208.     int nbits = countBits( lbits );
  209.     /*
  210.      * We now know where the high-order 1 bit is,
  211.      * and we know how many there are.
  212.      */
  213.     int lowOrderZeros = expShift+1-nbits;
  214.     lbits >>>= lowOrderZeros;
  215.  
  216.     bigIntExp = binexp+1-nbits;
  217.     bigIntNBits = nbits;
  218.     return new FDBigInt( lbits );
  219.     }
  220.  
  221.     /*
  222.      * Compute a number that is the ULP of the given value,
  223.      * for purposes of addition/subtraction. Generally easy.
  224.      * More difficult if subtracting and the argument
  225.      * is a normalized a power of 2, as the ULP changes at these points.
  226.      */
  227.     private static double
  228.     ulp( double dval, boolean subtracting ){
  229.     long lbits = Double.doubleToLongBits( dval ) & ~signMask;
  230.     int binexp = (int)(lbits >>> expShift);
  231.     double ulpval;
  232.     if ( subtracting && ( binexp >= expShift ) && ((lbits&fractMask) == 0L) ){
  233.         // for subtraction from normalized, powers of 2,
  234.         // use next-smaller exponent
  235.         binexp -= 1;
  236.     }
  237.     if ( binexp > expShift ){
  238.         ulpval = Double.longBitsToDouble( ((long)(binexp-expShift))<<expShift );
  239.     } else if ( binexp == 0 ){
  240.         ulpval = Double.MIN_VALUE;
  241.     } else {
  242.         ulpval = Double.longBitsToDouble( 1L<<(binexp-1) );
  243.     }
  244.     if ( subtracting ) ulpval = - ulpval;
  245.  
  246.     return ulpval;
  247.     }
  248.  
  249.     /*
  250.      * Round a double to a float.
  251.      * In addition to the fraction bits of the double,
  252.      * look at the class instance variable roundDir,
  253.      * which should help us avoid double-rounding error.
  254.      * roundDir was set in hardValueOf if the estimate was
  255.      * close enough, but not exact. It tells us which direction
  256.      * of rounding is preferred.
  257.      */
  258.     float
  259.     stickyRound( double dval ){
  260.     long lbits = Double.doubleToLongBits( dval );
  261.     long binexp = lbits & expMask;
  262.     if ( binexp == 0L || binexp == expMask ){
  263.         // what we have here is special.
  264.         // don't worry, the right thing will happen.
  265.         return (float) dval;
  266.     }
  267.     lbits += (long)roundDir; // hack-o-matic.
  268.     return (float)Double.longBitsToDouble( lbits );
  269.     }
  270.  
  271.  
  272.     /*
  273.      * This is the easy subcase --
  274.      * all the significant bits, after scaling, are held in lvalue.
  275.      * negSign and decExponent tell us what processing and scaling
  276.      * has already been done. Exceptional cases have already been
  277.      * stripped out.
  278.      * In particular:
  279.      * lvalue is a finite number (not Inf, nor NaN)
  280.      * lvalue > 0L (not zero, nor negative).
  281.      *
  282.      * The only reason that we develop the digits here, rather than
  283.      * calling on Long.toString() is that we can do it a little faster,
  284.      * and besides want to treat trailing 0s specially. If Long.toString
  285.      * changes, we should re-evaluate this strategy!
  286.      */
  287.     private void
  288.     developLongDigits( int decExponent, long lvalue, long insignificant ){
  289.     char digits[];
  290.     int  ndigits;
  291.     int  digitno;
  292.     int  c;
  293.     //
  294.     // Discard non-significant low-order bits, while rounding,
  295.     // up to insignificant value.
  296.     int i;
  297.     for ( i = 0; insignificant >= 10L; i++ )
  298.         insignificant /= 10L;
  299.     if ( i != 0 ){
  300.         long pow10 = long5pow[i] << i; // 10^i == 5^i * 2^i;
  301.         long residue = lvalue % pow10;
  302.         lvalue /= pow10;
  303.         decExponent += i;
  304.         if ( residue >= (pow10>>1) ){
  305.         // round up based on the low-order bits we're discarding
  306.         lvalue++;
  307.         }
  308.     }
  309.     if ( lvalue <= Integer.MAX_VALUE ){
  310.         if ( lvalue <= 0L )
  311.         throw new RuntimeException("Assertion botch: value "+lvalue+" <= 0");
  312.  
  313.         // even easier subcase!
  314.         // can do int arithmetic rather than long!
  315.         int  ivalue = (int)lvalue;
  316.         digits = new char[ ndigits=10 ];
  317.         digitno = ndigits-1;
  318.         c = ivalue%10;
  319.         ivalue /= 10;
  320.         while ( c == 0 ){
  321.         decExponent++;
  322.         c = ivalue%10;
  323.         ivalue /= 10;
  324.         }
  325.         while ( ivalue != 0){
  326.         digits[digitno--] = (char)(c+'0');
  327.         decExponent++;
  328.         c = ivalue%10;
  329.         ivalue /= 10;
  330.         }
  331.         digits[digitno] = (char)(c+'0');
  332.     } else {
  333.         // same algorithm as above (same bugs, too )
  334.         // but using long arithmetic.
  335.         digits = new char[ ndigits=20 ];
  336.         digitno = ndigits-1;
  337.         c = (int)(lvalue%10L);
  338.         lvalue /= 10L;
  339.         while ( c == 0 ){
  340.         decExponent++;
  341.         c = (int)(lvalue%10L);
  342.         lvalue /= 10L;
  343.         }
  344.         while ( lvalue != 0L ){
  345.         digits[digitno--] = (char)(c+'0');
  346.         decExponent++;
  347.         c = (int)(lvalue%10L);
  348.         lvalue /= 10;
  349.         }
  350.         digits[digitno] = (char)(c+'0');
  351.     }
  352.     char result [];
  353.     ndigits -= digitno;
  354.     if ( digitno == 0 )
  355.         result = digits;
  356.     else {
  357.         result = new char[ ndigits ];
  358.         System.arraycopy( digits, digitno, result, 0, ndigits );
  359.     }
  360.     this.digits = result;
  361.     this.decExponent = decExponent+1;
  362.     this.nDigits = ndigits;
  363.     }
  364.  
  365.     //
  366.     // add one to the least significant digit.
  367.     // in the unlikely event there is a carry out,
  368.     // deal with it.
  369.     // assert that this will only happen where there
  370.     // is only one digit, e.g. (float)1e-44 seems to do it.
  371.     //
  372.     private void
  373.     roundup(){
  374.     int i;
  375.     int q = digits[ i = (nDigits-1)];
  376.     if ( q == '9' ){
  377.         while ( q == '9' && i > 0 ){
  378.         digits[i] = '0';
  379.         q = digits[--i];
  380.         }
  381.         if ( q == '9' ){
  382.         // carryout! High-order 1, rest 0s, larger exp.
  383.         decExponent += 1;
  384.         digits[0] = '1';
  385.         return;
  386.         }
  387.         // else fall through.
  388.     }
  389.     digits[i] = (char)(q+1);
  390.     }
  391.  
  392.     /*
  393.      * FIRST IMPORTANT CONSTRUCTOR: DOUBLE
  394.      */
  395.     public FloatingDecimal( double d )
  396.     {
  397.     long    dBits = Double.doubleToLongBits( d );
  398.     long    fractBits;
  399.     int    binExp;
  400.     int    nSignificantBits;
  401.  
  402.     // discover and delete sign
  403.     if ( (dBits&signMask) != 0 ){
  404.         isNegative = true;
  405.         dBits ^= signMask;
  406.     } else {
  407.         isNegative = false;
  408.     }
  409.     // Begin to unpack
  410.     // Discover obvious special cases of NaN and Infinity.
  411.     binExp = (int)( (dBits&expMask) >> expShift );
  412.     fractBits = dBits&fractMask;
  413.     if ( binExp == (int)(expMask>>expShift) ) {
  414.         isExceptional = true;
  415.         if ( fractBits == 0L ){
  416.         digits =  infinity;
  417.         } else {
  418.         digits = notANumber;
  419.         isNegative = false; // NaN has no sign!
  420.         }
  421.         nDigits = digits.length;
  422.         return;
  423.     }
  424.     isExceptional = false;
  425.     // Finish unpacking
  426.     // Normalize denormalized numbers.
  427.     // Insert assumed high-order bit for normalized numbers.
  428.     // Subtract exponent bias.
  429.     if ( binExp == 0 ){
  430.         if ( fractBits == 0L ){
  431.         // not a denorm, just a 0!
  432.         decExponent = 0;
  433.         digits = zero;
  434.         nDigits = 1;
  435.         return;
  436.         }
  437.         while ( (fractBits&fractHOB) == 0L ){
  438.         fractBits <<= 1;
  439.         binExp -= 1;
  440.         }
  441.         nSignificantBits = expShift + binExp +1; // recall binExp is  - shift count.
  442.         binExp += 1;
  443.     } else {
  444.         fractBits |= fractHOB;
  445.         nSignificantBits = expShift+1;
  446.     }
  447.     binExp -= expBias;
  448.     // call the routine that actually does all the hard work.
  449.     dtoa( binExp, fractBits, nSignificantBits );
  450.     }
  451.  
  452.     /*
  453.      * SECOND IMPORTANT CONSTRUCTOR: SINGLE
  454.      */
  455.     public FloatingDecimal( float f )
  456.     {
  457.     int    fBits = Float.floatToIntBits( f );
  458.     int    fractBits;
  459.     int    binExp;
  460.     int    nSignificantBits;
  461.  
  462.     // discover and delete sign
  463.     if ( (fBits&singleSignMask) != 0 ){
  464.         isNegative = true;
  465.         fBits ^= singleSignMask;
  466.     } else {
  467.         isNegative = false;
  468.     }
  469.     // Begin to unpack
  470.     // Discover obvious special cases of NaN and Infinity.
  471.     binExp = (int)( (fBits&singleExpMask) >> singleExpShift );
  472.     fractBits = fBits&singleFractMask;
  473.     if ( binExp == (int)(singleExpMask>>singleExpShift) ) {
  474.         isExceptional = true;
  475.         if ( fractBits == 0L ){
  476.         digits =  infinity;
  477.         } else {
  478.         digits = notANumber;
  479.         isNegative = false; // NaN has no sign!
  480.         }
  481.         nDigits = digits.length;
  482.         return;
  483.     }
  484.     isExceptional = false;
  485.     // Finish unpacking
  486.     // Normalize denormalized numbers.
  487.     // Insert assumed high-order bit for normalized numbers.
  488.     // Subtract exponent bias.
  489.     if ( binExp == 0 ){
  490.         if ( fractBits == 0 ){
  491.         // not a denorm, just a 0!
  492.         decExponent = 0;
  493.         digits = zero;
  494.         nDigits = 1;
  495.         return;
  496.         }
  497.         while ( (fractBits&singleFractHOB) == 0 ){
  498.         fractBits <<= 1;
  499.         binExp -= 1;
  500.         }
  501.         nSignificantBits = singleExpShift + binExp +1; // recall binExp is  - shift count.
  502.         binExp += 1;
  503.     } else {
  504.         fractBits |= singleFractHOB;
  505.         nSignificantBits = singleExpShift+1;
  506.     }
  507.     binExp -= singleExpBias;
  508.     // call the routine that actually does all the hard work.
  509.     dtoa( binExp, ((long)fractBits)<<(expShift-singleExpShift), nSignificantBits );
  510.     }
  511.  
  512.     private void
  513.     dtoa( int binExp, long fractBits, int nSignificantBits )
  514.     {
  515.     int    nFractBits; // number of significant bits of fractBits;
  516.     int    nTinyBits;  // number of these to the right of the point.
  517.     int    decExp;
  518.  
  519.     // Examine number. Determine if it is an easy case,
  520.     // which we can do pretty trivially using float/long conversion,
  521.     // or whether we must do real work.
  522.     nFractBits = countBits( fractBits );
  523.     nTinyBits = Math.max( 0, nFractBits - binExp - 1 );
  524.     if ( binExp <= maxSmallBinExp && binExp >= minSmallBinExp ){
  525.         // Look more closely at the number to decide if,
  526.         // with scaling by 10^nTinyBits, the result will fit in
  527.         // a long.
  528.         if ( (nTinyBits < long5pow.length) && ((nFractBits + n5bits[nTinyBits]) < 64 ) ){
  529.         /*
  530.          * We can do this:
  531.          * take the fraction bits, which are normalized.
  532.          * (a) nTinyBits == 0: Shift left or right appropriately
  533.          *     to align the binary point at the extreme right, i.e.
  534.          *     where a long int point is expected to be. The integer
  535.          *     result is easily converted to a string.
  536.          * (b) nTinyBits > 0: Shift right by expShift-nFractBits,
  537.          *     which effectively converts to long and scales by
  538.          *     2^nTinyBits. Then multiply by 5^nTinyBits to
  539.          *     complete the scaling. We know this won't overflow
  540.          *     because we just counted the number of bits necessary
  541.          *     in the result. The integer you get from this can
  542.          *     then be converted to a string pretty easily.
  543.          */
  544.         long halfULP;
  545.         if ( nTinyBits == 0 ) {
  546.             if ( binExp > nSignificantBits ){
  547.             halfULP = 1L << ( binExp-nSignificantBits-1);
  548.             } else {
  549.             halfULP = 0L;
  550.             }
  551.             if ( binExp >= expShift ){
  552.             fractBits <<= (binExp-expShift);
  553.             } else {
  554.             fractBits >>>= (expShift-binExp) ;
  555.             }
  556.             developLongDigits( 0, fractBits, halfULP );
  557.             return;
  558.         }
  559.         /*
  560.          * The following causes excess digits to be printed
  561.          * out in the single-float case. Our manipulation of
  562.          * halfULP here is apparently not correct. If we
  563.          * better understand how this works, perhaps we can
  564.          * use this special case again. But for the time being,
  565.          * we do not.
  566.          * else {
  567.          *     fractBits >>>= expShift+1-nFractBits;
  568.          *     fractBits *= long5pow[ nTinyBits ];
  569.          *     halfULP = long5pow[ nTinyBits ] >> (1+nSignificantBits-nFractBits);
  570.          *     developLongDigits( -nTinyBits, fractBits, halfULP );
  571.          *     return;
  572.          * }
  573.          */
  574.         }
  575.     }
  576.     /*
  577.      * This is the hard case. We are going to compute large positive
  578.      * integers B and S and integer decExp, s.t.
  579.      *    d = ( B / S ) * 10^decExp
  580.      *    1 <= B / S < 10
  581.      * Obvious choices are:
  582.      *    decExp = floor( log10(d) )
  583.      *     B      = d * 2^nTinyBits * 10^max( 0, -decExp )
  584.      *    S      = 10^max( 0, decExp) * 2^nTinyBits
  585.      * (noting that nTinyBits has already been forced to non-negative)
  586.      * I am also going to compute a large positive integer
  587.      *    M      = (1/2^nSignificantBits) * 2^nTinyBits * 10^max( 0, -decExp )
  588.      * i.e. M is (1/2) of the ULP of d, scaled like B.
  589.      * When we iterate through dividing B/S and picking off the
  590.      * quotient bits, we will know when to stop when the remainder
  591.      * is <= M.
  592.      *
  593.      * We keep track of powers of 2 and powers of 5.
  594.      */
  595.  
  596.     /*
  597.      * Estimate decimal exponent. (If it is small-ish,
  598.      * we could double-check.)
  599.      *
  600.      * First, scale the mantissa bits such that 1 <= d2 < 2.
  601.      * We are then going to estimate
  602.      *        log10(d2) ~=~  (d2-1.5)/1.5 + log(1.5)
  603.      * and so we can estimate
  604.      *      log10(d) ~=~ log10(d2) + binExp * log10(2)
  605.      * take the floor and call it decExp.
  606.      * FIXME -- use more precise constants here. It costs no more.
  607.      */
  608.     double d2 = Double.longBitsToDouble(
  609.         expOne | ( fractBits &~ fractHOB ) );
  610.     decExp = (int)Math.floor(
  611.         (d2-1.5D)*0.289529654D + 0.176091259 + (double)binExp * 0.301029995663981 );
  612.     int B2, B5; // powers of 2 and powers of 5, respectively, in B
  613.     int S2, S5; // powers of 2 and powers of 5, respectively, in S
  614.     int M2, M5; // powers of 2 and powers of 5, respectively, in M
  615.     int Bbits; // binary digits needed to represent B, approx.
  616.     int tenSbits; // binary digits needed to represent 10*S, approx.
  617.     FDBigInt Sval, Bval, Mval;
  618.  
  619.     B5 = Math.max( 0, -decExp );
  620.     B2 = B5 + nTinyBits + binExp;
  621.  
  622.     S5 = Math.max( 0, decExp );
  623.     S2 = S5 + nTinyBits;
  624.  
  625.     M5 = B5;
  626.     M2 = B2 - nSignificantBits;
  627.  
  628.     /*
  629.      * the long integer fractBits contains the (nFractBits) interesting
  630.      * bits from the mantissa of d ( hidden 1 added if necessary) followed
  631.      * by (expShift+1-nFractBits) zeros. In the interest of compactness,
  632.      * I will shift out those zeros before turning fractBits into a
  633.      * FDBigInt. The resulting whole number will be
  634.      *     d * 2^(nFractBits-1-binExp).
  635.      */
  636.     fractBits >>>= (expShift+1-nFractBits);
  637.     B2 -= nFractBits-1;
  638.     int common2factor = Math.min( B2, S2 );
  639.     B2 -= common2factor;
  640.     S2 -= common2factor;
  641.     M2 -= common2factor;
  642.  
  643.     /*
  644.      * HACK!! For exact powers of two, the next smallest number
  645.      * is only half as far away as we think (because the meaning of
  646.      * ULP changes at power-of-two bounds) for this reason, we
  647.      * hack M2. Hope this works.
  648.      */
  649.     if ( nFractBits == 1 )
  650.         M2 -= 1;
  651.  
  652.     if ( M2 < 0 ){
  653.         // oops.
  654.         // since we cannot scale M down far enough,
  655.         // we must scale the other values up.
  656.         B2 -= M2;
  657.         S2 -= M2;
  658.         M2 =  0;
  659.     }
  660.     /*
  661.      * Construct, Scale, iterate.
  662.      * Some day, we'll write a stopping test that takes
  663.      * account of the assymetry of the spacing of floating-point
  664.      * numbers below perfect powers of 2
  665.      * 26 Sept 96 is not that day.
  666.      * So we use a symmetric test.
  667.      */
  668.     char digits[] = this.digits = new char[18];
  669.     int  ndigit = 0;
  670.     boolean low, high;
  671.     long lowDigitDifference;
  672.     int  q;
  673.  
  674.     /*
  675.      * Detect the special cases where all the numbers we are about
  676.      * to compute will fit in int or long integers.
  677.      * In these cases, we will avoid doing FDBigInt arithmetic.
  678.      * We use the same algorithms, except that we "normalize"
  679.      * our FDBigInts before iterating. This is to make division easier,
  680.      * as it makes our fist guess (quotient of high-order words)
  681.      * more accurate!
  682.      *
  683.      * Some day, we'll write a stopping test that takes
  684.      * account of the assymetry of the spacing of floating-point
  685.      * numbers below perfect powers of 2
  686.      * 26 Sept 96 is not that day.
  687.      * So we use a symmetric test.
  688.      */
  689.     Bbits = nFractBits + B2 + (( B5 < n5bits.length )? n5bits[B5] : ( B5*3 ));
  690.     tenSbits = S2+1 + (( (S5+1) < n5bits.length )? n5bits[(S5+1)] : ( (S5+1)*3 ));
  691.     if ( Bbits < 64 && tenSbits < 64){
  692.         if ( Bbits < 32 && tenSbits < 32){
  693.         // wa-hoo! They're all ints!
  694.         int b = ((int)fractBits * small5pow[B5] ) << B2;
  695.         int s = small5pow[S5] << S2;
  696.         int m = small5pow[M5] << M2;
  697.         int tens = s * 10;
  698.         /*
  699.          * Unroll the first iteration. If our decExp estimate
  700.          * was too high, our first quotient will be zero. In this
  701.          * case, we discard it and decrement decExp.
  702.          */
  703.         ndigit = 0;
  704.         q = (int) ( b / s );
  705.         b = 10 * ( b % s );
  706.         m *= 10;
  707.         low  = (b <  m );
  708.         high = (b+m > tens );
  709.         if ( q >= 10 ){
  710.             // bummer, dude
  711.             throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
  712.         } else if ( (q == 0) && ! high ){
  713.             // oops. Usually ignore leading zero.
  714.             decExp--;
  715.         } else {
  716.             digits[ndigit++] = (char)('0' + q);
  717.         }
  718.         /*
  719.          * HACK! Java spec sez that we always have at least
  720.          * one digit after the . in either F- or E-form output.
  721.          * Thus we will need more than one digit if we're using
  722.          * E-form
  723.          */
  724.         if ( decExp <= -3 || decExp >= 8 ){
  725.             high = low = false;
  726.         }
  727.         while( ! low && ! high ){
  728.             q = (int) ( b / s );
  729.             b = 10 * ( b % s );
  730.             m *= 10;
  731.             if ( q >= 10 ){
  732.             // bummer, dude
  733.             throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
  734.             }
  735.             if ( m > 0L ){
  736.             low  = (b <  m );
  737.             high = (b+m > tens );
  738.             } else {
  739.             // hack -- m might overflow!
  740.             // in this case, it is certainly > b,
  741.             // which won't
  742.             // and b+m > tens, too, since that has overflowed
  743.             // either!
  744.             low = true;
  745.             high = true;
  746.             }
  747.             digits[ndigit++] = (char)('0' + q);
  748.         }
  749.         lowDigitDifference = (b<<1) - tens;
  750.         } else {
  751.         // still good! they're all longs!
  752.         long b = (fractBits * long5pow[B5] ) << B2;
  753.         long s = long5pow[S5] << S2;
  754.         long m = long5pow[M5] << M2;
  755.         long tens = s * 10L;
  756.         /*
  757.          * Unroll the first iteration. If our decExp estimate
  758.          * was too high, our first quotient will be zero. In this
  759.          * case, we discard it and decrement decExp.
  760.          */
  761.         ndigit = 0;
  762.         q = (int) ( b / s );
  763.         b = 10L * ( b % s );
  764.         m *= 10L;
  765.         low  = (b <  m );
  766.         high = (b+m > tens );
  767.         if ( q >= 10 ){
  768.             // bummer, dude
  769.             throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
  770.         } else if ( (q == 0) && ! high ){
  771.             // oops. Usually ignore leading zero.
  772.             decExp--;
  773.         } else {
  774.             digits[ndigit++] = (char)('0' + q);
  775.         }
  776.         /*
  777.          * HACK! Java spec sez that we always have at least
  778.          * one digit after the . in either F- or E-form output.
  779.          * Thus we will need more than one digit if we're using
  780.          * E-form
  781.          */
  782.         if ( decExp <= -3 || decExp >= 8 ){
  783.             high = low = false;
  784.         }
  785.         while( ! low && ! high ){
  786.             q = (int) ( b / s );
  787.             b = 10 * ( b % s );
  788.             m *= 10;
  789.             if ( q >= 10 ){
  790.             // bummer, dude
  791.             throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
  792.             }
  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.         if ( q >= 10 ){
  837.         // bummer, dude
  838.         throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
  839.         } else if ( (q == 0) && ! high ){
  840.         // oops. Usually ignore leading zero.
  841.         decExp--;
  842.         } else {
  843.         digits[ndigit++] = (char)('0' + q);
  844.         }
  845.         /*
  846.          * HACK! Java spec sez that we always have at least
  847.          * one digit after the . in either F- or E-form output.
  848.          * Thus we will need more than one digit if we're using
  849.          * E-form
  850.          */
  851.         if ( decExp <= -3 || decExp >= 8 ){
  852.         high = low = false;
  853.         }
  854.         while( ! low && ! high ){
  855.         q = Bval.quoRemIteration( Sval );
  856.         Mval = Mval.mult( 10 );
  857.         if ( q >= 10 ){
  858.             // bummer, dude
  859.             throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
  860.         }
  861.         low  = (Bval.cmp( Mval ) < 0);
  862.         high = (Bval.add( Mval ).cmp( tenSval ) > 0 );
  863.         digits[ndigit++] = (char)('0' + q);
  864.         }
  865.         if ( high && low ){
  866.         Bval.lshiftMe(1);
  867.         lowDigitDifference = Bval.cmp(tenSval);
  868.         } else
  869.         lowDigitDifference = 0L; // this here only for flow analysis!
  870.     }
  871.     this.decExponent = decExp+1;
  872.     this.digits = digits;
  873.     this.nDigits = ndigit;
  874.     /*
  875.      * Last digit gets rounded based on stopping condition.
  876.      */
  877.     if ( high ){
  878.         if ( low ){
  879.         if ( lowDigitDifference == 0L ){
  880.             // it's a tie!
  881.             // choose based on which digits we like.
  882.             if ( (digits[nDigits-1]&1) != 0 ) roundup();
  883.         } else if ( lowDigitDifference > 0 ){
  884.             roundup();
  885.         }
  886.         } else {
  887.         roundup();
  888.         }
  889.     }
  890.     }
  891.  
  892.     public String
  893.     toString(){
  894.     // most brain-dead version
  895.     StringBuffer result = new StringBuffer( nDigits+8 );
  896.     if ( isNegative ){ result.append( '-' ); }
  897.     if ( isExceptional ){
  898.         result.append( digits, 0, nDigits );
  899.     } else {
  900.         result.append( "0.");
  901.         result.append( digits, 0, nDigits );
  902.         result.append('e');
  903.         result.append( decExponent );
  904.     }
  905.     return new String(result);
  906.     }
  907.  
  908.     public String
  909.     toJavaFormatString(){
  910.     char result[] = new char[ nDigits + 10 ];
  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 new String(result, 0, i);
  979.     }
  980.  
  981.     public static FloatingDecimal
  982.     readJavaFormatString( String in ) throws NumberFormatException {
  983.     boolean isNegative = false;
  984.     boolean signSeen   = false;
  985.     int     decExp;
  986.     char    c;
  987.  
  988.     parseNumber:
  989.     try{
  990.         in = in.trim(); // don't fool around with white space.
  991.                // throws NullPointerException if null
  992.         int    l = in.length();
  993.         if ( l == 0 ) throw new NumberFormatException("empty String");
  994.         int    i = 0;
  995.         switch ( c = in.charAt( i ) ){
  996.         case '-':
  997.         isNegative = true;
  998.         //FALLTHROUGH
  999.         case '+':
  1000.         i++;
  1001.         signSeen = true;
  1002.         }
  1003.         // Would handle NaN and Infinity here, but it isn't
  1004.         // part of the spec!
  1005.         //
  1006.         char[] digits = new char[ l ];
  1007.         int    nDigits= 0;
  1008.         boolean decSeen = false;
  1009.         int    decPt = 0;
  1010.         int    nLeadZero = 0;
  1011.         int    nTrailZero= 0;
  1012.     digitLoop:
  1013.         while ( i < l ){
  1014.         switch ( c = in.charAt( i ) ){
  1015.         case '0':
  1016.             if ( nDigits > 0 ){
  1017.             nTrailZero += 1;
  1018.             } else {
  1019.             nLeadZero += 1;
  1020.             }
  1021.             break; // out of switch.
  1022.         case '1':
  1023.         case '2':
  1024.         case '3':
  1025.         case '4':
  1026.         case '5':
  1027.         case '6':
  1028.         case '7':
  1029.         case '8':
  1030.         case '9':
  1031.             while ( nTrailZero > 0 ){
  1032.             digits[nDigits++] = '0';
  1033.             nTrailZero -= 1;
  1034.             }
  1035.             digits[nDigits++] = c;
  1036.             break; // out of switch.
  1037.         case '.':
  1038.             if ( decSeen ){
  1039.             // already saw one ., this is the 2nd.
  1040.             throw new NumberFormatException("multiple points");
  1041.             }
  1042.             decPt = i;
  1043.             if ( signSeen ){
  1044.             decPt -= 1;
  1045.             }
  1046.             decSeen = true;
  1047.             break; // out of switch.
  1048.         default:
  1049.             break digitLoop;
  1050.         }
  1051.         i++;
  1052.         }
  1053.         /*
  1054.          * At this point, we've scanned all the digits and decimal
  1055.          * point we're going to see. Trim off leading and trailing
  1056.          * zeros, which will just confuse us later, and adjust
  1057.          * our initial decimal exponent accordingly.
  1058.          * To review:
  1059.          * we have seen i total characters.
  1060.          * nLeadZero of them were zeros before any other digits.
  1061.          * nTrailZero of them were zeros after any other digits.
  1062.          * if ( decSeen ), then a . was seen after decPt characters
  1063.          * ( including leading zeros which have been discarded )
  1064.          * nDigits characters were neither lead nor trailing
  1065.          * zeros, nor point
  1066.          */
  1067.         /*
  1068.          * special hack: if we saw no non-zero digits, then the
  1069.          * answer is zero!
  1070.          * Unfortunately, we feel honor-bound to keep parsing!
  1071.          */
  1072.         if ( nDigits == 0 ){
  1073.         digits = zero;
  1074.         nDigits = 1;
  1075.         if ( nLeadZero == 0 ){
  1076.             // we saw NO DIGITS AT ALL,
  1077.             // not even a crummy 0!
  1078.             // this is not allowed.
  1079.             break parseNumber; // go throw exception
  1080.         }
  1081.  
  1082.         }
  1083.  
  1084.         /* Our initial exponent is decPt, adjusted by the number of
  1085.          * discarded zeros. Or, if there was no decPt,
  1086.          * then its just nDigits adjusted by discarded trailing zeros.
  1087.          */
  1088.         if ( decSeen ){
  1089.         decExp = decPt - nLeadZero;
  1090.         } else {
  1091.         decExp = nDigits+nTrailZero;
  1092.         }
  1093.  
  1094.         /*
  1095.          * Look for 'e' or 'E' and an optionally signed integer.
  1096.          */
  1097.         if ( (i < l) &&  ((c = in.charAt(i) )=='e') || (c == 'E') ){
  1098.         int expSign = 1;
  1099.         int expVal  = 0;
  1100.         int reallyBig = Integer.MAX_VALUE / 10;
  1101.         boolean expOverflow = false;
  1102.         switch( in.charAt(++i) ){
  1103.         case '-':
  1104.             expSign = -1;
  1105.             //FALLTHROUGH
  1106.         case '+':
  1107.             i++;
  1108.         }
  1109.         int expAt = i;
  1110.         expLoop:
  1111.         while ( i < l  ){
  1112.             if ( expVal >= reallyBig ){
  1113.             // the next character will cause integer
  1114.             // overflow.
  1115.             expOverflow = true;
  1116.             }
  1117.             switch ( c = in.charAt(i++) ){
  1118.             case '0':
  1119.             case '1':
  1120.             case '2':
  1121.             case '3':
  1122.             case '4':
  1123.             case '5':
  1124.             case '6':
  1125.             case '7':
  1126.             case '8':
  1127.             case '9':
  1128.             expVal = expVal*10 + ( (int)c - (int)'0' );
  1129.             continue;
  1130.             default:
  1131.             i--;           // back up.
  1132.             break expLoop; // stop parsing exponent.
  1133.             }
  1134.         }
  1135.         int expLimit = bigDecimalExponent+nDigits+nTrailZero;
  1136.         if ( expOverflow || ( expVal > expLimit ) ){
  1137.             //
  1138.             // The intent here is to end up with
  1139.             // infinity or zero, as appropriate.
  1140.             // The reason for yielding such a small decExponent,
  1141.             // rather than something intuitive such as
  1142.             // expSign*Integer.MAX_VALUE, is that this value
  1143.             // is subject to further manipulation in
  1144.             // doubleValue() and floatValue(), and I don't want
  1145.             // it to be able to cause overflow there!
  1146.             // (The only way we can get into trouble here is for
  1147.             // really outrageous nDigits+nTrailZero, such as 2 billion. )
  1148.             //
  1149.             decExp = expSign*expLimit;
  1150.         } else {
  1151.             // this should not overflow, since we tested
  1152.             // for expVal > (MAX+N), where N >= abs(decExp)
  1153.             decExp = decExp + expSign*expVal;
  1154.         }
  1155.  
  1156.         // if we saw something not a digit ( or end of string )
  1157.         // after the [Ee][+-], without seeing any digits at all
  1158.         // this is certainly an error. If we saw some digits,
  1159.         // but then some trailing garbage, that might be ok.
  1160.         // so we just fall through in that case.
  1161.         // HUMBUG
  1162.                 if ( i == expAt ) 
  1163.             break parseNumber; // certainly bad
  1164.         }
  1165.         /*
  1166.          * We parsed everything we could.
  1167.          * If there are leftovers, then this is not good input!
  1168.          */
  1169.         if ( i < l &&
  1170.                 ((i != l - 1) ||
  1171.                 (in.charAt(i) != 'f' &&
  1172.                  in.charAt(i) != 'F' &&
  1173.                  in.charAt(i) != 'd' &&
  1174.                  in.charAt(i) != 'D'))) {
  1175.                 break parseNumber; // go throw exception
  1176.             }
  1177.  
  1178.         return new FloatingDecimal( isNegative, decExp, digits, nDigits,  false );
  1179.     } catch ( StringIndexOutOfBoundsException e ){ }
  1180.     throw new NumberFormatException( in );
  1181.     }
  1182.  
  1183.     /*
  1184.      * Take a FloatingDecimal, which we presumably just scanned in,
  1185.      * and find out what its value is, as a double.
  1186.      *
  1187.      * AS A SIDE EFFECT, SET roundDir TO INDICATE PREFERRED
  1188.      * ROUNDING DIRECTION in case the result is really destined
  1189.      * for a single-precision float.
  1190.      */
  1191.  
  1192.     public double
  1193.     doubleValue(){
  1194.     int    kDigits = Math.min( nDigits, maxDecimalDigits+1 );
  1195.     long    lValue;
  1196.     double    dValue;
  1197.     double  rValue, tValue;
  1198.  
  1199.     roundDir = 0;
  1200.  
  1201.     /*
  1202.      * convert the lead kDigits to a long integer.
  1203.      */
  1204.     // (special performance hack: start to do it using int)
  1205.     int iValue = (int)digits[0]-(int)'0';
  1206.     int iDigits = Math.min( kDigits, intDecimalDigits );
  1207.     for ( int i=1; i < iDigits; i++ ){
  1208.         iValue = iValue*10 + (int)digits[i]-(int)'0';
  1209.     }
  1210.     lValue = (long)iValue;
  1211.     for ( int i=iDigits; i < kDigits; i++ ){
  1212.         lValue = lValue*10L + (long)((int)digits[i]-(int)'0');
  1213.     }
  1214.     dValue = (double)lValue;
  1215.     int exp = decExponent-kDigits;
  1216.     /*
  1217.      * lValue now contains a long integer with the value of
  1218.      * the first kDigits digits of the number.
  1219.      * dValue contains the (double) of the same.
  1220.      */
  1221.  
  1222.     if ( nDigits <= maxDecimalDigits ){
  1223.         /*
  1224.          * possibly an easy case.
  1225.          * We know that the digits can be represented
  1226.          * exactly. And if the exponent isn't too outrageous,
  1227.          * the whole thing can be done with one operation,
  1228.          * thus one rounding error.
  1229.          */
  1230.         if ( exp == 0 ) return (isNegative)? -dValue : dValue; // small floating integer
  1231.         else if ( exp >= 0 ){
  1232.         if ( exp <= maxSmallTen ){
  1233.             /*
  1234.              * Can get the answer with one operation,
  1235.              * thus one roundoff.
  1236.              */
  1237.             rValue = dValue * small10pow[exp];
  1238.             if ( mustSetRoundDir ){
  1239.             tValue = rValue / small10pow[exp];
  1240.             roundDir = ( tValue ==  dValue ) ? 0
  1241.                       :( tValue < dValue ) ? 1
  1242.                       : -1;
  1243.             }
  1244.             return (isNegative)? -rValue : rValue;
  1245.         }
  1246.         int slop = maxDecimalDigits - kDigits;
  1247.         if ( exp <= maxSmallTen+slop ){
  1248.             /*
  1249.              * We can multiply dValue by 10^(slop)
  1250.              * and it is still "small" and exact.
  1251.              * Then we can multiply by 10^(exp-slop)
  1252.              * with one rounding.
  1253.              */
  1254.             dValue *= small10pow[slop];
  1255.             rValue = dValue * small10pow[exp-slop];
  1256.  
  1257.             if ( mustSetRoundDir ){
  1258.             tValue = rValue / small10pow[exp-slop];
  1259.             roundDir = ( tValue ==  dValue ) ? 0
  1260.                   :( tValue < dValue ) ? 1
  1261.                   : -1;
  1262.             }
  1263.             return (isNegative)? -rValue : rValue;
  1264.         }
  1265.         /*
  1266.          * Else we have a hard case with a positive exp.
  1267.          */
  1268.         } else {
  1269.         if ( exp >= -maxSmallTen ){
  1270.             /*
  1271.              * Can get the answer in one division.
  1272.              */
  1273.             rValue = dValue / small10pow[-exp];
  1274.             tValue = rValue * small10pow[-exp];
  1275.             if ( mustSetRoundDir ){
  1276.             roundDir = ( tValue ==  dValue ) ? 0
  1277.                   :( tValue < dValue ) ? 1
  1278.                   : -1;
  1279.             }
  1280.             return (isNegative)? -rValue : rValue;
  1281.         }
  1282.         /*
  1283.          * Else we have a hard case with a negative exp.
  1284.          */
  1285.         }
  1286.     }
  1287.  
  1288.     /*
  1289.      * Harder cases:
  1290.      * The sum of digits plus exponent is greater than
  1291.      * what we think we can do with one error.
  1292.      *
  1293.      * Start by approximating the right answer by,
  1294.      * naively, scaling by powers of 10.
  1295.      */
  1296.     if ( exp > 0 ){
  1297.         if ( decExponent > maxDecimalExponent+1 ){
  1298.         /*
  1299.          * Lets face it. This is going to be
  1300.          * Infinity. Cut to the chase.
  1301.          */
  1302.         return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
  1303.         }
  1304.         if ( (exp&15) != 0 ){
  1305.         dValue *= small10pow[exp&15];
  1306.         }
  1307.         if ( (exp>>=4) != 0 ){
  1308.         int j;
  1309.         for( j = 0; exp > 1; j++, exp>>=1 ){
  1310.             if ( (exp&1)!=0)
  1311.             dValue *= big10pow[j];
  1312.         }
  1313.         /*
  1314.          * The reason for the weird exp > 1 condition
  1315.          * in the above loop was so that the last multiply
  1316.          * would get unrolled. We handle it here.
  1317.          * It could overflow.
  1318.          */
  1319.         double t = dValue * big10pow[j];
  1320.         if ( Double.isInfinite( t ) ){
  1321.             /*
  1322.              * It did overflow.
  1323.              * Look more closely at the result.
  1324.              * If the exponent is just one too large,
  1325.              * then use the maximum finite as our estimate
  1326.              * value. Else call the result infinity
  1327.              * and punt it.
  1328.              * ( I presume this could happen because
  1329.              * rounding forces the result here to be
  1330.              * an ULP or two larger than
  1331.              * Double.MAX_VALUE ).
  1332.              */
  1333.             t = dValue / 2.0;
  1334.             t *= big10pow[j];
  1335.             if ( Double.isInfinite( t ) ){
  1336.             return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
  1337.             }
  1338.             t = Double.MAX_VALUE;
  1339.         }
  1340.         dValue = t;
  1341.         }
  1342.     } else if ( exp < 0 ){
  1343.         exp = -exp;
  1344.         if ( decExponent < minDecimalExponent-1 ){
  1345.         /*
  1346.          * Lets face it. This is going to be
  1347.          * zero. Cut to the chase.
  1348.          */
  1349.         return (isNegative)? -0.0 : 0.0;
  1350.         }
  1351.         if ( (exp&15) != 0 ){
  1352.         dValue /= small10pow[exp&15];
  1353.         }
  1354.         if ( (exp>>=4) != 0 ){
  1355.         int j;
  1356.         for( j = 0; exp > 1; j++, exp>>=1 ){
  1357.             if ( (exp&1)!=0)
  1358.             dValue *= tiny10pow[j];
  1359.         }
  1360.         /*
  1361.          * The reason for the weird exp > 1 condition
  1362.          * in the above loop was so that the last multiply
  1363.          * would get unrolled. We handle it here.
  1364.          * It could underflow.
  1365.          */
  1366.         double t = dValue * tiny10pow[j];
  1367.         if ( t == 0.0 ){
  1368.             /*
  1369.              * It did underflow.
  1370.              * Look more closely at the result.
  1371.              * If the exponent is just one too small,
  1372.              * then use the minimum finite as our estimate
  1373.              * value. Else call the result 0.0
  1374.              * and punt it.
  1375.              * ( I presume this could happen because
  1376.              * rounding forces the result here to be
  1377.              * an ULP or two less than
  1378.              * Double.MIN_VALUE ).
  1379.              */
  1380.             t = dValue * 2.0;
  1381.             t *= tiny10pow[j];
  1382.             if ( t == 0.0 ){
  1383.             return (isNegative)? -0.0 : 0.0;
  1384.             }
  1385.             t = Double.MIN_VALUE;
  1386.         }
  1387.         dValue = t;
  1388.         }
  1389.     }
  1390.  
  1391.     /*
  1392.      * dValue is now approximately the result.
  1393.      * The hard part is adjusting it, by comparison
  1394.      * with FDBigInt arithmetic.
  1395.      * Formulate the EXACT big-number result as
  1396.      * bigD0 * 10^exp
  1397.      */
  1398.     FDBigInt bigD0 = new FDBigInt( lValue, digits, kDigits, nDigits );
  1399.     exp   = decExponent - nDigits;
  1400.  
  1401.     correctionLoop:
  1402.     while(true){
  1403.         /* AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES
  1404.          * bigIntExp and bigIntNBits
  1405.          */
  1406.         FDBigInt bigB = doubleToBigInt( dValue );
  1407.  
  1408.         /*
  1409.          * Scale bigD, bigB appropriately for
  1410.          * big-integer operations.
  1411.          * Naively, we multipy by powers of ten
  1412.          * and powers of two. What we actually do
  1413.          * is keep track of the powers of 5 and
  1414.          * powers of 2 we would use, then factor out
  1415.          * common divisors before doing the work.
  1416.          */
  1417.         int B2, B5; // powers of 2, 5 in bigB
  1418.         int    D2, D5;    // powers of 2, 5 in bigD
  1419.         int Ulp2;   // powers of 2 in halfUlp.
  1420.         if ( exp >= 0 ){
  1421.         B2 = B5 = 0;
  1422.         D2 = D5 = exp;
  1423.         } else {
  1424.         B2 = B5 = -exp;
  1425.         D2 = D5 = 0;
  1426.         }
  1427.         if ( bigIntExp >= 0 ){
  1428.         B2 += bigIntExp;
  1429.         } else {
  1430.         D2 -= bigIntExp;
  1431.         }
  1432.         Ulp2 = B2;
  1433.         // shift bigB and bigD left by a number s. t.
  1434.         // halfUlp is still an integer.
  1435.         int hulpbias;
  1436.         if ( bigIntExp+bigIntNBits <= -expBias+1 ){
  1437.         // This is going to be a denormalized number
  1438.         // (if not actually zero).
  1439.         // half an ULP is at 2^-(expBias+expShift+1)
  1440.         hulpbias = bigIntExp+ expBias + expShift;
  1441.         } else {
  1442.         hulpbias = expShift + 2 - bigIntNBits;
  1443.         }
  1444.         B2 += hulpbias;
  1445.         D2 += hulpbias;
  1446.         // if there are common factors of 2, we might just as well
  1447.         // factor them out, as they add nothing useful.
  1448.         int common2 = Math.min( B2, Math.min( D2, Ulp2 ) );
  1449.         B2 -= common2;
  1450.         D2 -= common2;
  1451.         Ulp2 -= common2;
  1452.         // do multiplications by powers of 5 and 2
  1453.         bigB = multPow52( bigB, B5, B2 );
  1454.         FDBigInt bigD = multPow52( new FDBigInt( bigD0 ), D5, D2 );
  1455.         //
  1456.         // to recap:
  1457.         // bigB is the scaled-big-int version of our floating-point
  1458.         // candidate.
  1459.         // bigD is the scaled-big-int version of the exact value
  1460.         // as we understand it.
  1461.         // halfUlp is 1/2 an ulp of bigB, except for special cases
  1462.         // of exact powers of 2
  1463.         //
  1464.         // the plan is to compare bigB with bigD, and if the difference
  1465.         // is less than halfUlp, then we're satisfied. Otherwise,
  1466.         // use the ratio of difference to halfUlp to calculate a fudge
  1467.         // factor to add to the floating value, then go 'round again.
  1468.         //
  1469.         FDBigInt diff;
  1470.         int cmpResult;
  1471.         boolean overvalue;
  1472.         if ( (cmpResult = bigB.cmp( bigD ) ) > 0 ){
  1473.         overvalue = true; // our candidate is too big.
  1474.         diff = bigB.sub( bigD );
  1475.         if ( (bigIntNBits == 1) && (bigIntExp > -expBias) ){
  1476.             // candidate is a normalized exact power of 2 and
  1477.             // is too big. We will be subtracting.
  1478.             // For our purposes, ulp is the ulp of the
  1479.             // next smaller range.
  1480.             Ulp2 -= 1;
  1481.             if ( Ulp2 < 0 ){
  1482.             // rats. Cannot de-scale ulp this far.
  1483.             // must scale diff in other direction.
  1484.             Ulp2 = 0;
  1485.             diff.lshiftMe( 1 );
  1486.             }
  1487.         }
  1488.         } else if ( cmpResult < 0 ){
  1489.         overvalue = false; // our candidate is too small.
  1490.         diff = bigD.sub( bigB );
  1491.         } else {
  1492.         // the candidate is exactly right!
  1493.         // this happens with surprising fequency
  1494.         break correctionLoop;
  1495.         }
  1496.         FDBigInt halfUlp = constructPow52( B5, Ulp2 );
  1497.         if ( (cmpResult = diff.cmp( halfUlp ) ) < 0 ){
  1498.         // difference is small.
  1499.         // this is close enough
  1500.         roundDir = overvalue ? -1 : 1;
  1501.         break correctionLoop;
  1502.         } else if ( cmpResult == 0 ){
  1503.         // difference is exactly half an ULP
  1504.         // round to some other value maybe, then finish
  1505.         dValue += 0.5*ulp( dValue, overvalue );
  1506.         // should check for bigIntNBits == 1 here??
  1507.         roundDir = overvalue ? -1 : 1;
  1508.         break correctionLoop;
  1509.         } else {
  1510.         // difference is non-trivial.
  1511.         // could scale addend by ratio of difference to
  1512.         // halfUlp here, if we bothered to compute that difference.
  1513.         // Most of the time ( I hope ) it is about 1 anyway.
  1514.         dValue += ulp( dValue, overvalue );
  1515.         if ( dValue == 0.0 || dValue == Double.POSITIVE_INFINITY )
  1516.              break correctionLoop; // oops. Fell off end of range.
  1517.         continue; // try again.
  1518.         }
  1519.  
  1520.     }
  1521.     return (isNegative)? -dValue : dValue;
  1522.     }
  1523.  
  1524.     /*
  1525.      * Take a FloatingDecimal, which we presumably just scanned in,
  1526.      * and find out what its value is, as a float.
  1527.      * This is distinct from doubleValue() to avoid the extremely
  1528.      * unlikely case of a double rounding error, wherein the converstion
  1529.      * to double has one rounding error, and the conversion of that double
  1530.      * to a float has another rounding error, IN THE WRONG DIRECTION,
  1531.      * ( because of the preference to a zero low-order bit ).
  1532.      */
  1533.  
  1534.     public float
  1535.     floatValue(){
  1536.     int    kDigits = Math.min( nDigits, singleMaxDecimalDigits+1 );
  1537.     int    iValue;
  1538.     float    fValue;
  1539.  
  1540.     /*
  1541.      * convert the lead kDigits to an integer.
  1542.      */
  1543.     iValue = (int)digits[0]-(int)'0';
  1544.     for ( int i=1; i < kDigits; i++ ){
  1545.         iValue = iValue*10 + (int)digits[i]-(int)'0';
  1546.     }
  1547.     fValue = (float)iValue;
  1548.     int exp = decExponent-kDigits;
  1549.     /*
  1550.      * iValue now contains an integer with the value of
  1551.      * the first kDigits digits of the number.
  1552.      * fValue contains the (float) of the same.
  1553.      */
  1554.  
  1555.     if ( nDigits <= singleMaxDecimalDigits ){
  1556.         /*
  1557.          * possibly an easy case.
  1558.          * We know that the digits can be represented
  1559.          * exactly. And if the exponent isn't too outrageous,
  1560.          * the whole thing can be done with one operation,
  1561.          * thus one rounding error.
  1562.          */
  1563.         if ( exp == 0 ) return (isNegative)? -fValue : fValue; // small floating integer
  1564.         else if ( exp >= 0 ){
  1565.         if ( exp <= singleMaxSmallTen ){
  1566.             /*
  1567.              * Can get the answer with one operation,
  1568.              * thus one roundoff.
  1569.              */
  1570.             fValue *= singleSmall10pow[exp];
  1571.             return (isNegative)? -fValue : fValue;
  1572.         }
  1573.         int slop = singleMaxDecimalDigits - kDigits;
  1574.         if ( exp <= singleMaxSmallTen+slop ){
  1575.             /*
  1576.              * We can multiply dValue by 10^(slop)
  1577.              * and it is still "small" and exact.
  1578.              * Then we can multiply by 10^(exp-slop)
  1579.              * with one rounding.
  1580.              */
  1581.             fValue *= singleSmall10pow[slop];
  1582.             fValue *= singleSmall10pow[exp-slop];
  1583.             return (isNegative)? -fValue : fValue;
  1584.         }
  1585.         /*
  1586.          * Else we have a hard case with a positive exp.
  1587.          */
  1588.         } else {
  1589.         if ( exp >= -singleMaxSmallTen ){
  1590.             /*
  1591.              * Can get the answer in one division.
  1592.              */
  1593.             fValue /= singleSmall10pow[-exp];
  1594.             return (isNegative)? -fValue : fValue;
  1595.         }
  1596.         /*
  1597.          * Else we have a hard case with a negative exp.
  1598.          */
  1599.         }
  1600.     } else if ( (decExponent >= nDigits) && (nDigits+decExponent <= maxDecimalDigits) ){
  1601.         /*
  1602.          * In double-precision, this is an exact floating integer.
  1603.          * So we can compute to double, then shorten to float
  1604.          * with one round, and get the right answer.
  1605.          *
  1606.          * First, finish accumulating digits.
  1607.          * Then convert that integer to a double, multiply
  1608.          * by the appropriate power of ten, and convert to float.
  1609.          */
  1610.         long lValue = (long)iValue;
  1611.         for ( int i=kDigits; i < nDigits; i++ ){
  1612.         lValue = lValue*10L + (long)((int)digits[i]-(int)'0');
  1613.         }
  1614.         double dValue = (double)lValue;
  1615.         exp = decExponent-nDigits;
  1616.         dValue *= small10pow[exp];
  1617.         fValue = (float)dValue;
  1618.         return (isNegative)? -fValue : fValue;
  1619.  
  1620.     }
  1621.     /*
  1622.      * Harder cases:
  1623.      * The sum of digits plus exponent is greater than
  1624.      * what we think we can do with one error.
  1625.      *
  1626.      * Start by weeding out obviously out-of-range
  1627.      * results, then convert to double and go to
  1628.      * common hard-case code.
  1629.      */
  1630.     if ( decExponent > singleMaxDecimalExponent+1 ){
  1631.         /*
  1632.          * Lets face it. This is going to be
  1633.          * Infinity. Cut to the chase.
  1634.          */
  1635.         return (isNegative)? Float.NEGATIVE_INFINITY : Float.POSITIVE_INFINITY;
  1636.     } else if ( decExponent < singleMinDecimalExponent-1 ){
  1637.         /*
  1638.          * Lets face it. This is going to be
  1639.          * zero. Cut to the chase.
  1640.          */
  1641.         return (isNegative)? -0.0f : 0.0f;
  1642.     }
  1643.  
  1644.     /*
  1645.      * Here, we do 'way too much work, but throwing away
  1646.      * our partial results, and going and doing the whole
  1647.      * thing as double, then throwing away half the bits that computes
  1648.      * when we convert back to float.
  1649.      *
  1650.      * The alternative is to reproduce the whole multiple-precision
  1651.      * algorythm for float precision, or to try to parameterize it
  1652.      * for common usage. The former will take about 400 lines of code,
  1653.      * and the latter I tried without success. Thus the semi-hack
  1654.      * answer here.
  1655.      */
  1656.     mustSetRoundDir = true;
  1657.     double dValue = doubleValue();
  1658.     return stickyRound( dValue );
  1659.     }
  1660.  
  1661.  
  1662.     /*
  1663.      * All the positive powers of 10 that can be
  1664.      * represented exactly in double/float.
  1665.      */
  1666.     private static final double small10pow[] = {
  1667.     1.0e0,
  1668.     1.0e1, 1.0e2, 1.0e3, 1.0e4, 1.0e5,
  1669.     1.0e6, 1.0e7, 1.0e8, 1.0e9, 1.0e10,
  1670.     1.0e11, 1.0e12, 1.0e13, 1.0e14, 1.0e15,
  1671.     1.0e16, 1.0e17, 1.0e18, 1.0e19, 1.0e20,
  1672.     1.0e21, 1.0e22
  1673.     };
  1674.  
  1675.     private static final float singleSmall10pow[] = {
  1676.     1.0e0f,
  1677.     1.0e1f, 1.0e2f, 1.0e3f, 1.0e4f, 1.0e5f,
  1678.     1.0e6f, 1.0e7f, 1.0e8f, 1.0e9f, 1.0e10f
  1679.     };
  1680.  
  1681.     private static final double big10pow[] = {
  1682.     1e16, 1e32, 1e64, 1e128, 1e256 };
  1683.     private static final double tiny10pow[] = {
  1684.     1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
  1685.  
  1686.     private static final int maxSmallTen = small10pow.length-1;
  1687.     private static final int singleMaxSmallTen = singleSmall10pow.length-1;
  1688.  
  1689.     private static final int small5pow[] = {
  1690.     1,
  1691.     5,
  1692.     5*5,
  1693.     5*5*5,
  1694.     5*5*5*5,
  1695.     5*5*5*5*5,
  1696.     5*5*5*5*5*5,
  1697.     5*5*5*5*5*5*5,
  1698.     5*5*5*5*5*5*5*5,
  1699.     5*5*5*5*5*5*5*5*5,
  1700.     5*5*5*5*5*5*5*5*5*5,
  1701.     5*5*5*5*5*5*5*5*5*5*5,
  1702.     5*5*5*5*5*5*5*5*5*5*5*5,
  1703.     5*5*5*5*5*5*5*5*5*5*5*5*5
  1704.     };
  1705.  
  1706.  
  1707.     private static final long long5pow[] = {
  1708.     1L,
  1709.     5L,
  1710.     5L*5,
  1711.     5L*5*5,
  1712.     5L*5*5*5,
  1713.     5L*5*5*5*5,
  1714.     5L*5*5*5*5*5,
  1715.     5L*5*5*5*5*5*5,
  1716.     5L*5*5*5*5*5*5*5,
  1717.     5L*5*5*5*5*5*5*5*5,
  1718.     5L*5*5*5*5*5*5*5*5*5,
  1719.     5L*5*5*5*5*5*5*5*5*5*5,
  1720.     5L*5*5*5*5*5*5*5*5*5*5*5,
  1721.     5L*5*5*5*5*5*5*5*5*5*5*5*5,
  1722.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5,
  1723.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  1724.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  1725.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  1726.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  1727.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  1728.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  1729.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  1730.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  1731.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  1732.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  1733.     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,
  1734.     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,
  1735.     };
  1736.  
  1737.     // approximately ceil( log2( long5pow[i] ) )
  1738.     private static final int n5bits[] = {
  1739.     0,
  1740.     3,
  1741.     5,
  1742.     7,
  1743.     10,
  1744.     12,
  1745.     14,
  1746.     17,
  1747.     19,
  1748.     21,
  1749.     24,
  1750.     26,
  1751.     28,
  1752.     31,
  1753.     33,
  1754.     35,
  1755.     38,
  1756.     40,
  1757.     42,
  1758.     45,
  1759.     47,
  1760.     49,
  1761.     52,
  1762.     54,
  1763.     56,
  1764.     59,
  1765.     61,
  1766.     };
  1767.  
  1768.     private static final char infinity[] = { 'I', 'n', 'f', 'i', 'n', 'i', 't', 'y' };
  1769.     private static final char notANumber[] = { 'N', 'a', 'N' };
  1770.     private static final char zero[] = { '0', '0', '0', '0', '0', '0', '0', '0' };
  1771.  
  1772. }
  1773.  
  1774. /*
  1775.  * A really, really simple bigint package
  1776.  * tailored to the needs of floating base conversion.
  1777.  */
  1778. class FDBigInt {
  1779.     int    nWords; // number of words used
  1780.     int data[]; // value: data[0] is least significant
  1781.  
  1782.  
  1783.     public FDBigInt( int v ){
  1784.     nWords = 1;
  1785.     data = new int[1];
  1786.     data[0] = v;
  1787.     }
  1788.  
  1789.     public FDBigInt( long v ){
  1790.     data = new int[2];
  1791.     data[0] = (int)v;
  1792.     data[1] = (int)(v>>>32);
  1793.     nWords = (data[1]==0) ? 1 : 2;
  1794.     }
  1795.  
  1796.     public FDBigInt( FDBigInt other ){
  1797.     data = new int[nWords = other.nWords];
  1798.     System.arraycopy( other.data, 0, data, 0, nWords );
  1799.     }
  1800.  
  1801.     private FDBigInt( int [] d, int n ){
  1802.     data = d;
  1803.     nWords = n;
  1804.     }
  1805.  
  1806.     public FDBigInt( long seed, char digit[], int nd0, int nd ){
  1807.     int n= (nd+8)/9;    // estimate size needed.
  1808.     if ( n < 2 ) n = 2;
  1809.     data = new int[n];    // allocate enough space
  1810.     data[0] = (int)seed;    // starting value
  1811.     data[1] = (int)(seed>>>32);
  1812.     nWords = (data[1]==0) ? 1 : 2;
  1813.     int i = nd0;
  1814.     int limit = nd-5;    // slurp digits 5 at a time.
  1815.     int v;
  1816.     while ( i < limit ){
  1817.         int ilim = i+5;
  1818.         v = (int)digit[i++]-(int)'0';
  1819.         while( i <ilim ){
  1820.         v = 10*v + (int)digit[i++]-(int)'0';
  1821.         }
  1822.         multaddMe( 100000, v); // ... where 100000 is 10^5.
  1823.     }
  1824.     int factor = 1;
  1825.     v = 0;
  1826.     while ( i < nd ){
  1827.         v = 10*v + (int)digit[i++]-(int)'0';
  1828.         factor *= 10;
  1829.     }
  1830.     if ( factor != 1 ){
  1831.         multaddMe( factor, v );
  1832.     }
  1833.     }
  1834.  
  1835.     /*
  1836.      * Left shift by c bits.
  1837.      * Shifts this in place.
  1838.      */
  1839.     public void
  1840.     lshiftMe( int c )throws IllegalArgumentException {
  1841.     if ( c <= 0 ){
  1842.         if ( c == 0 )
  1843.         return; // silly.
  1844.         else
  1845.         throw new IllegalArgumentException("negative shift count");
  1846.     }
  1847.     int wordcount = c>>5;
  1848.     int bitcount  = c & 0x1f;
  1849.     int anticount = 32-bitcount;
  1850.     int t[] = data;
  1851.     int s[] = data;
  1852.     if ( nWords+wordcount+1 > t.length ){
  1853.         // reallocate.
  1854.         t = new int[ nWords+wordcount+1 ];
  1855.     }
  1856.     int target = nWords+wordcount;
  1857.     int src    = nWords-1;
  1858.     if ( bitcount == 0 ){
  1859.         // special hack, since an anticount of 32 won't go!
  1860.         System.arraycopy( s, 0, t, wordcount, nWords );
  1861.         target = wordcount-1;
  1862.     } else {
  1863.         t[target--] = s[src]>>>anticount;
  1864.         while ( src >= 1 ){
  1865.         t[target--] = (s[src]<<bitcount) | (s[--src]>>>anticount);
  1866.         }
  1867.         t[target--] = s[src]<<bitcount;
  1868.     }
  1869.     while( target >= 0 ){
  1870.         t[target--] = 0;
  1871.     }
  1872.     data = t;
  1873.     nWords += wordcount + 1;
  1874.     // may have constructed high-order word of 0.
  1875.     // if so, trim it
  1876.     while ( nWords > 1 && data[nWords-1] == 0 )
  1877.         nWords--;
  1878.     }
  1879.  
  1880.     /*
  1881.      * normalize this number by shifting until
  1882.      * the MSB of the number is at 0x08000000.
  1883.      * This is in preparation for quoRemIteration, below.
  1884.      * The idea is that, to make division easier, we want the
  1885.      * divisor to be "normalized" -- usually this means shifting
  1886.      * the MSB into the high words sign bit. But because we know that
  1887.      * the quotient will be 0 < q < 10, we would like to arrange that
  1888.      * the dividend not span up into another word of precision.
  1889.      * (This needs to be explained more clearly!)
  1890.      */
  1891.     public int
  1892.     normalizeMe() throws IllegalArgumentException {
  1893.     int src;
  1894.     int wordcount = 0;
  1895.     int bitcount  = 0;
  1896.     int v = 0;
  1897.     for ( src= nWords-1 ; src >= 0 && (v=data[src]) == 0 ; src--){
  1898.         wordcount += 1;
  1899.     }
  1900.     if ( src < 0 ){
  1901.         // oops. Value is zero. Cannot normalize it!
  1902.         throw new IllegalArgumentException("zero value");
  1903.     }
  1904.     /*
  1905.      * In most cases, we assume that wordcount is zero. This only
  1906.      * makes sense, as we try not to maintain any high-order
  1907.      * words full of zeros. In fact, if there are zeros, we will
  1908.      * simply SHORTEN our number at this point. Watch closely...
  1909.      */
  1910.     nWords -= wordcount;
  1911.     /*
  1912.      * Compute how far left we have to shift v s.t. its highest-
  1913.      * order bit is in the right place. Then call lshiftMe to
  1914.      * do the work.
  1915.      */
  1916.     if ( (v & 0xf0000000) != 0 ){
  1917.         // will have to shift up into the next word.
  1918.         // too bad.
  1919.         for( bitcount = 32 ; (v & 0xf0000000) != 0 ; bitcount-- )
  1920.         v >>>= 1;
  1921.     } else {
  1922.         while ( v <= 0x000fffff ){
  1923.         // hack: byte-at-a-time shifting
  1924.         v <<= 8;
  1925.         bitcount += 8;
  1926.         }
  1927.         while ( v <= 0x07ffffff ){
  1928.         v <<= 1;
  1929.         bitcount += 1;
  1930.         }
  1931.     }
  1932.     if ( bitcount != 0 )
  1933.         lshiftMe( bitcount );
  1934.     return bitcount;
  1935.     }
  1936.  
  1937.     /*
  1938.      * Multiply a FDBigInt by an int.
  1939.      * Result is a new FDBigInt.
  1940.      */
  1941.     public FDBigInt
  1942.     mult( int iv ) {
  1943.     long v = iv;
  1944.     int r[];
  1945.     long p;
  1946.  
  1947.     // guess adequate size of r.
  1948.     r = new int[ ( v * ((long)data[nWords-1]&0xffffffffL) > 0xfffffffL ) ? nWords+1 : nWords ];
  1949.     p = 0L;
  1950.     for( int i=0; i < nWords; i++ ) {
  1951.         p += v * ((long)data[i]&0xffffffffL);
  1952.         r[i] = (int)p;
  1953.         p >>>= 32;
  1954.     }
  1955.     if ( p == 0L){
  1956.         return new FDBigInt( r, nWords );
  1957.     } else {
  1958.         r[nWords] = (int)p;
  1959.         return new FDBigInt( r, nWords+1 );
  1960.     }
  1961.     }
  1962.  
  1963.     /*
  1964.      * Multiply a FDBigInt by an int and add another int.
  1965.      * Result is computed in place.
  1966.      * Hope it fits!
  1967.      */
  1968.     public void
  1969.     multaddMe( int iv, int addend ) {
  1970.     long v = iv;
  1971.     long p;
  1972.  
  1973.     // unroll 0th iteration, doing addition.
  1974.     p = v * ((long)data[0]&0xffffffffL) + ((long)addend&0xffffffffL);
  1975.     data[0] = (int)p;
  1976.     p >>>= 32;
  1977.     for( int i=1; i < nWords; i++ ) {
  1978.         p += v * ((long)data[i]&0xffffffffL);
  1979.         data[i] = (int)p;
  1980.         p >>>= 32;
  1981.     }
  1982.     if ( p != 0L){
  1983.         data[nWords] = (int)p; // will fail noisily if illegal!
  1984.         nWords++;
  1985.     }
  1986.     }
  1987.  
  1988.     /*
  1989.      * Multiply a FDBigInt by another FDBigInt.
  1990.      * Result is a new FDBigInt.
  1991.      */
  1992.     public FDBigInt
  1993.     mult( FDBigInt other ){
  1994.     // crudely guess adequate size for r
  1995.     int r[] = new int[ nWords + other.nWords ];
  1996.     int i;
  1997.     // I think I am promised zeros...
  1998.  
  1999.     for( i = 0; i < this.nWords; i++ ){
  2000.         long v = (long)this.data[i] & 0xffffffffL; // UNSIGNED CONVERSION
  2001.         long p = 0L;
  2002.         int j;
  2003.         for( j = 0; j < other.nWords; j++ ){
  2004.         p += ((long)r[i+j]&0xffffffffL) + v*((long)other.data[j]&0xffffffffL); // UNSIGNED CONVERSIONS ALL 'ROUND.
  2005.         r[i+j] = (int)p;
  2006.         p >>>= 32;
  2007.         }
  2008.         r[i+j] = (int)p;
  2009.     }
  2010.     // compute how much of r we actually needed for all that.
  2011.     for ( i = r.length-1; i> 0; i--)
  2012.         if ( r[i] != 0 )
  2013.         break;
  2014.     return new FDBigInt( r, i+1 );
  2015.     }
  2016.  
  2017.     /*
  2018.      * Add one FDBigInt to another. Return a FDBigInt
  2019.      */
  2020.     public FDBigInt
  2021.     add( FDBigInt other ){
  2022.     int i;
  2023.     int a[], b[];
  2024.     int n, m;
  2025.     long c = 0L;
  2026.     // arrange such that a.nWords >= b.nWords;
  2027.     // n = a.nWords, m = b.nWords
  2028.     if ( this.nWords >= other.nWords ){
  2029.         a = this.data;
  2030.         n = this.nWords;
  2031.         b = other.data;
  2032.         m = other.nWords;
  2033.     } else {
  2034.         a = other.data;
  2035.         n = other.nWords;
  2036.         b = this.data;
  2037.         m = this.nWords;
  2038.     }
  2039.     int r[] = new int[ n ];
  2040.     for ( i = 0; i < n; i++ ){
  2041.         c += (long)a[i] & 0xffffffffL;
  2042.         if ( i < m ){
  2043.         c += (long)b[i] & 0xffffffffL;
  2044.         }
  2045.         r[i] = (int) c;
  2046.         c >>= 32; // signed shift.
  2047.     }
  2048.     if ( c != 0L ){
  2049.         // oops -- carry out -- need longer result.
  2050.         int s[] = new int[ r.length+1 ];
  2051.         System.arraycopy( r, 0, s, 0, r.length );
  2052.         s[i++] = (int)c;
  2053.         return new FDBigInt( s, i );
  2054.     }
  2055.     return new FDBigInt( r, i );
  2056.     }
  2057.  
  2058.     /*
  2059.      * Subtract one FDBigInt from another. Return a FDBigInt
  2060.      * Assert that the result is positive.
  2061.      */
  2062.     public FDBigInt
  2063.     sub( FDBigInt other ){
  2064.     int r[] = new int[ this.nWords ];
  2065.     int i;
  2066.     int n = this.nWords;
  2067.     int m = other.nWords;
  2068.     int nzeros = 0;
  2069.     long c = 0L;
  2070.     for ( i = 0; i < n; i++ ){
  2071.         c += (long)this.data[i] & 0xffffffffL;
  2072.         if ( i < m ){
  2073.         c -= (long)other.data[i] & 0xffffffffL;
  2074.         }
  2075.         if ( ( r[i] = (int) c ) == 0 )
  2076.         nzeros++;
  2077.         else
  2078.         nzeros = 0;
  2079.         c >>= 32; // signed shift.
  2080.     }
  2081.     if ( c != 0L )
  2082.         throw new RuntimeException("Assertion botch: borrow out of subtract");
  2083.     while ( i < m )
  2084.         if ( other.data[i++] != 0 )
  2085.         throw new RuntimeException("Assertion botch: negative result of subtract");
  2086.     return new FDBigInt( r, n-nzeros );
  2087.     }
  2088.  
  2089.     /*
  2090.      * Compare FDBigInt with another FDBigInt. Return an integer
  2091.      * >0: this > other
  2092.      *  0: this == other
  2093.      * <0: this < other
  2094.      */
  2095.     public int
  2096.     cmp( FDBigInt other ){
  2097.         int i;
  2098.     if ( this.nWords > other.nWords ){
  2099.         // if any of my high-order words is non-zero,
  2100.         // then the answer is evident
  2101.         int j = other.nWords-1;
  2102.         for ( i = this.nWords-1; i > j ; i-- )
  2103.         if ( this.data[i] != 0 ) return 1;
  2104.     }else if ( this.nWords < other.nWords ){
  2105.         // if any of other's high-order words is non-zero,
  2106.         // then the answer is evident
  2107.         int j = this.nWords-1;
  2108.         for ( i = other.nWords-1; i > j ; i-- )
  2109.         if ( other.data[i] != 0 ) return -1;
  2110.     } else{
  2111.         i = this.nWords-1;
  2112.     }
  2113.     for ( ; i > 0 ; i-- )
  2114.         if ( this.data[i] != other.data[i] )
  2115.         break;
  2116.     // careful! want unsigned compare!
  2117.     // use brute force here.
  2118.     int a = this.data[i];
  2119.     int b = other.data[i];
  2120.     if ( a < 0 ){
  2121.         // a is really big, unsigned
  2122.         if ( b < 0 ){
  2123.         return a-b; // both big, negative
  2124.         } else {
  2125.         return 1; // b not big, answer is obvious;
  2126.         }
  2127.     } else {
  2128.         // a is not really big
  2129.         if ( b < 0 ) {
  2130.         // but b is really big
  2131.         return -1;
  2132.         } else {
  2133.         return a - b;
  2134.         }
  2135.     }
  2136.     }
  2137.  
  2138.     /*
  2139.      * Compute
  2140.      * q = (int)( this / S )
  2141.      * this = 10 * ( this mod S )
  2142.      * Return q.
  2143.      * This is the iteration step of digit development for output.
  2144.      * We assume that S has been normalized, as above, and that
  2145.      * "this" has been lshift'ed accordingly.
  2146.      * Also assume, of course, that the result, q, can be expressed
  2147.      * as an integer, 0 <= q < 10.
  2148.      */
  2149.     public int
  2150.     quoRemIteration( FDBigInt S )throws IllegalArgumentException {
  2151.     // ensure that this and S have the same number of
  2152.     // digits. If S is properly normalized and q < 10 then
  2153.     // this must be so.
  2154.     if ( nWords != S.nWords ){
  2155.         throw new IllegalArgumentException("disparate values");
  2156.     }
  2157.     // estimate q the obvious way. We will usually be
  2158.     // right. If not, then we're only off by a little and
  2159.     // will re-add.
  2160.     int n = nWords-1;
  2161.     long q = ((long)data[n]&0xffffffffL) / (long)S.data[n];
  2162.     long diff = 0L;
  2163.     for ( int i = 0; i <= n ; i++ ){
  2164.         diff += ((long)data[i]&0xffffffffL) -  q*((long)S.data[i]&0xffffffffL);
  2165.         data[i] = (int)diff;
  2166.         diff >>= 32; // N.B. SIGNED shift.
  2167.     }
  2168.     if ( diff != 0L ) {
  2169.         // damn, damn, damn. q is too big.
  2170.         // add S back in until this turns +. This should
  2171.         // not be very many times!
  2172.         long sum = 0L;
  2173.         while ( sum ==  0L ){
  2174.         sum = 0L;
  2175.         for ( int i = 0; i <= n; i++ ){
  2176.             sum += ((long)data[i]&0xffffffffL) +  ((long)S.data[i]&0xffffffffL);
  2177.             data[i] = (int) sum;
  2178.             sum >>= 32; // Signed or unsigned, answer is 0 or 1
  2179.         }
  2180.         /*
  2181.          * Originally the following line read
  2182.          * "if ( sum !=0 && sum != -1 )"
  2183.          * but that would be wrong, because of the
  2184.          * treatment of the two values as entirely unsigned,
  2185.          * it would be impossible for a carry-out to be interpreted
  2186.          * as -1 -- it would have to be a single-bit carry-out, or
  2187.          * +1.
  2188.          */
  2189.         if ( sum !=0 && sum != 1 )
  2190.             throw new RuntimeException("Assertion botch: "+sum+" carry out of division correction");
  2191.         q -= 1;
  2192.         }
  2193.     }
  2194.     // finally, we can multiply this by 10.
  2195.     // it cannot overflow, right, as the high-order word has
  2196.     // at least 4 high-order zeros!
  2197.     long p = 0L;
  2198.     for ( int i = 0; i <= n; i++ ){
  2199.         p += 10*((long)data[i]&0xffffffffL);
  2200.         data[i] = (int)p;
  2201.         p >>= 32; // SIGNED shift.
  2202.     }
  2203.     if ( p != 0L )
  2204.         throw new RuntimeException("Assertion botch: carry out of *10");
  2205.  
  2206.     return (int)q;
  2207.     }
  2208.  
  2209.     public long
  2210.     longValue(){
  2211.     // if this can be represented as a long,
  2212.     // return the value
  2213.     int i;
  2214.     for ( i = this.nWords-1; i > 1 ; i-- ){
  2215.         if ( data[i] != 0 ){
  2216.         throw new RuntimeException("Assertion botch: value too big");
  2217.         }
  2218.     }
  2219.     switch(i){
  2220.     case 1:
  2221.         if ( data[1] < 0 )
  2222.         throw new RuntimeException("Assertion botch: value too big");
  2223.         return ((long)(data[1]) << 32) | ((long)data[0]&0xffffffffL);
  2224.     case 0:
  2225.         return ((long)data[0]&0xffffffffL);
  2226.     default:
  2227.         throw new RuntimeException("Assertion botch: longValue confused");
  2228.     }
  2229.     }
  2230.  
  2231.     public String
  2232.     toString() {
  2233.     StringBuffer r = new StringBuffer(30);
  2234.     r.append('[');
  2235.     int i = Math.min( nWords-1, data.length-1) ;
  2236.     if ( nWords > data.length ){
  2237.         r.append( "("+data.length+"<"+nWords+"!)" );
  2238.     }
  2239.     for( ; i> 0 ; i-- ){
  2240.         r.append( Integer.toHexString( data[i] ) );
  2241.         r.append( (char) ' ' );
  2242.     }
  2243.     r.append( Integer.toHexString( data[0] ) );
  2244.     r.append( (char) ']' );
  2245.     return new String( r );
  2246.     }
  2247. }
  2248.