home *** CD-ROM | disk | FTP | other *** search
/ Chip 1998 November / Chip_1998-11_cd.bin / tema / Cafe / main.bin / FloatingDecimal.java < prev    next >
Text File  |  1997-05-20  |  39KB  |  1,376 lines
  1. /*
  2.  * @(#)FloatingDecimal.java    1.9 97/01/22
  3.  * 
  4.  * Copyright (c) 1995, 1996 Sun Microsystems, Inc. All Rights Reserved.
  5.  * 
  6.  * This software is the confidential and proprietary information of Sun
  7.  * Microsystems, Inc. ("Confidential Information").  You shall not
  8.  * disclose such Confidential Information and shall use it only in
  9.  * accordance with the terms of the license agreement you entered into
  10.  * with Sun.
  11.  * 
  12.  * SUN MAKES NO REPRESENTATIONS OR WARRANTIES ABOUT THE SUITABILITY OF THE
  13.  * SOFTWARE, EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
  14.  * IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
  15.  * PURPOSE, OR NON-INFRINGEMENT. SUN SHALL NOT BE LIABLE FOR ANY DAMAGES
  16.  * SUFFERED BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
  17.  * THIS SOFTWARE OR ITS DERIVATIVES.
  18.  * 
  19.  * CopyrightVersion 1.1_beta
  20.  * 
  21.  */
  22.  
  23. package java.lang;
  24.  
  25. class FloatingDecimal{
  26.     boolean    isExceptional;
  27.     boolean    isNegative;
  28.     int        decExponent;
  29.     char    digits[];
  30.     int        nDigits;
  31.  
  32.     private    FloatingDecimal( boolean negSign, int decExponent, char []digits, int n,  boolean e )
  33.     {
  34.     isNegative = negSign;
  35.     isExceptional = e;
  36.     this.decExponent = decExponent;
  37.     this.digits = digits;
  38.     this.nDigits = n;
  39.     }
  40.  
  41.     /*
  42.      * Constants of the implementation
  43.      * Most are IEEE-754 related.
  44.      * (There are more really boring constants at the end.)
  45.      */
  46.     static final long    signMask = 0x8000000000000000L;
  47.     static final long    expMask  = 0x7ff0000000000000L;
  48.     static final long    fractMask= ~(signMask|expMask);
  49.     static final int    expShift = 52;
  50.     static final int    expBias  = 1023;
  51.     static final long    fractHOB = ( 1L<<expShift ); // assumed High-Order bit
  52.     static final long    expOne     = ((long)expBias)<<expShift; // exponent of 1.0
  53.     static final int    maxSmallBinExp = 62;
  54.     static final int    minSmallBinExp = -( 63 / 3 );
  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.  
  67.     /*
  68.      * count number of bits from high-order 1 bit to low-order 1 bit,
  69.      * inclusive.
  70.      */
  71.     private static int
  72.     countBits( long v ){
  73.     // 
  74.     // the strategy is to shift until we get a non-zero sign bit
  75.     // then shift until we have no bits left, counting the difference.
  76.     // we do byte shifting as a hack. Hope it helps.
  77.     //
  78.     if ( v == 0L ) return 0;
  79.  
  80.     while ( ( v & highbyte ) == 0L ){
  81.         v <<= 8;
  82.     }
  83.     while ( v > 0L ) { // i.e. while ((v&highbit) == 0L )
  84.         v <<= 1;
  85.     }
  86.  
  87.     int n = 0;
  88.     while (( v & lowbytes ) != 0L ){
  89.         v <<= 8;
  90.         n += 8;
  91.     }
  92.     while ( v != 0L ){
  93.         v <<= 1;
  94.         n += 1;
  95.     }
  96.     return n;
  97.     }
  98.  
  99.     /*
  100.      * Keep big powers of 5 handy for future reference.
  101.      */
  102.     private static FDBigInt b5p[];
  103.  
  104.     private static FDBigInt
  105.     big5pow( int p ){
  106.     if ( p < 0 )
  107.         throw new RuntimeException( "Assertion botch: negative power of 5");
  108.     if ( b5p == null ){
  109.         b5p = new FDBigInt[ p+1 ];
  110.     }else if (b5p.length <= p ){
  111.         FDBigInt t[] = new FDBigInt[ p+1 ];
  112.         System.arraycopy( b5p, 0, t, 0, b5p.length );
  113.         b5p = t;
  114.     }
  115.     if ( b5p[p] != null )
  116.         return b5p[p];
  117.     else if ( p < small5pow.length )
  118.         return b5p[p] = new FDBigInt( small5pow[p] );
  119.     else if ( p < long5pow.length )
  120.         return b5p[p] = new FDBigInt( long5pow[p] );
  121.     else {
  122.         // construct the damn thing.
  123.         // recursively.
  124.         int q, r;
  125.         // in order to compute 5^p,
  126.         // compute its square root, 5^(p/2) and square.
  127.         // or, let q = p / 2, r = p -q, then
  128.         // 5^p = 5^(q+r) = 5^q * 5^r
  129.         q = p >> 1;
  130.         r = p - q;
  131.         FDBigInt bigq =  b5p[q];
  132.         if ( bigq == null ) 
  133.         bigq = big5pow ( q );
  134.         if ( r < small5pow.length ){
  135.         return (b5p[p] = bigq.mult( small5pow[r] ) );
  136.         }else{
  137.         FDBigInt bigr = b5p[ r ];
  138.         if ( bigr == null ) 
  139.             bigr = big5pow( r );
  140.         return (b5p[p] = bigq.mult( bigr ) );
  141.         }
  142.     }
  143.     }
  144.  
  145.     /*
  146.      * This is the easy subcase -- 
  147.      * all the significant bits, after scaling, are held in lvalue.
  148.      * negSign and decExponent tell us what processing and scaling
  149.      * has already been done. Exceptional cases have already been
  150.      * stripped out. 
  151.      * In particular:
  152.      * lvalue is a finite number (not Inf, nor NaN)
  153.      * lvalue > 0L (not zero, nor negative).
  154.      *
  155.      * The only reason that we develop the digits here, rather than
  156.      * calling on Long.toString() is that we can do it a little faster,
  157.      * and besides want to treat trailing 0s specially. If Long.toString
  158.      * changes, we should re-evaluate this strategy!
  159.      */
  160.     private void
  161.     developLongDigits( int decExponent, long lvalue, long insignificant ){
  162.     char digits[];
  163.     int  ndigits;
  164.     int  digitno;
  165.     int  c;
  166.     //
  167.     // Discard non-significant low-order bits, while rounding,
  168.     // up to insignificant value.
  169.     int i;
  170.     for ( i = 0; insignificant >= 10L; i++ )
  171.         insignificant /= 10L;
  172.     if ( i != 0 ){
  173.         long pow10 = long5pow[i] << i; // 10^i == 5^i * 2^i;
  174.         long residue = lvalue % pow10;
  175.         lvalue /= pow10;
  176.         decExponent += i;
  177.         if ( residue >= (pow10>>1) ){
  178.         // round up based on the low-order bits we're discarding
  179.         lvalue++;
  180.         }
  181.     }
  182.     if ( lvalue <= Integer.MAX_VALUE ){
  183.         if ( lvalue <= 0L )
  184.         throw new RuntimeException("Assertion botch: value "+lvalue+" <= 0");
  185.  
  186.         // even easier subcase!
  187.         // can do int arithmetic rather than long!
  188.         int  ivalue = (int)lvalue;
  189.         digits = new char[ ndigits=10 ];
  190.         digitno = ndigits-1;
  191.         c = ivalue%10;
  192.         ivalue /= 10;
  193.         while ( c == 0 ){
  194.         decExponent++;
  195.         c = ivalue%10;
  196.         ivalue /= 10;
  197.         }
  198.         while ( ivalue != 0){
  199.         digits[digitno--] = (char)(c+'0');
  200.         decExponent++;
  201.         c = ivalue%10;
  202.         ivalue /= 10;
  203.         }
  204.         digits[digitno] = (char)(c+'0');
  205.     } else {
  206.         // same algorithm as above (same bugs, too )
  207.         // but using long arithmetic.
  208.         digits = new char[ ndigits=20 ];
  209.         digitno = ndigits-1;
  210.         c = (int)(lvalue%10L);
  211.         lvalue /= 10L;
  212.         while ( c == 0 ){
  213.         decExponent++;
  214.         c = (int)(lvalue%10L);
  215.         lvalue /= 10L;
  216.         }
  217.         while ( lvalue != 0L ){
  218.         digits[digitno--] = (char)(c+'0');
  219.         decExponent++;
  220.         c = (int)(lvalue%10L);
  221.         lvalue /= 10;
  222.         }
  223.         digits[digitno] = (char)(c+'0');
  224.     }
  225.     char result [];
  226.     ndigits -= digitno;
  227.     if ( digitno == 0 )
  228.         result = digits;
  229.     else {
  230.         result = new char[ ndigits ];
  231.         System.arraycopy( digits, digitno, result, 0, ndigits );
  232.     }
  233.     this.digits = result;
  234.     this.decExponent = decExponent+1;
  235.     this.nDigits = ndigits;
  236.     }
  237.  
  238.     //
  239.     // add one to the least significant digit.
  240.     // in the unlikely event there is a carry out,
  241.     // deal with it.
  242.     // assert that this will only happen where there
  243.     // is only one digit, e.g. (float)1e-44 seems to do it.
  244.     //
  245.     private void
  246.     roundup(){
  247.     int i;
  248.     int q = digits[ i = (nDigits-1)];
  249.     if ( q == '9' ){
  250.         while ( q == '9' && i > 0 ){
  251.         digits[i] = '0';
  252.         q = digits[--i];
  253.         }
  254.         if ( q == '9' ){
  255.         // carryout! High-order 1, rest 0s, larger exp.
  256.         decExponent += 1;
  257.         digits[0] = '1';
  258.         return;
  259.         }
  260.         // else fall through.
  261.     }
  262.     digits[i] = (char)(q+1);
  263.     }
  264.  
  265.     /*
  266.      * FIRST IMPORTANT CONSTRUCTOR: DOUBLE
  267.      */
  268.     public FloatingDecimal( double d )
  269.     {
  270.     long    dBits = Double.doubleToLongBits( d );
  271.     long    fractBits;
  272.     int    binExp;
  273.     int    nSignificantBits;
  274.  
  275.     // discover and delete sign
  276.     if ( (dBits&signMask) != 0 ){
  277.         isNegative = true;
  278.         dBits ^= signMask;
  279.     } else {
  280.         isNegative = false;
  281.     }
  282.     // Begin to unpack
  283.     // Discover obvious special cases of NaN and Infinity.
  284.     binExp = (int)( (dBits&expMask) >> expShift );
  285.     fractBits = dBits&fractMask;
  286.     if ( binExp == (int)(expMask>>expShift) ) {
  287.         isExceptional = true;
  288.         if ( fractBits == 0L ){
  289.         digits =  infinity;
  290.         } else {
  291.         digits = notANumber;
  292.         isNegative = false; // NaN has no sign!
  293.         }
  294.         nDigits = digits.length;
  295.         return;
  296.     }
  297.     isExceptional = false;
  298.     // Finish unpacking
  299.     // Normalize denormalized numbers.
  300.     // Insert assumed high-order bit for normalized numbers.
  301.     // Subtract exponent bias.
  302.     if ( binExp == 0 ){
  303.         if ( fractBits == 0L ){
  304.         // not a denorm, just a 0!
  305.         decExponent = 0;
  306.         digits = zero;
  307.         nDigits = 1;
  308.         return;
  309.         }
  310.         while ( (fractBits&fractHOB) == 0L ){
  311.         fractBits <<= 1;
  312.         binExp -= 1;
  313.         }
  314.         nSignificantBits = expShift + binExp; // recall binExp is  - shift count.
  315.         binExp += 1;
  316.     } else {
  317.         fractBits |= fractHOB;
  318.         nSignificantBits = expShift+1;
  319.     }
  320.     binExp -= expBias;
  321.     // call the routine that actually does all the hard work.
  322.     dtoa( binExp, fractBits, nSignificantBits );
  323.     }
  324.  
  325.     /*
  326.      * SECOND IMPORTANT CONSTRUCTOR: SINGLE
  327.      */
  328.     public FloatingDecimal( float f )
  329.     {
  330.     int    fBits = Float.floatToIntBits( f );
  331.     int    fractBits;
  332.     int    binExp;
  333.     int    nSignificantBits;
  334.  
  335.     // discover and delete sign
  336.     if ( (fBits&singleSignMask) != 0 ){
  337.         isNegative = true;
  338.         fBits ^= singleSignMask;
  339.     } else {
  340.         isNegative = false;
  341.     }
  342.     // Begin to unpack
  343.     // Discover obvious special cases of NaN and Infinity.
  344.     binExp = (int)( (fBits&singleExpMask) >> singleExpShift );
  345.     fractBits = fBits&singleFractMask;
  346.     if ( binExp == (int)(singleExpMask>>singleExpShift) ) {
  347.         isExceptional = true;
  348.         if ( fractBits == 0L ){
  349.         digits =  infinity;
  350.         } else {
  351.         digits = notANumber;
  352.         isNegative = false; // NaN has no sign!
  353.         }
  354.         nDigits = digits.length;
  355.         return;
  356.     }
  357.     isExceptional = false;
  358.     // Finish unpacking
  359.     // Normalize denormalized numbers.
  360.     // Insert assumed high-order bit for normalized numbers.
  361.     // Subtract exponent bias.
  362.     if ( binExp == 0 ){
  363.         if ( fractBits == 0 ){
  364.         // not a denorm, just a 0!
  365.         decExponent = 0;
  366.         digits = zero;
  367.         nDigits = 1;
  368.         return;
  369.         }
  370.         while ( (fractBits&singleFractHOB) == 0 ){
  371.         fractBits <<= 1;
  372.         binExp -= 1;
  373.         }
  374.         nSignificantBits = singleExpShift + binExp; // recall binExp is  - shift count.
  375.         binExp += 1;
  376.     } else {
  377.         fractBits |= singleFractHOB;
  378.         nSignificantBits = singleExpShift+1;
  379.     }
  380.     binExp -= singleExpBias;
  381.     // call the routine that actually does all the hard work.
  382.     dtoa( binExp, ((long)fractBits)<<(expShift-singleExpShift), nSignificantBits );
  383.     }
  384.  
  385.     private void
  386.     dtoa( int binExp, long fractBits, int nSignificantBits )
  387.     {
  388.     int    nFractBits; // number of significant bits of fractBits;
  389.     int    nTinyBits;  // number of these to the right of the point.
  390.     int    decExp;
  391.  
  392.     // Examine number. Determine if it is an easy case,
  393.     // which we can do pretty trivially using float/long conversion,
  394.     // or whether we must do real work.
  395.     nFractBits = countBits( fractBits );
  396.     nTinyBits = Math.max( 0, nFractBits - binExp - 1 );
  397.     if ( binExp <= maxSmallBinExp && binExp >= minSmallBinExp ){
  398.         // Look more closely at the number to decide if,
  399.         // with scaling by 10^nTinyBits, the result will fit in
  400.         // a long.
  401.         if ( (nTinyBits < long5pow.length) && ((nFractBits + n5bits[nTinyBits]) < 64 ) ){
  402.         /*
  403.          * We can do this:
  404.          * take the fraction bits, which are normalized.
  405.          * (a) nTinyBits == 0: Shift left or right appropriately
  406.          *     to align the binary point at the extreme right, i.e.
  407.          *     where a long int point is expected to be. The integer
  408.          *     result is easily converted to a string.
  409.          * (b) nTinyBits > 0: Shift right by expShift-nFractBits,
  410.          *     which effectively converts to long and scales by
  411.          *     2^nTinyBits. Then multiply by 5^nTinyBits to
  412.          *     complete the scaling. We know this won't overflow
  413.          *     because we just counted the number of bits necessary
  414.          *     in the result. The integer you get from this can
  415.          *     then be converted to a string pretty easily.
  416.          */
  417.         long halfULP;
  418.         if ( nTinyBits == 0 ) {
  419.             if ( binExp > nSignificantBits ){
  420.             halfULP = 1L << ( binExp-nSignificantBits-1);
  421.             } else {
  422.             halfULP = 0L;
  423.             }
  424.             if ( binExp >= expShift ){
  425.             fractBits <<= (binExp-expShift);
  426.             } else {
  427.             fractBits >>>= (expShift-binExp) ;
  428.             }
  429.             developLongDigits( 0, fractBits, halfULP );
  430.             return;
  431.         }
  432.         /*
  433.          * The following causes excess digits to be printed
  434.          * out in the single-float case. Our manipulation of
  435.          * halfULP here is apparently not correct. If we
  436.          * better understand how this works, perhaps we can
  437.          * use this special case again. But for the time being,
  438.          * we do not.
  439.          * else {
  440.          *     fractBits >>>= expShift+1-nFractBits;
  441.          *     fractBits *= long5pow[ nTinyBits ];
  442.          *     halfULP = long5pow[ nTinyBits ] >> (1+nSignificantBits-nFractBits);
  443.          *     developLongDigits( -nTinyBits, fractBits, halfULP );
  444.          *     return;
  445.          * }
  446.          */
  447.         }
  448.     }
  449.     /*
  450.      * This is the hard case. We are going to compute large positive
  451.      * integers B and S and integer decExp, s.t.
  452.      *    d = ( B / S ) * 10^decExp
  453.      *    1 <= B / S < 10
  454.      * Obvious choices are:
  455.      *    decExp = floor( log10(d) )
  456.      *     B      = d * 2^nTinyBits * 10^max( 0, -decExp )
  457.      *    S      = 10^max( 0, decExp) * 2^nTinyBits
  458.      * (noting that nTinyBits has already been forced to non-negative)
  459.      * I am also going to compute a large positive integer
  460.      *    M      = (1/2^nSignificantBits) * 2^nTinyBits * 10^max( 0, -decExp )
  461.      * i.e. M is (1/2) of the ULP of d, scaled like B.
  462.      * When we iterate through dividing B/S and picking off the
  463.      * quotient bits, we will know when to stop when the remainder
  464.      * is <= M.
  465.      *
  466.      * We keep track of powers of 2 and powers of 5.
  467.      */
  468.  
  469.     /*
  470.      * Estimate decimal exponent. (If it is small-ish,
  471.      * we could double-check.)
  472.      *
  473.      * First, scale the mantissa bits such that 1 <= d2 < 2.
  474.      * We are then going to estimate
  475.      *        log10(d2) ~=~  (d2-1.5)/1.5 + log(1.5) 
  476.      * and so we can estimate 
  477.      *      log10(d) ~=~ log10(d2) + binExp * log10(2)
  478.      * take the floor and call it decExp.
  479.      * FIXME -- use more precise constants here. It costs no more.
  480.      */
  481.     double d2 = Double.longBitsToDouble( 
  482.         expOne | ( fractBits &~ fractHOB ) );
  483.     decExp = (int)Math.floor(
  484.         (d2-1.5D)*0.289529654D + 0.176091259 + (double)binExp * 0.301029995663981 );
  485.     int B2, B5; // powers of 2 and powers of 5, respectively, in B
  486.     int S2, S5; // powers of 2 and powers of 5, respectively, in S
  487.     int M2, M5; // powers of 2 and powers of 5, respectively, in M
  488.     int Bbits; // binary digits needed to represent B, approx.
  489.     int tenSbits; // binary digits needed to represent 10*S, approx.
  490.     FDBigInt Sval, Bval, Mval;
  491.  
  492.     B5 = Math.max( 0, -decExp );
  493.     B2 = B5 + nTinyBits + binExp;
  494.  
  495.     S5 = Math.max( 0, decExp );
  496.     S2 = S5 + nTinyBits;
  497.  
  498.     M5 = B5;
  499.     M2 = B2 - nSignificantBits;
  500.  
  501.     /*
  502.      * the long integer fractBits contains the (nFractBits) interesting
  503.      * bits from the mantissa of d ( hidden 1 added if necessary) followed
  504.      * by (expShift+1-nFractBits) zeros. In the interest of compactness,
  505.      * I will shift out those zeros before turning fractBits into a
  506.      * FDBigInt. The resulting whole number will be 
  507.      *     d * 2^(nFractBits-1-binExp).
  508.      */
  509.     fractBits >>>= (expShift+1-nFractBits);
  510.     B2 -= nFractBits-1;
  511.     int common2factor = Math.min( B2, S2 );
  512.     B2 -= common2factor;
  513.     S2 -= common2factor;
  514.     M2 -= common2factor;
  515.  
  516.     /*
  517.      * HACK!! For exact powers of two, the next smallest number
  518.      * is only half as far away as we think (because the meaning of
  519.      * ULP changes at power-of-two bounds) for this reason, we
  520.      * hack M2. Hope this works.
  521.      */
  522.     if ( nFractBits == 1 )
  523.         M2 -= 1;
  524.  
  525.     if ( M2 < 0 ){
  526.         // oops. 
  527.         // since we cannot scale M down far enough,
  528.         // we must scale the other values up.
  529.         B2 -= M2;
  530.         S2 -= M2;
  531.         M2 =  0;
  532.     }
  533.     /*
  534.      * Construct, Scale, iterate.
  535.      * Some day, we'll write a stopping test that takes
  536.      * account of the assymetry of the spacing of floating-point
  537.      * numbers below perfect powers of 2
  538.      * 26 Sept 96 is not that day.
  539.      * So we use a symmetric test.
  540.      */
  541.     char digits[] = this.digits = new char[18];
  542.     int  ndigit = 0;
  543.     boolean low, high;
  544.     long lowDigitDifference;
  545.     int  q;
  546.  
  547.     /*
  548.      * Detect the special cases where all the numbers we are about
  549.      * to compute will fit in int or long integers.
  550.      * In these cases, we will avoid doing FDBigInt arithmetic.
  551.      * We use the same algorithms, except that we "normalize"
  552.      * our FDBigInts before iterating. This is to make division easier,
  553.      * as it makes our fist guess (quotient of high-order words)
  554.      * more accurate!
  555.      *
  556.      * Some day, we'll write a stopping test that takes
  557.      * account of the assymetry of the spacing of floating-point
  558.      * numbers below perfect powers of 2
  559.      * 26 Sept 96 is not that day.
  560.      * So we use a symmetric test.
  561.      */
  562.     Bbits = nFractBits + B2 + (( B5 < n5bits.length )? n5bits[B5] : ( B5*3 ));
  563.     tenSbits = S2+1 + (( (S5+1) < n5bits.length )? n5bits[(S5+1)] : ( (S5+1)*3 ));
  564.     if ( Bbits < 64 && tenSbits < 64){
  565.         if ( Bbits < 32 && tenSbits < 32){
  566.         // wa-hoo! They're all ints!
  567.         int b = ((int)fractBits * small5pow[B5] ) << B2;
  568.         int s = small5pow[S5] << S2;
  569.         int m = small5pow[M5] << M2;
  570.         int tens = s * 10;
  571.         /*
  572.          * Unroll the first iteration. If our decExp estimate
  573.          * was too high, our first quotient will be zero. In this
  574.          * case, we discard it and decrement decExp.
  575.          */
  576.         ndigit = 0;
  577.         q = (int) ( b / s );
  578.         b = 10 * ( b % s );
  579.         m *= 10;
  580.         low  = (b <  m );
  581.         high = (b+m > tens );
  582.         if ( q >= 10 ){
  583.             // bummer, dude
  584.             throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
  585.         } else if ( (q == 0) && ! high ){
  586.             // oops. Usually ignore leading zero.
  587.             decExp--;
  588.         } else {
  589.             digits[ndigit++] = (char)('0' + q);
  590.         }
  591.         /*
  592.          * HACK! Java spec sez that we always have at least
  593.          * one digit after the . in either F- or E-form output.
  594.          * Thus we will need more than one digit if we're using
  595.          * E-form
  596.          */
  597.         if ( decExp <= -3 || decExp >= 8 ){
  598.             high = low = false;
  599.         }
  600.         while( ! low && ! high ){
  601.             q = (int) ( b / s );
  602.             b = 10 * ( b % s );
  603.             m *= 10;
  604.             if ( q >= 10 ){
  605.             // bummer, dude
  606.             throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
  607.             }
  608.             if ( m > 0L ){
  609.             low  = (b <  m );
  610.             high = (b+m > tens );
  611.             } else {
  612.             // hack -- m might overflow!
  613.             // in this case, it is certainly > b,
  614.             // which won't
  615.             // and b+m > tens, too, since that has overflowed
  616.             // either!
  617.             low = true;
  618.             high = true;
  619.             }
  620.             digits[ndigit++] = (char)('0' + q);
  621.         }
  622.         lowDigitDifference = (b<<1) - tens;
  623.         } else {
  624.         // still good! they're all longs!
  625.         long b = (fractBits * long5pow[B5] ) << B2;
  626.         long s = long5pow[S5] << S2;
  627.         long m = long5pow[M5] << M2;
  628.         long tens = s * 10L;
  629.         /*
  630.          * Unroll the first iteration. If our decExp estimate
  631.          * was too high, our first quotient will be zero. In this
  632.          * case, we discard it and decrement decExp.
  633.          */
  634.         ndigit = 0;
  635.         q = (int) ( b / s );
  636.         b = 10L * ( b % s );
  637.         m *= 10L;
  638.         low  = (b <  m );
  639.         high = (b+m > tens );
  640.         if ( q >= 10 ){
  641.             // bummer, dude
  642.             throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
  643.         } else if ( (q == 0) && ! high ){
  644.             // oops. Usually ignore leading zero.
  645.             decExp--;
  646.         } else {
  647.             digits[ndigit++] = (char)('0' + q);
  648.         }
  649.         /*
  650.          * HACK! Java spec sez that we always have at least
  651.          * one digit after the . in either F- or E-form output.
  652.          * Thus we will need more than one digit if we're using
  653.          * E-form
  654.          */
  655.         if ( decExp <= -3 || decExp >= 8 ){
  656.             high = low = false;
  657.         }
  658.         while( ! low && ! high ){
  659.             q = (int) ( b / s );
  660.             b = 10 * ( b % s );
  661.             m *= 10;
  662.             if ( q >= 10 ){
  663.             // bummer, dude
  664.             throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
  665.             }
  666.             if ( m > 0L ){
  667.             low  = (b <  m );
  668.             high = (b+m > tens );
  669.             } else {
  670.             // hack -- m might overflow!
  671.             // in this case, it is certainly > b,
  672.             // which won't
  673.             // and b+m > tens, too, since that has overflowed
  674.             // either!
  675.             low = true;
  676.             high = true;
  677.             }
  678.             digits[ndigit++] = (char)('0' + q);
  679.         }
  680.         lowDigitDifference = (b<<1) - tens;
  681.         }
  682.     } else {
  683.         FDBigInt tenSval;
  684.         int  shiftBias;
  685.  
  686.         /*
  687.          * We really must do FDBigInt arithmetic.
  688.          * Fist, construct our FDBigInt initial values.
  689.          */
  690.         Bval = new FDBigInt( fractBits  );
  691.         if ( B5 != 0 ){
  692.         if ( B5 < small5pow.length ){
  693.             Bval = Bval.mult( small5pow[B5] );
  694.         } else {
  695.             Bval = Bval.mult( big5pow( B5 ) );
  696.         }
  697.         }
  698.         if ( B2 != 0 ){
  699.         Bval.lshiftMe( B2 );
  700.         }
  701.         Sval = new FDBigInt( big5pow( S5 ) );
  702.         if ( S2 != 0 ){
  703.         Sval.lshiftMe( S2 );
  704.         }
  705.         Mval = new FDBigInt( big5pow( M5 ) );
  706.         if ( M2 != 0 ){
  707.         Mval.lshiftMe( M2 );
  708.         }
  709.  
  710.  
  711.         // normalize so that division works better
  712.         Bval.lshiftMe( shiftBias = Sval.normalizeMe() );
  713.         Mval.lshiftMe( shiftBias );
  714.         tenSval = Sval.mult( 10 );
  715.         /*
  716.          * Unroll the first iteration. If our decExp estimate
  717.          * was too high, our first quotient will be zero. In this
  718.          * case, we discard it and decrement decExp.
  719.          */
  720.         ndigit = 0;
  721.         q = Bval.quoRemIteration( Sval );
  722.         Mval = Mval.mult( 10 );
  723.         low  = (Bval.cmp( Mval ) < 0);
  724.         high = (Bval.add( Mval ).cmp( tenSval ) > 0 );
  725.         if ( q >= 10 ){
  726.         // bummer, dude
  727.         throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
  728.         } else if ( (q == 0) && ! high ){
  729.         // oops. Usually ignore leading zero.
  730.         decExp--;
  731.         } else {
  732.         digits[ndigit++] = (char)('0' + q);
  733.         }
  734.         /*
  735.          * HACK! Java spec sez that we always have at least
  736.          * one digit after the . in either F- or E-form output.
  737.          * Thus we will need more than one digit if we're using
  738.          * E-form
  739.          */
  740.         if ( decExp <= -3 || decExp >= 8 ){
  741.         high = low = false;
  742.         }
  743.         while( ! low && ! high ){
  744.         q = Bval.quoRemIteration( Sval );
  745.         Mval = Mval.mult( 10 );
  746.         if ( q >= 10 ){
  747.             // bummer, dude
  748.             throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
  749.         }
  750.         low  = (Bval.cmp( Mval ) < 0);
  751.         high = (Bval.add( Mval ).cmp( tenSval ) > 0 );
  752.         digits[ndigit++] = (char)('0' + q);
  753.         }
  754.         if ( high && low ){
  755.         Bval.lshiftMe(1);
  756.         lowDigitDifference = Bval.cmp(tenSval);
  757.         } else
  758.         lowDigitDifference = 0L; // this here only for flow analysis!
  759.     }
  760.     this.decExponent = decExp+1;
  761.     this.digits = digits;
  762.     this.nDigits = ndigit;
  763.     /*
  764.      * Last digit gets rounded based on stopping condition.
  765.      */
  766.     if ( high ){
  767.         if ( low ){
  768.         if ( lowDigitDifference == 0L ){
  769.             // it's a tie!
  770.             // choose based on which digits we like.
  771.             if ( (digits[nDigits-1]&1) != 0 ) roundup();
  772.         } else if ( lowDigitDifference > 0 ){
  773.             roundup();
  774.         }
  775.         } else {
  776.         roundup();
  777.         }
  778.     }
  779.     }
  780.  
  781.     public String
  782.     toString(){
  783.     // most brain-dead version
  784.     StringBuffer result = new StringBuffer( nDigits+8 );
  785.     if ( isNegative ){ result.append( '-' ); }
  786.     if ( isExceptional ){
  787.         result.append( digits, 0, nDigits );
  788.     } else {
  789.         result.append( "0.");
  790.         result.append( digits, 0, nDigits );
  791.         result.append('e');
  792.         result.append( decExponent );
  793.     }
  794.     return new String(result);
  795.     }
  796.  
  797.     public String
  798.     toJavaFormatString(){
  799.     char result[] = new char[ nDigits + 10 ];
  800.     int  i = 0;
  801.     if ( isNegative ){ result[0] = '-'; i = 1; }
  802.     if ( isExceptional ){
  803.         System.arraycopy( digits, 0, result, i, nDigits );
  804.         i += nDigits;
  805.     } else {
  806.         if ( decExponent > 0 && decExponent < 8 ){
  807.         // print digits.digits.
  808.         int charLength = Math.min( nDigits, decExponent );
  809.         System.arraycopy( digits, 0, result, i, charLength );
  810.         i += charLength;
  811.         if ( charLength < decExponent ){
  812.             charLength = decExponent-charLength;
  813.             System.arraycopy( zero, 0, result, i, charLength );
  814.             i += charLength;
  815.             result[i++] = '.';
  816.             result[i++] = '0';
  817.         } else {
  818.             result[i++] = '.';
  819.             if ( charLength < nDigits ){
  820.             int t = nDigits - charLength;
  821.             System.arraycopy( digits, charLength, result, i, t );
  822.             i += t;
  823.             } else{
  824.             result[i++] = '0';
  825.             }
  826.         }
  827.         } else if ( decExponent <=0 && decExponent > -3 ){
  828.         result[i++] = '0';
  829.         result[i++] = '.';
  830.         if ( decExponent != 0 ){
  831.             System.arraycopy( zero, 0, result, i, -decExponent );
  832.             i -= decExponent;
  833.         }
  834.         System.arraycopy( digits, 0, result, i, nDigits );
  835.         i += nDigits;
  836.         } else {
  837.         result[i++] = digits[0];
  838.         result[i++] = '.';
  839.         if ( nDigits > 1 ){
  840.             System.arraycopy( digits, 1, result, i, nDigits-1 );
  841.             i += nDigits-1;
  842.         } else {
  843.             result[i++] = '0';
  844.         }
  845.         result[i++] = 'E';
  846.         int e;
  847.         if ( decExponent <= 0 ){
  848.             result[i++] = '-';
  849.             e = -decExponent+1;
  850.         } else {
  851.             e = decExponent-1;
  852.         }
  853.         // decExponent has 1, 2, or 3, digits
  854.         if ( e <= 9 ) {
  855.             result[i++] = (char)( e+'0' );
  856.         } else if ( e <= 99 ){
  857.             result[i++] = (char)( e/10 +'0' );
  858.             result[i++] = (char)( e%10 + '0' );
  859.         } else {
  860.             result[i++] = (char)(e/100+'0');
  861.             e %= 100;
  862.             result[i++] = (char)(e/10+'0');
  863.             result[i++] = (char)( e%10 + '0' );
  864.         }
  865.         }
  866.     }
  867.     return new String(result, 0, i);
  868.     }
  869.  
  870.     private static final int small5pow[] = {
  871.     1,
  872.     5,
  873.     5*5,
  874.     5*5*5,
  875.     5*5*5*5,
  876.     5*5*5*5*5,
  877.     5*5*5*5*5*5,
  878.     5*5*5*5*5*5*5,
  879.     5*5*5*5*5*5*5*5,
  880.     5*5*5*5*5*5*5*5*5,
  881.     5*5*5*5*5*5*5*5*5*5,
  882.     5*5*5*5*5*5*5*5*5*5*5,
  883.     5*5*5*5*5*5*5*5*5*5*5*5,
  884.     5*5*5*5*5*5*5*5*5*5*5*5*5
  885.     };
  886.  
  887.     private static final long long5pow[] = {
  888.     1L,
  889.     5L,
  890.     5L*5,
  891.     5L*5*5,
  892.     5L*5*5*5,
  893.     5L*5*5*5*5,
  894.     5L*5*5*5*5*5,
  895.     5L*5*5*5*5*5*5,
  896.     5L*5*5*5*5*5*5*5,
  897.     5L*5*5*5*5*5*5*5*5,
  898.     5L*5*5*5*5*5*5*5*5*5,
  899.     5L*5*5*5*5*5*5*5*5*5*5,
  900.     5L*5*5*5*5*5*5*5*5*5*5*5,
  901.     5L*5*5*5*5*5*5*5*5*5*5*5*5,
  902.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5,
  903.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  904.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  905.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  906.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  907.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  908.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  909.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  910.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  911.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  912.     5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
  913.     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,
  914.     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,
  915.     };
  916.  
  917.     // approximately ceil( log2( long5pow[i] ) )
  918.     private static final int n5bits[] = {
  919.     0,
  920.     3,
  921.     5,
  922.     7,
  923.     10,
  924.     12,
  925.     14,
  926.     17,
  927.     19,
  928.     21,
  929.     24,
  930.     26,
  931.     28,
  932.     31,
  933.     33,
  934.     35,
  935.     38,
  936.     40,
  937.     42,
  938.     45,
  939.     47,
  940.     49,
  941.     52,
  942.     54,
  943.     56,
  944.     59,
  945.     61,
  946.     };
  947.  
  948.     private static final char infinity[] = { 'I', 'n', 'f', 'i', 'n', 'i', 't', 'y' };
  949.     private static final char notANumber[] = { 'N', 'a', 'N' };
  950.     private static final char zero[] = { '0', '0', '0', '0', '0', '0', '0', '0' };
  951. }
  952.  
  953. /*
  954.  * A really, really simple bigint package
  955.  * tailored to the needs of floating base conversion.
  956.  */
  957. class FDBigInt {
  958.     int    nWords; // number of words used
  959.     int data[]; // value: data[0] is least significant
  960.  
  961.     private static boolean debugging = false;
  962.  
  963.     public static void setDebugging( boolean d ) { debugging = d; }
  964.  
  965.     public FDBigInt( int v ){
  966.     nWords = 1;
  967.     data = new int[1];
  968.     data[0] = v;
  969.     }
  970.  
  971.     public FDBigInt( long v ){
  972.     data = new int[2];
  973.     data[0] = (int)v;
  974.     data[1] = (int)(v>>>32);
  975.     nWords = (data[1]==0) ? 1 : 2;
  976.     }
  977.  
  978.     public FDBigInt( FDBigInt other ){
  979.     data = new int[nWords = other.nWords];
  980.     System.arraycopy( other.data, 0, data, 0, nWords );
  981.     }
  982.  
  983.     private FDBigInt( int [] d, int n ){
  984.     data = d;
  985.     nWords = n;
  986.     }
  987.  
  988.     /*
  989.      * Left shift by c bits.
  990.      * Shifts this in place.
  991.      */
  992.     public void
  993.     lshiftMe( int c )throws IllegalArgumentException {
  994.     if ( c <= 0 ){
  995.         if ( c == 0 )
  996.         return; // silly.
  997.         else
  998.         throw new IllegalArgumentException("negative shift count");
  999.     }
  1000.     int wordcount = c>>5;
  1001.     int bitcount  = c & 0x1f;
  1002.     int anticount = 32-bitcount;
  1003.     int t[] = data;
  1004.     int s[] = data;
  1005.     if ( nWords+wordcount+1 > t.length ){
  1006.         // reallocate.
  1007.         t = new int[ nWords+wordcount+1 ];
  1008.     }
  1009.     int target = nWords+wordcount;
  1010.     int src    = nWords-1;
  1011.     if ( bitcount == 0 ){
  1012.         // special hack, since an anticount of 32 won't go!
  1013.         System.arraycopy( s, 0, t, wordcount, nWords );
  1014.         target = wordcount-1;
  1015.     } else {
  1016.         t[target--] = s[src]>>>anticount;
  1017.         while ( src >= 1 ){
  1018.         t[target--] = (s[src]<<bitcount) | (s[--src]>>>anticount);
  1019.         }
  1020.         t[target--] = s[src]<<bitcount;
  1021.     }
  1022.     while( target >= 0 ){
  1023.         t[target--] = 0;
  1024.     }
  1025.     data = t;
  1026.     nWords += wordcount + 1;
  1027.     // may have constructed high-order word of 0.
  1028.     // if so, trim it
  1029.     while ( nWords > 1 && data[nWords-1] == 0 )
  1030.         nWords--;
  1031.     }
  1032.  
  1033.     /*
  1034.      * normalize this number by shifting until
  1035.      * the MSB of the number is at 0x08000000.
  1036.      * This is in preparation for quoRemIteration, below.
  1037.      * The idea is that, to make division easier, we want the
  1038.      * divisor to be "normalized" -- usually this means shifting
  1039.      * the MSB into the high words sign bit. But because we know that
  1040.      * the quotient will be 0 < q < 10, we would like to arrange that
  1041.      * the dividend not span up into another word of precision.
  1042.      * (This needs to be explained more clearly!)
  1043.      */
  1044.     public int
  1045.     normalizeMe() throws IllegalArgumentException {
  1046.     int src;
  1047.     int wordcount = 0;
  1048.     int bitcount  = 0;
  1049.     int v = 0;
  1050.     for ( src= nWords-1 ; src >= 0 && (v=data[src]) == 0 ; src--){
  1051.         wordcount += 1;
  1052.     }
  1053.     if ( src < 0 ){
  1054.         // oops. Value is zero. Cannot normalize it!
  1055.         throw new IllegalArgumentException("zero value");
  1056.     }
  1057.     /*
  1058.      * In most cases, we assume that wordcount is zero. This only
  1059.      * makes sense, as we try not to maintain any high-order
  1060.      * words full of zeros. In fact, if there are zeros, we will
  1061.      * simply SHORTEN our number at this point. Watch closely...
  1062.      */
  1063.     nWords -= wordcount;
  1064.     /*
  1065.      * Compute how far left we have to shift v s.t. its highest-
  1066.      * order bit is in the right place. Then call lshiftMe to
  1067.      * do the work.
  1068.      */
  1069.     if ( (v & 0xf0000000) != 0 ){
  1070.         // will have to shift up into the next word.
  1071.         // too bad.
  1072.         for( bitcount = 32 ; (v & 0xf0000000) != 0 ; bitcount-- )
  1073.         v >>>= 1;
  1074.     } else {
  1075.         while ( v <= 0x000fffff ){
  1076.         // hack: byte-at-a-time shifting
  1077.         v <<= 8;
  1078.         bitcount += 8;
  1079.         }
  1080.         while ( v <= 0x07ffffff ){
  1081.         v <<= 1;
  1082.         bitcount += 1;
  1083.         }
  1084.     }
  1085.     if ( bitcount != 0 )
  1086.         lshiftMe( bitcount );
  1087.     return bitcount;
  1088.     }
  1089.  
  1090.     /*
  1091.      * Multiply a FDBigInt by an int.
  1092.      * Result is a new FDBigInt.
  1093.      */
  1094.     public FDBigInt
  1095.     mult( int iv ) {
  1096.     long v = iv;
  1097.     int r[];
  1098.     long p;
  1099.  
  1100.     // guess adequate size of r.
  1101.     r = new int[ ( v * ((long)data[nWords-1]&0xffffffffL) > 0xfffffffL ) ? nWords+1 : nWords ];
  1102.     p = 0L;
  1103.     for( int i=0; i < nWords; i++ ) {
  1104.         p += v * ((long)data[i]&0xffffffffL);
  1105.         r[i] = (int)p;
  1106.         p >>>= 32;
  1107.     }
  1108.     if ( p == 0L){
  1109.         return new FDBigInt( r, nWords );
  1110.     } else {
  1111.         r[nWords] = (int)p;
  1112.         return new FDBigInt( r, nWords+1 );
  1113.     }
  1114.     }
  1115.  
  1116.     /*
  1117.      * Multiply a FDBigInt by another FDBigInt.
  1118.      * Result is a new FDBigInt.
  1119.      */
  1120.     public FDBigInt
  1121.     mult( FDBigInt other ){
  1122.     // crudely guess adequate size for r
  1123.     int r[] = new int[ nWords + other.nWords ];
  1124.     int i;
  1125.     // I think I am promised zeros...
  1126.  
  1127.     for( i = 0; i < this.nWords; i++ ){
  1128.         long v = (long)this.data[i] & 0xffffffffL; // UNSIGNED CONVERSION
  1129.         long p = 0L;
  1130.         int j;
  1131.         for( j = 0; j < other.nWords; j++ ){
  1132.         p += ((long)r[i+j]&0xffffffffL) + v*((long)other.data[j]&0xffffffffL); // UNSIGNED CONVERSIONS ALL 'ROUND.
  1133.         r[i+j] = (int)p;
  1134.         p >>>= 32;
  1135.         }
  1136.         r[i+j] = (int)p;
  1137.     }
  1138.     // compute how much of r we actually needed for all that.
  1139.     for ( i = r.length-1; i> 0; i--)
  1140.         if ( r[i] != 0 )
  1141.         break;
  1142.     return new FDBigInt( r, i+1 );
  1143.     }
  1144.  
  1145.     /*
  1146.      * Add one FDBigInt to another. Return a FDBigInt
  1147.      */
  1148.     public FDBigInt
  1149.     add( FDBigInt other ){
  1150.     int i;
  1151.     int a[], b[];
  1152.     int n, m;
  1153.     long c = 0L;
  1154.     // arrange such that a.nWords >= b.nWords;
  1155.     // n = a.nWords, m = b.nWords
  1156.     if ( this.nWords >= other.nWords ){
  1157.         a = this.data;
  1158.         n = this.nWords;
  1159.         b = other.data;
  1160.         m = other.nWords;
  1161.     } else {
  1162.         a = other.data;
  1163.         n = other.nWords;
  1164.         b = this.data;
  1165.         m = this.nWords;
  1166.     }
  1167.     int r[] = new int[ n ];
  1168.     for ( i = 0; i < n; i++ ){
  1169.         c += (long)a[i] & 0xffffffffL;
  1170.         if ( i < m ){
  1171.         c += (long)b[i] & 0xffffffffL;
  1172.         }
  1173.         r[i] = (int) c;
  1174.         c >>= 32; // signed shift.
  1175.     }
  1176.     if ( c != 0L ){
  1177.         // oops -- carry out -- need longer result.
  1178.         int s[] = new int[ r.length+1 ];
  1179.         System.arraycopy( r, 0, s, 0, r.length );
  1180.         s[i++] = (int)c;
  1181.         return new FDBigInt( s, i );
  1182.     }
  1183.     return new FDBigInt( r, i );
  1184.     }
  1185.  
  1186.     /*
  1187.      * Subtract one FDBigInt from another. Return a FDBigInt
  1188.      * Assert that the result is positive.
  1189.      */
  1190.     public FDBigInt
  1191.     sub( FDBigInt other ){
  1192.     int r[] = new int[ this.nWords ];
  1193.     int i;
  1194.     int n = this.nWords;
  1195.     int m = other.nWords;
  1196.     int nzeros = 0;
  1197.     long c = 0L;
  1198.     for ( i = 0; i < n; i++ ){
  1199.         c += (long)this.data[i] & 0xffffffffL;
  1200.         if ( i < m ){
  1201.         c -= (long)other.data[i] & 0xffffffffL;
  1202.         }
  1203.         if ( ( r[i] = (int) c ) == 0 )
  1204.         nzeros++;
  1205.         else
  1206.         nzeros = 0;
  1207.         c >>= 32; // signed shift.
  1208.     }
  1209.     if ( c != 0L )
  1210.         throw new RuntimeException("Assertion botch: borrow out of subtract");
  1211.     while ( i < m )
  1212.         if ( other.data[i++] != 0 )
  1213.         throw new RuntimeException("Assertion botch: negative result of subtract");
  1214.     return new FDBigInt( r, n-nzeros );
  1215.     }
  1216.  
  1217.     /*
  1218.      * Compare FDBigInt with another FDBigInt. Return an integer
  1219.      * >0: this > other
  1220.      *  0: this == other
  1221.      * <0: this < other
  1222.      */
  1223.     public int
  1224.     cmp( FDBigInt other ){
  1225.         int i;
  1226.     if ( this.nWords > other.nWords ){
  1227.         // if any of my high-order words is non-zero,
  1228.         // then the answer is evident
  1229.         int j = other.nWords-1;
  1230.         for ( i = this.nWords-1; i > j ; i-- )
  1231.         if ( this.data[i] != 0 ) return 1;
  1232.     }else if ( this.nWords < other.nWords ){
  1233.         // if any of other's high-order words is non-zero,
  1234.         // then the answer is evident
  1235.         int j = this.nWords-1;
  1236.         for ( i = other.nWords-1; i > j ; i-- )
  1237.         if ( other.data[i] != 0 ) return -1;
  1238.     } else{
  1239.         i = this.nWords-1;
  1240.     }
  1241.     for ( ; i > 0 ; i-- )
  1242.         if ( this.data[i] != other.data[i] )
  1243.         break;
  1244.     // careful! want unsigned compare!
  1245.     // use brute force here.
  1246.     int a = this.data[i];
  1247.     int b = other.data[i];
  1248.     if ( a < 0 ){
  1249.         // a is really big, unsigned
  1250.         if ( b < 0 ){
  1251.         return a-b; // both big, negative
  1252.         } else {
  1253.         return 1; // b not big, answer is obvious;
  1254.         }
  1255.     } else {
  1256.         // a is not really big
  1257.         if ( b < 0 ) {
  1258.         // but b is really big
  1259.         return -1;
  1260.         } else {
  1261.         return a - b;
  1262.         }
  1263.     }
  1264.     }
  1265.  
  1266.     /*
  1267.      * Compute
  1268.      * q = (int)( this / S )
  1269.      * this = 10 * ( this mod S )
  1270.      * Return q.
  1271.      * This is the iteration step of digit development for output.
  1272.      * We assume that S has been normalized, as above, and that
  1273.      * "this" has been lshift'ed accordingly.
  1274.      * Also assume, of course, that the result, q, can be expressed
  1275.      * as an integer, 0 <= q < 10.
  1276.      */
  1277.     public int
  1278.     quoRemIteration( FDBigInt S )throws IllegalArgumentException {
  1279.     // ensure that this and S have the same number of
  1280.     // digits. If S is properly normalized and q < 10 then
  1281.     // this must be so.
  1282.     if ( nWords != S.nWords ){
  1283.         throw new IllegalArgumentException("disparate values");
  1284.     }
  1285.     // estimate q the obvious way. We will usually be
  1286.     // right. If not, then we're only off by a little and
  1287.     // will re-add.
  1288.     int n = nWords-1;
  1289.     long q = ((long)data[n]&0xffffffffL) / (long)S.data[n];
  1290.     long diff = 0L;
  1291.     for ( int i = 0; i <= n ; i++ ){
  1292.         diff += ((long)data[i]&0xffffffffL) -  q*((long)S.data[i]&0xffffffffL);
  1293.         data[i] = (int)diff;
  1294.         diff >>= 32; // N.B. SIGNED shift.
  1295.     }
  1296.     if ( diff != 0L ) {
  1297.         // damn, damn, damn. q is too big.
  1298.         // add S back in until this turns +. This should
  1299.         // not be very many times!
  1300.         long sum = 0L;
  1301.         while ( sum ==  0L ){
  1302.         sum = 0L;
  1303.         for ( int i = 0; i <= n; i++ ){
  1304.             sum += ((long)data[i]&0xffffffffL) +  ((long)S.data[i]&0xffffffffL);
  1305.             data[i] = (int) sum;
  1306.             sum >>= 32; // Signed or unsigned, answer is 0 or 1
  1307.         }
  1308.         /*
  1309.          * Originally the following line read
  1310.          * "if ( sum !=0 && sum != -1 )"
  1311.          * but that would be wrong, because of the 
  1312.          * treatment of the two values as entirely unsigned,
  1313.          * it would be impossible for a carry-out to be interpreted
  1314.          * as -1 -- it would have to be a single-bit carry-out, or
  1315.          * +1.
  1316.          */
  1317.         if ( sum !=0 && sum != 1 )
  1318.             throw new RuntimeException("Assertion botch: "+sum+" carry out of division correction");
  1319.         q -= 1;
  1320.         }
  1321.     }
  1322.     // finally, we can multiply this by 10.
  1323.     // it cannot overflow, right, as the high-order word has
  1324.     // at least 4 high-order zeros!
  1325.     long p = 0L;
  1326.     for ( int i = 0; i <= n; i++ ){
  1327.         p += 10*((long)data[i]&0xffffffffL);
  1328.         data[i] = (int)p;
  1329.         p >>= 32; // SIGNED shift.
  1330.     }
  1331.     if ( p != 0L )
  1332.         throw new RuntimeException("Assertion botch: carry out of *10");
  1333.  
  1334.     return (int)q;
  1335.     }
  1336.  
  1337.     public long
  1338.     longValue(){
  1339.     // if this can be represented as a long,
  1340.     // return the value
  1341.     int i;
  1342.     for ( i = this.nWords-1; i > 1 ; i-- ){
  1343.         if ( data[i] != 0 ){
  1344.         throw new RuntimeException("Assertion botch: value too big");
  1345.         }
  1346.     }
  1347.     switch(i){
  1348.     case 1:
  1349.         if ( data[1] < 0 )
  1350.         throw new RuntimeException("Assertion botch: value too big");
  1351.         return ((long)(data[1]) << 32) | ((long)data[0]&0xffffffffL);
  1352.     case 0:
  1353.         return ((long)data[0]&0xffffffffL);
  1354.     default:
  1355.         throw new RuntimeException("Assertion botch: longValue confused");
  1356.     }
  1357.     }
  1358.  
  1359.     public String
  1360.     toString() {
  1361.     StringBuffer r = new StringBuffer(30);
  1362.     r.append('[');
  1363.     int i = Math.min( nWords-1, data.length-1) ;
  1364.     if ( nWords > data.length ){
  1365.         r.append( "("+data.length+"<"+nWords+"!)" );
  1366.     }
  1367.     for( ; i> 0 ; i-- ){
  1368.         r.append( Integer.toHexString( data[i] ) );
  1369.         r.append( (char) ' ' );
  1370.     }
  1371.     r.append( Integer.toHexString( data[0] ) );
  1372.     r.append( (char) ']' );
  1373.     return new String( r );
  1374.     }
  1375. }
  1376.