FloatingDecimalpublic class FloatingDecimal extends Object
| Fields Summary |
|---|
| boolean | isExceptional | | boolean | isNegative | | int | decExponent | | char[] | digits | | int | nDigits | | int | bigIntExp | | int | bigIntNBits | | boolean | mustSetRoundDir | | int | roundDir | | static final long | signMask | | static final long | expMask | | static final long | fractMask | | static final int | expShift | | static final int | expBias | | static final long | fractHOB | | static final long | expOne | | static final int | maxSmallBinExp | | static final int | minSmallBinExp | | static final int | maxDecimalDigits | | static final int | maxDecimalExponent | | static final int | minDecimalExponent | | static final int | bigDecimalExponent | | static final long | highbyte | | static final long | highbit | | static final long | lowbytes | | static final int | singleSignMask | | static final int | singleExpMask | | static final int | singleFractMask | | static final int | singleExpShift | | static final int | singleFractHOB | | static final int | singleExpBias | | static final int | singleMaxDecimalDigits | | static final int | singleMaxDecimalExponent | | static final int | singleMinDecimalExponent | | static final int | intDecimalDigits | | private static FDBigInt[] | b5p | | private static final double[] | small10pow | | private static final float[] | singleSmall10pow | | private static final double[] | big10pow | | private static final double[] | tiny10pow | | private static final int | maxSmallTen | | private static final int | singleMaxSmallTen | | private static final int[] | small5pow | | private static final long[] | long5pow | | private static final int[] | n5bits | | private static final char[] | infinity | | private static final char[] | notANumber | | private static final char[] | zero |
| Constructors Summary |
|---|
private FloatingDecimal(boolean negSign, int decExponent, char[] digits, int n, boolean e) // set by doubleValue
isNegative = negSign;
isExceptional = e;
this.decExponent = decExponent;
this.digits = digits;
this.nDigits = n;
| public FloatingDecimal(double d)
long dBits = Double.doubleToLongBits( d );
long fractBits;
int binExp;
int nSignificantBits;
// discover and delete sign
if ( (dBits&signMask) != 0 ){
isNegative = true;
dBits ^= signMask;
} else {
isNegative = false;
}
// Begin to unpack
// Discover obvious special cases of NaN and Infinity.
binExp = (int)( (dBits&expMask) >> expShift );
fractBits = dBits&fractMask;
if ( binExp == (int)(expMask>>expShift) ) {
isExceptional = true;
if ( fractBits == 0L ){
digits = infinity;
} else {
digits = notANumber;
isNegative = false; // NaN has no sign!
}
nDigits = digits.length;
return;
}
isExceptional = false;
// Finish unpacking
// Normalize denormalized numbers.
// Insert assumed high-order bit for normalized numbers.
// Subtract exponent bias.
if ( binExp == 0 ){
if ( fractBits == 0L ){
// not a denorm, just a 0!
decExponent = 0;
digits = zero;
nDigits = 1;
return;
}
while ( (fractBits&fractHOB) == 0L ){
fractBits <<= 1;
binExp -= 1;
}
nSignificantBits = expShift + binExp +1; // recall binExp is - shift count.
binExp += 1;
} else {
fractBits |= fractHOB;
nSignificantBits = expShift+1;
}
binExp -= expBias;
// call the routine that actually does all the hard work.
dtoa( binExp, fractBits, nSignificantBits );
| public FloatingDecimal(float f)
int fBits = Float.floatToIntBits( f );
int fractBits;
int binExp;
int nSignificantBits;
// discover and delete sign
if ( (fBits&singleSignMask) != 0 ){
isNegative = true;
fBits ^= singleSignMask;
} else {
isNegative = false;
}
// Begin to unpack
// Discover obvious special cases of NaN and Infinity.
binExp = (int)( (fBits&singleExpMask) >> singleExpShift );
fractBits = fBits&singleFractMask;
if ( binExp == (int)(singleExpMask>>singleExpShift) ) {
isExceptional = true;
if ( fractBits == 0L ){
digits = infinity;
} else {
digits = notANumber;
isNegative = false; // NaN has no sign!
}
nDigits = digits.length;
return;
}
isExceptional = false;
// Finish unpacking
// Normalize denormalized numbers.
// Insert assumed high-order bit for normalized numbers.
// Subtract exponent bias.
if ( binExp == 0 ){
if ( fractBits == 0 ){
// not a denorm, just a 0!
decExponent = 0;
digits = zero;
nDigits = 1;
return;
}
while ( (fractBits&singleFractHOB) == 0 ){
fractBits <<= 1;
binExp -= 1;
}
nSignificantBits = singleExpShift + binExp +1; // recall binExp is - shift count.
binExp += 1;
} else {
fractBits |= singleFractHOB;
nSignificantBits = singleExpShift+1;
}
binExp -= singleExpBias;
// call the routine that actually does all the hard work.
dtoa( binExp, ((long)fractBits)<<(expShift-singleExpShift), nSignificantBits );
|
| Methods Summary |
|---|
| private static synchronized java.lang.FDBigInt | big5pow(int p)
if ( p < 0 )
throw new RuntimeException( "Assertion botch: negative power of 5");
if ( b5p == null ){
b5p = new FDBigInt[ p+1 ];
}else if (b5p.length <= p ){
FDBigInt t[] = new FDBigInt[ p+1 ];
System.arraycopy( b5p, 0, t, 0, b5p.length );
b5p = t;
}
if ( b5p[p] != null )
return b5p[p];
else if ( p < small5pow.length )
return b5p[p] = new FDBigInt( small5pow[p] );
else if ( p < long5pow.length )
return b5p[p] = new FDBigInt( long5pow[p] );
else {
// construct the value.
// recursively.
int q, r;
// in order to compute 5^p,
// compute its square root, 5^(p/2) and square.
// or, let q = p / 2, r = p -q, then
// 5^p = 5^(q+r) = 5^q * 5^r
q = p >> 1;
r = p - q;
FDBigInt bigq = b5p[q];
if ( bigq == null )
bigq = big5pow ( q );
if ( r < small5pow.length ){
return (b5p[p] = bigq.mult( small5pow[r] ) );
}else{
FDBigInt bigr = b5p[ r ];
if ( bigr == null )
bigr = big5pow( r );
return (b5p[p] = bigq.mult( bigr ) );
}
}
| | private static java.lang.FDBigInt | constructPow52(int p5, int p2)
FDBigInt v = new FDBigInt( big5pow( p5 ) );
if ( p2 != 0 ){
v.lshiftMe( p2 );
}
return v;
| | private static int | countBits(long v)
/*
* count number of bits from high-order 1 bit to low-order 1 bit,
* inclusive.
*/
//
// the strategy is to shift until we get a non-zero sign bit
// then shift until we have no bits left, counting the difference.
// we do byte shifting as a hack. Hope it helps.
//
if ( v == 0L ) return 0;
while ( ( v & highbyte ) == 0L ){
v <<= 8;
}
while ( v > 0L ) { // i.e. while ((v&highbit) == 0L )
v <<= 1;
}
int n = 0;
while (( v & lowbytes ) != 0L ){
v <<= 8;
n += 8;
}
while ( v != 0L ){
v <<= 1;
n += 1;
}
return n;
| | private void | developLongDigits(int decExponent, long lvalue, long insignificant)
char digits[];
int ndigits;
int digitno;
int c;
//
// Discard non-significant low-order bits, while rounding,
// up to insignificant value.
int i;
for ( i = 0; insignificant >= 10L; i++ )
insignificant /= 10L;
if ( i != 0 ){
long pow10 = long5pow[i] << i; // 10^i == 5^i * 2^i;
long residue = lvalue % pow10;
lvalue /= pow10;
decExponent += i;
if ( residue >= (pow10>>1) ){
// round up based on the low-order bits we're discarding
lvalue++;
}
}
if ( lvalue <= Integer.MAX_VALUE ){
if ( lvalue <= 0L )
throw new RuntimeException("Assertion botch: value "+lvalue+" <= 0");
// even easier subcase!
// can do int arithmetic rather than long!
int ivalue = (int)lvalue;
digits = new char[ ndigits=10 ];
digitno = ndigits-1;
c = ivalue%10;
ivalue /= 10;
while ( c == 0 ){
decExponent++;
c = ivalue%10;
ivalue /= 10;
}
while ( ivalue != 0){
digits[digitno--] = (char)(c+'0");
decExponent++;
c = ivalue%10;
ivalue /= 10;
}
digits[digitno] = (char)(c+'0");
} else {
// same algorithm as above (same bugs, too )
// but using long arithmetic.
digits = new char[ ndigits=20 ];
digitno = ndigits-1;
c = (int)(lvalue%10L);
lvalue /= 10L;
while ( c == 0 ){
decExponent++;
c = (int)(lvalue%10L);
lvalue /= 10L;
}
while ( lvalue != 0L ){
digits[digitno--] = (char)(c+'0");
decExponent++;
c = (int)(lvalue%10L);
lvalue /= 10;
}
digits[digitno] = (char)(c+'0");
}
char result [];
ndigits -= digitno;
if ( digitno == 0 )
result = digits;
else {
result = new char[ ndigits ];
System.arraycopy( digits, digitno, result, 0, ndigits );
}
this.digits = result;
this.decExponent = decExponent+1;
this.nDigits = ndigits;
| | private java.lang.FDBigInt | doubleToBigInt(double dval)
long lbits = Double.doubleToLongBits( dval ) & ~signMask;
int binexp = (int)(lbits >>> expShift);
lbits &= fractMask;
if ( binexp > 0 ){
lbits |= fractHOB;
} else {
if ( lbits == 0L )
throw new RuntimeException("Assertion botch: doubleToBigInt(0.0)");
binexp +=1;
while ( (lbits & fractHOB ) == 0L){
lbits <<= 1;
binexp -= 1;
}
}
binexp -= expBias;
int nbits = countBits( lbits );
/*
* We now know where the high-order 1 bit is,
* and we know how many there are.
*/
int lowOrderZeros = expShift+1-nbits;
lbits >>>= lowOrderZeros;
bigIntExp = binexp+1-nbits;
bigIntNBits = nbits;
return new FDBigInt( lbits );
| | public double | doubleValue()
int kDigits = Math.min( nDigits, maxDecimalDigits+1 );
long lValue;
double dValue;
double rValue, tValue;
roundDir = 0;
/*
* convert the lead kDigits to a long integer.
*/
// (special performance hack: start to do it using int)
int iValue = (int)digits[0]-(int)'0";
int iDigits = Math.min( kDigits, intDecimalDigits );
for ( int i=1; i < iDigits; i++ ){
iValue = iValue*10 + (int)digits[i]-(int)'0";
}
lValue = (long)iValue;
for ( int i=iDigits; i < kDigits; i++ ){
lValue = lValue*10L + (long)((int)digits[i]-(int)'0");
}
dValue = (double)lValue;
int exp = decExponent-kDigits;
/*
* lValue now contains a long integer with the value of
* the first kDigits digits of the number.
* dValue contains the (double) of the same.
*/
if ( nDigits <= maxDecimalDigits ){
/*
* possibly an easy case.
* We know that the digits can be represented
* exactly. And if the exponent isn't too outrageous,
* the whole thing can be done with one operation,
* thus one rounding error.
* Note that all our constructors trim all leading and
* trailing zeros, so simple values (including zero)
* will always end up here
*/
if (exp == 0 || dValue == 0.0)
return (isNegative)? -dValue : dValue; // small floating integer
else if ( exp >= 0 ){
if ( exp <= maxSmallTen ){
/*
* Can get the answer with one operation,
* thus one roundoff.
*/
rValue = dValue * small10pow[exp];
if ( mustSetRoundDir ){
tValue = rValue / small10pow[exp];
roundDir = ( tValue == dValue ) ? 0
:( tValue < dValue ) ? 1
: -1;
}
return (isNegative)? -rValue : rValue;
}
int slop = maxDecimalDigits - kDigits;
if ( exp <= maxSmallTen+slop ){
/*
* We can multiply dValue by 10^(slop)
* and it is still "small" and exact.
* Then we can multiply by 10^(exp-slop)
* with one rounding.
*/
dValue *= small10pow[slop];
rValue = dValue * small10pow[exp-slop];
if ( mustSetRoundDir ){
tValue = rValue / small10pow[exp-slop];
roundDir = ( tValue == dValue ) ? 0
:( tValue < dValue ) ? 1
: -1;
}
return (isNegative)? -rValue : rValue;
}
/*
* Else we have a hard case with a positive exp.
*/
} else {
if ( exp >= -maxSmallTen ){
/*
* Can get the answer in one division.
*/
rValue = dValue / small10pow[-exp];
tValue = rValue * small10pow[-exp];
if ( mustSetRoundDir ){
roundDir = ( tValue == dValue ) ? 0
:( tValue < dValue ) ? 1
: -1;
}
return (isNegative)? -rValue : rValue;
}
/*
* Else we have a hard case with a negative exp.
*/
}
}
/*
* Harder cases:
* The sum of digits plus exponent is greater than
* what we think we can do with one error.
*
* Start by approximating the right answer by,
* naively, scaling by powers of 10.
*/
if ( exp > 0 ){
if ( decExponent > maxDecimalExponent+1 ){
/*
* Lets face it. This is going to be
* Infinity. Cut to the chase.
*/
return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
}
if ( (exp&15) != 0 ){
dValue *= small10pow[exp&15];
}
if ( (exp>>=4) != 0 ){
int j;
for( j = 0; exp > 1; j++, exp>>=1 ){
if ( (exp&1)!=0)
dValue *= big10pow[j];
}
/*
* The reason for the weird exp > 1 condition
* in the above loop was so that the last multiply
* would get unrolled. We handle it here.
* It could overflow.
*/
double t = dValue * big10pow[j];
if ( Double.isInfinite( t ) ){
/*
* It did overflow.
* Look more closely at the result.
* If the exponent is just one too large,
* then use the maximum finite as our estimate
* value. Else call the result infinity
* and punt it.
* ( I presume this could happen because
* rounding forces the result here to be
* an ULP or two larger than
* Double.MAX_VALUE ).
*/
t = dValue / 2.0;
t *= big10pow[j];
if ( Double.isInfinite( t ) ){
return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
}
t = Double.MAX_VALUE;
}
dValue = t;
}
} else if ( exp < 0 ){
exp = -exp;
if ( decExponent < minDecimalExponent-1 ){
/*
* Lets face it. This is going to be
* zero. Cut to the chase.
*/
return (isNegative)? -0.0 : 0.0;
}
if ( (exp&15) != 0 ){
dValue /= small10pow[exp&15];
}
if ( (exp>>=4) != 0 ){
int j;
for( j = 0; exp > 1; j++, exp>>=1 ){
if ( (exp&1)!=0)
dValue *= tiny10pow[j];
}
/*
* The reason for the weird exp > 1 condition
* in the above loop was so that the last multiply
* would get unrolled. We handle it here.
* It could underflow.
*/
double t = dValue * tiny10pow[j];
if ( t == 0.0 ){
/*
* It did underflow.
* Look more closely at the result.
* If the exponent is just one too small,
* then use the minimum finite as our estimate
* value. Else call the result 0.0
* and punt it.
* ( I presume this could happen because
* rounding forces the result here to be
* an ULP or two less than
* Double.MIN_VALUE ).
*/
t = dValue * 2.0;
t *= tiny10pow[j];
if ( t == 0.0 ){
return (isNegative)? -0.0 : 0.0;
}
t = Double.MIN_VALUE;
}
dValue = t;
}
}
/*
* dValue is now approximately the result.
* The hard part is adjusting it, by comparison
* with FDBigInt arithmetic.
* Formulate the EXACT big-number result as
* bigD0 * 10^exp
*/
FDBigInt bigD0 = new FDBigInt( lValue, digits, kDigits, nDigits );
exp = decExponent - nDigits;
correctionLoop:
while(true){
/* AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES
* bigIntExp and bigIntNBits
*/
FDBigInt bigB = doubleToBigInt( dValue );
/*
* Scale bigD, bigB appropriately for
* big-integer operations.
* Naively, we multiply by powers of ten
* and powers of two. What we actually do
* is keep track of the powers of 5 and
* powers of 2 we would use, then factor out
* common divisors before doing the work.
*/
int B2, B5; // powers of 2, 5 in bigB
int D2, D5; // powers of 2, 5 in bigD
int Ulp2; // powers of 2 in halfUlp.
if ( exp >= 0 ){
B2 = B5 = 0;
D2 = D5 = exp;
} else {
B2 = B5 = -exp;
D2 = D5 = 0;
}
if ( bigIntExp >= 0 ){
B2 += bigIntExp;
} else {
D2 -= bigIntExp;
}
Ulp2 = B2;
// shift bigB and bigD left by a number s. t.
// halfUlp is still an integer.
int hulpbias;
if ( bigIntExp+bigIntNBits <= -expBias+1 ){
// This is going to be a denormalized number
// (if not actually zero).
// half an ULP is at 2^-(expBias+expShift+1)
hulpbias = bigIntExp+ expBias + expShift;
} else {
hulpbias = expShift + 2 - bigIntNBits;
}
B2 += hulpbias;
D2 += hulpbias;
// if there are common factors of 2, we might just as well
// factor them out, as they add nothing useful.
int common2 = Math.min( B2, Math.min( D2, Ulp2 ) );
B2 -= common2;
D2 -= common2;
Ulp2 -= common2;
// do multiplications by powers of 5 and 2
bigB = multPow52( bigB, B5, B2 );
FDBigInt bigD = multPow52( new FDBigInt( bigD0 ), D5, D2 );
//
// to recap:
// bigB is the scaled-big-int version of our floating-point
// candidate.
// bigD is the scaled-big-int version of the exact value
// as we understand it.
// halfUlp is 1/2 an ulp of bigB, except for special cases
// of exact powers of 2
//
// the plan is to compare bigB with bigD, and if the difference
// is less than halfUlp, then we're satisfied. Otherwise,
// use the ratio of difference to halfUlp to calculate a fudge
// factor to add to the floating value, then go 'round again.
//
FDBigInt diff;
int cmpResult;
boolean overvalue;
if ( (cmpResult = bigB.cmp( bigD ) ) > 0 ){
overvalue = true; // our candidate is too big.
diff = bigB.sub( bigD );
if ( (bigIntNBits == 1) && (bigIntExp > -expBias) ){
// candidate is a normalized exact power of 2 and
// is too big. We will be subtracting.
// For our purposes, ulp is the ulp of the
// next smaller range.
Ulp2 -= 1;
if ( Ulp2 < 0 ){
// rats. Cannot de-scale ulp this far.
// must scale diff in other direction.
Ulp2 = 0;
diff.lshiftMe( 1 );
}
}
} else if ( cmpResult < 0 ){
overvalue = false; // our candidate is too small.
diff = bigD.sub( bigB );
} else {
// the candidate is exactly right!
// this happens with surprising frequency
break correctionLoop;
}
FDBigInt halfUlp = constructPow52( B5, Ulp2 );
if ( (cmpResult = diff.cmp( halfUlp ) ) < 0 ){
// difference is small.
// this is close enough
roundDir = overvalue ? -1 : 1;
break correctionLoop;
} else if ( cmpResult == 0 ){
// difference is exactly half an ULP
// round to some other value maybe, then finish
dValue += 0.5*ulp( dValue, overvalue );
// should check for bigIntNBits == 1 here??
roundDir = overvalue ? -1 : 1;
break correctionLoop;
} else {
// difference is non-trivial.
// could scale addend by ratio of difference to
// halfUlp here, if we bothered to compute that difference.
// Most of the time ( I hope ) it is about 1 anyway.
dValue += ulp( dValue, overvalue );
if ( dValue == 0.0 || dValue == Double.POSITIVE_INFINITY )
break correctionLoop; // oops. Fell off end of range.
continue; // try again.
}
}
return (isNegative)? -dValue : dValue;
| | private void | dtoa(int binExp, long fractBits, int nSignificantBits)
int nFractBits; // number of significant bits of fractBits;
int nTinyBits; // number of these to the right of the point.
int decExp;
// Examine number. Determine if it is an easy case,
// which we can do pretty trivially using float/long conversion,
// or whether we must do real work.
nFractBits = countBits( fractBits );
nTinyBits = Math.max( 0, nFractBits - binExp - 1 );
if ( binExp <= maxSmallBinExp && binExp >= minSmallBinExp ){
// Look more closely at the number to decide if,
// with scaling by 10^nTinyBits, the result will fit in
// a long.
if ( (nTinyBits < long5pow.length) && ((nFractBits + n5bits[nTinyBits]) < 64 ) ){
/*
* We can do this:
* take the fraction bits, which are normalized.
* (a) nTinyBits == 0: Shift left or right appropriately
* to align the binary point at the extreme right, i.e.
* where a long int point is expected to be. The integer
* result is easily converted to a string.
* (b) nTinyBits > 0: Shift right by expShift-nFractBits,
* which effectively converts to long and scales by
* 2^nTinyBits. Then multiply by 5^nTinyBits to
* complete the scaling. We know this won't overflow
* because we just counted the number of bits necessary
* in the result. The integer you get from this can
* then be converted to a string pretty easily.
*/
long halfULP;
if ( nTinyBits == 0 ) {
if ( binExp > nSignificantBits ){
halfULP = 1L << ( binExp-nSignificantBits-1);
} else {
halfULP = 0L;
}
if ( binExp >= expShift ){
fractBits <<= (binExp-expShift);
} else {
fractBits >>>= (expShift-binExp) ;
}
developLongDigits( 0, fractBits, halfULP );
return;
}
/*
* The following causes excess digits to be printed
* out in the single-float case. Our manipulation of
* halfULP here is apparently not correct. If we
* better understand how this works, perhaps we can
* use this special case again. But for the time being,
* we do not.
* else {
* fractBits >>>= expShift+1-nFractBits;
* fractBits *= long5pow[ nTinyBits ];
* halfULP = long5pow[ nTinyBits ] >> (1+nSignificantBits-nFractBits);
* developLongDigits( -nTinyBits, fractBits, halfULP );
* return;
* }
*/
}
}
/*
* This is the hard case. We are going to compute large positive
* integers B and S and integer decExp, s.t.
* d = ( B / S ) * 10^decExp
* 1 <= B / S < 10
* Obvious choices are:
* decExp = floor( log10(d) )
* B = d * 2^nTinyBits * 10^max( 0, -decExp )
* S = 10^max( 0, decExp) * 2^nTinyBits
* (noting that nTinyBits has already been forced to non-negative)
* I am also going to compute a large positive integer
* M = (1/2^nSignificantBits) * 2^nTinyBits * 10^max( 0, -decExp )
* i.e. M is (1/2) of the ULP of d, scaled like B.
* When we iterate through dividing B/S and picking off the
* quotient bits, we will know when to stop when the remainder
* is <= M.
*
* We keep track of powers of 2 and powers of 5.
*/
/*
* Estimate decimal exponent. (If it is small-ish,
* we could double-check.)
*
* First, scale the mantissa bits such that 1 <= d2 < 2.
* We are then going to estimate
* log10(d2) ~=~ (d2-1.5)/1.5 + log(1.5)
* and so we can estimate
* log10(d) ~=~ log10(d2) + binExp * log10(2)
* take the floor and call it decExp.
* FIXME -- use more precise constants here. It costs no more.
*/
double d2 = Double.longBitsToDouble(
expOne | ( fractBits &~ fractHOB ) );
decExp = (int)Math.floor(
(d2-1.5D)*0.289529654D + 0.176091259 + (double)binExp * 0.301029995663981 );
int B2, B5; // powers of 2 and powers of 5, respectively, in B
int S2, S5; // powers of 2 and powers of 5, respectively, in S
int M2, M5; // powers of 2 and powers of 5, respectively, in M
int Bbits; // binary digits needed to represent B, approx.
int tenSbits; // binary digits needed to represent 10*S, approx.
FDBigInt Sval, Bval, Mval;
B5 = Math.max( 0, -decExp );
B2 = B5 + nTinyBits + binExp;
S5 = Math.max( 0, decExp );
S2 = S5 + nTinyBits;
M5 = B5;
M2 = B2 - nSignificantBits;
/*
* the long integer fractBits contains the (nFractBits) interesting
* bits from the mantissa of d ( hidden 1 added if necessary) followed
* by (expShift+1-nFractBits) zeros. In the interest of compactness,
* I will shift out those zeros before turning fractBits into a
* FDBigInt. The resulting whole number will be
* d * 2^(nFractBits-1-binExp).
*/
fractBits >>>= (expShift+1-nFractBits);
B2 -= nFractBits-1;
int common2factor = Math.min( B2, S2 );
B2 -= common2factor;
S2 -= common2factor;
M2 -= common2factor;
/*
* HACK!! For exact powers of two, the next smallest number
* is only half as far away as we think (because the meaning of
* ULP changes at power-of-two bounds) for this reason, we
* hack M2. Hope this works.
*/
if ( nFractBits == 1 )
M2 -= 1;
if ( M2 < 0 ){
// oops.
// since we cannot scale M down far enough,
// we must scale the other values up.
B2 -= M2;
S2 -= M2;
M2 = 0;
}
/*
* Construct, Scale, iterate.
* Some day, we'll write a stopping test that takes
* account of the assymetry of the spacing of floating-point
* numbers below perfect powers of 2
* 26 Sept 96 is not that day.
* So we use a symmetric test.
*/
char digits[] = this.digits = new char[18];
int ndigit = 0;
boolean low, high;
long lowDigitDifference;
int q;
/*
* Detect the special cases where all the numbers we are about
* to compute will fit in int or long integers.
* In these cases, we will avoid doing FDBigInt arithmetic.
* We use the same algorithms, except that we "normalize"
* our FDBigInts before iterating. This is to make division easier,
* as it makes our fist guess (quotient of high-order words)
* more accurate!
*
* Some day, we'll write a stopping test that takes
* account of the assymetry of the spacing of floating-point
* numbers below perfect powers of 2
* 26 Sept 96 is not that day.
* So we use a symmetric test.
*/
Bbits = nFractBits + B2 + (( B5 < n5bits.length )? n5bits[B5] : ( B5*3 ));
tenSbits = S2+1 + (( (S5+1) < n5bits.length )? n5bits[(S5+1)] : ( (S5+1)*3 ));
if ( Bbits < 64 && tenSbits < 64){
if ( Bbits < 32 && tenSbits < 32){
// wa-hoo! They're all ints!
int b = ((int)fractBits * small5pow[B5] ) << B2;
int s = small5pow[S5] << S2;
int m = small5pow[M5] << M2;
int tens = s * 10;
/*
* Unroll the first iteration. If our decExp estimate
* was too high, our first quotient will be zero. In this
* case, we discard it and decrement decExp.
*/
ndigit = 0;
q = (int) ( b / s );
b = 10 * ( b % s );
m *= 10;
low = (b < m );
high = (b+m > tens );
if ( q >= 10 ){
// bummer, dude
throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
} else if ( (q == 0) && ! high ){
// oops. Usually ignore leading zero.
decExp--;
} else {
digits[ndigit++] = (char)('0" + q);
}
/*
* HACK! Java spec sez that we always have at least
* one digit after the . in either F- or E-form output.
* Thus we will need more than one digit if we're using
* E-form
*/
if ( decExp <= -3 || decExp >= 8 ){
high = low = false;
}
while( ! low && ! high ){
q = (int) ( b / s );
b = 10 * ( b % s );
m *= 10;
if ( q >= 10 ){
// bummer, dude
throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
}
if ( m > 0L ){
low = (b < m );
high = (b+m > tens );
} else {
// hack -- m might overflow!
// in this case, it is certainly > b,
// which won't
// and b+m > tens, too, since that has overflowed
// either!
low = true;
high = true;
}
digits[ndigit++] = (char)('0" + q);
}
lowDigitDifference = (b<<1) - tens;
} else {
// still good! they're all longs!
long b = (fractBits * long5pow[B5] ) << B2;
long s = long5pow[S5] << S2;
long m = long5pow[M5] << M2;
long tens = s * 10L;
/*
* Unroll the first iteration. If our decExp estimate
* was too high, our first quotient will be zero. In this
* case, we discard it and decrement decExp.
*/
ndigit = 0;
q = (int) ( b / s );
b = 10L * ( b % s );
m *= 10L;
low = (b < m );
high = (b+m > tens );
if ( q >= 10 ){
// bummer, dude
throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
} else if ( (q == 0) && ! high ){
// oops. Usually ignore leading zero.
decExp--;
} else {
digits[ndigit++] = (char)('0" + q);
}
/*
* HACK! Java spec sez that we always have at least
* one digit after the . in either F- or E-form output.
* Thus we will need more than one digit if we're using
* E-form
*/
if ( decExp <= -3 || decExp >= 8 ){
high = low = false;
}
while( ! low && ! high ){
q = (int) ( b / s );
b = 10 * ( b % s );
m *= 10;
if ( q >= 10 ){
// bummer, dude
throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
}
if ( m > 0L ){
low = (b < m );
high = (b+m > tens );
} else {
// hack -- m might overflow!
// in this case, it is certainly > b,
// which won't
// and b+m > tens, too, since that has overflowed
// either!
low = true;
high = true;
}
digits[ndigit++] = (char)('0" + q);
}
lowDigitDifference = (b<<1) - tens;
}
} else {
FDBigInt tenSval;
int shiftBias;
/*
* We really must do FDBigInt arithmetic.
* Fist, construct our FDBigInt initial values.
*/
Bval = multPow52( new FDBigInt( fractBits ), B5, B2 );
Sval = constructPow52( S5, S2 );
Mval = constructPow52( M5, M2 );
// normalize so that division works better
Bval.lshiftMe( shiftBias = Sval.normalizeMe() );
Mval.lshiftMe( shiftBias );
tenSval = Sval.mult( 10 );
/*
* Unroll the first iteration. If our decExp estimate
* was too high, our first quotient will be zero. In this
* case, we discard it and decrement decExp.
*/
ndigit = 0;
q = Bval.quoRemIteration( Sval );
Mval = Mval.mult( 10 );
low = (Bval.cmp( Mval ) < 0);
high = (Bval.add( Mval ).cmp( tenSval ) > 0 );
if ( q >= 10 ){
// bummer, dude
throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
} else if ( (q == 0) && ! high ){
// oops. Usually ignore leading zero.
decExp--;
} else {
digits[ndigit++] = (char)('0" + q);
}
/*
* HACK! Java spec sez that we always have at least
* one digit after the . in either F- or E-form output.
* Thus we will need more than one digit if we're using
* E-form
*/
if ( decExp <= -3 || decExp >= 8 ){
high = low = false;
}
while( ! low && ! high ){
q = Bval.quoRemIteration( Sval );
Mval = Mval.mult( 10 );
if ( q >= 10 ){
// bummer, dude
throw new RuntimeException( "Assertion botch: excessivly large digit "+q);
}
low = (Bval.cmp( Mval ) < 0);
high = (Bval.add( Mval ).cmp( tenSval ) > 0 );
digits[ndigit++] = (char)('0" + q);
}
if ( high && low ){
Bval.lshiftMe(1);
lowDigitDifference = Bval.cmp(tenSval);
} else
lowDigitDifference = 0L; // this here only for flow analysis!
}
this.decExponent = decExp+1;
this.digits = digits;
this.nDigits = ndigit;
/*
* Last digit gets rounded based on stopping condition.
*/
if ( high ){
if ( low ){
if ( lowDigitDifference == 0L ){
// it's a tie!
// choose based on which digits we like.
if ( (digits[nDigits-1]&1) != 0 ) roundup();
} else if ( lowDigitDifference > 0 ){
roundup();
}
} else {
roundup();
}
}
| | public float | floatValue()
int kDigits = Math.min( nDigits, singleMaxDecimalDigits+1 );
int iValue;
float fValue;
/*
* convert the lead kDigits to an integer.
*/
iValue = (int)digits[0]-(int)'0";
for ( int i=1; i < kDigits; i++ ){
iValue = iValue*10 + (int)digits[i]-(int)'0";
}
fValue = (float)iValue;
int exp = decExponent-kDigits;
/*
* iValue now contains an integer with the value of
* the first kDigits digits of the number.
* fValue contains the (float) of the same.
*/
if ( nDigits <= singleMaxDecimalDigits ){
/*
* possibly an easy case.
* We know that the digits can be represented
* exactly. And if the exponent isn't too outrageous,
* the whole thing can be done with one operation,
* thus one rounding error.
* Note that all our constructors trim all leading and
* trailing zeros, so simple values (including zero)
* will always end up here.
*/
if (exp == 0 || fValue == 0.0f)
return (isNegative)? -fValue : fValue; // small floating integer
else if ( exp >= 0 ){
if ( exp <= singleMaxSmallTen ){
/*
* Can get the answer with one operation,
* thus one roundoff.
*/
fValue *= singleSmall10pow[exp];
return (isNegative)? -fValue : fValue;
}
int slop = singleMaxDecimalDigits - kDigits;
if ( exp <= singleMaxSmallTen+slop ){
/*
* We can multiply dValue by 10^(slop)
* and it is still "small" and exact.
* Then we can multiply by 10^(exp-slop)
* with one rounding.
*/
fValue *= singleSmall10pow[slop];
fValue *= singleSmall10pow[exp-slop];
return (isNegative)? -fValue : fValue;
}
/*
* Else we have a hard case with a positive exp.
*/
} else {
if ( exp >= -singleMaxSmallTen ){
/*
* Can get the answer in one division.
*/
fValue /= singleSmall10pow[-exp];
return (isNegative)? -fValue : fValue;
}
/*
* Else we have a hard case with a negative exp.
*/
}
} else if ( (decExponent >= nDigits) && (nDigits+decExponent <= maxDecimalDigits) ){
/*
* In double-precision, this is an exact floating integer.
* So we can compute to double, then shorten to float
* with one round, and get the right answer.
*
* First, finish accumulating digits.
* Then convert that integer to a double, multiply
* by the appropriate power of ten, and convert to float.
*/
long lValue = (long)iValue;
for ( int i=kDigits; i < nDigits; i++ ){
lValue = lValue*10L + (long)((int)digits[i]-(int)'0");
}
double dValue = (double)lValue;
exp = decExponent-nDigits;
dValue *= small10pow[exp];
fValue = (float)dValue;
return (isNegative)? -fValue : fValue;
}
/*
* Harder cases:
* The sum of digits plus exponent is greater than
* what we think we can do with one error.
*
* Start by weeding out obviously out-of-range
* results, then convert to double and go to
* common hard-case code.
*/
if ( decExponent > singleMaxDecimalExponent+1 ){
/*
* Lets face it. This is going to be
* Infinity. Cut to the chase.
*/
return (isNegative)? Float.NEGATIVE_INFINITY : Float.POSITIVE_INFINITY;
} else if ( decExponent < singleMinDecimalExponent-1 ){
/*
* Lets face it. This is going to be
* zero. Cut to the chase.
*/
return (isNegative)? -0.0f : 0.0f;
}
/*
* Here, we do 'way too much work, but throwing away
* our partial results, and going and doing the whole
* thing as double, then throwing away half the bits that computes
* when we convert back to float.
*
* The alternative is to reproduce the whole multiple-precision
* algorithm for float precision, or to try to parameterize it
* for common usage. The former will take about 400 lines of code,
* and the latter I tried without success. Thus the semi-hack
* answer here.
*/
mustSetRoundDir = true;
double dValue = doubleValue();
return stickyRound( dValue );
| | private static java.lang.FDBigInt | multPow52(java.lang.FDBigInt v, int p5, int p2)
if ( p5 != 0 ){
if ( p5 < small5pow.length ){
v = v.mult( small5pow[p5] );
} else {
v = v.mult( big5pow( p5 ) );
}
}
if ( p2 != 0 ){
v.lshiftMe( p2 );
}
return v;
| | public static java.lang.FloatingDecimal | readJavaFormatString(java.lang.String in)
boolean isNegative = false;
boolean signSeen = false;
int decExp;
char c;
parseNumber:
try{
in = in.trim(); // don't fool around with white space.
// throws NullPointerException if null
int l = in.length();
if ( l == 0 ) throw new NumberFormatException("empty String");
int i = 0;
switch ( c = in.charAt( i ) ){
case '-":
isNegative = true;
//FALLTHROUGH
case '+":
i++;
signSeen = true;
}
// Would handle NaN and Infinity here, but it isn't
// part of the spec!
//
char[] digits = new char[ l ];
int nDigits= 0;
boolean decSeen = false;
int decPt = 0;
int nLeadZero = 0;
int nTrailZero= 0;
digitLoop:
while ( i < l ){
switch ( c = in.charAt( i ) ){
case '0":
if ( nDigits > 0 ){
nTrailZero += 1;
} else {
nLeadZero += 1;
}
break; // out of switch.
case '1":
case '2":
case '3":
case '4":
case '5":
case '6":
case '7":
case '8":
case '9":
while ( nTrailZero > 0 ){
digits[nDigits++] = '0";
nTrailZero -= 1;
}
digits[nDigits++] = c;
break; // out of switch.
case '.":
if ( decSeen ){
// already saw one ., this is the 2nd.
throw new NumberFormatException("multiple points");
}
decPt = i;
if ( signSeen ){
decPt -= 1;
}
decSeen = true;
break; // out of switch.
default:
break digitLoop;
}
i++;
}
/*
* At this point, we've scanned all the digits and decimal
* point we're going to see. Trim off leading and trailing
* zeros, which will just confuse us later, and adjust
* our initial decimal exponent accordingly.
* To review:
* we have seen i total characters.
* nLeadZero of them were zeros before any other digits.
* nTrailZero of them were zeros after any other digits.
* if ( decSeen ), then a . was seen after decPt characters
* ( including leading zeros which have been discarded )
* nDigits characters were neither lead nor trailing
* zeros, nor point
*/
/*
* special hack: if we saw no non-zero digits, then the
* answer is zero!
* Unfortunately, we feel honor-bound to keep parsing!
*/
if ( nDigits == 0 ){
digits = zero;
nDigits = 1;
if ( nLeadZero == 0 ){
// we saw NO DIGITS AT ALL,
// not even a crummy 0!
// this is not allowed.
break parseNumber; // go throw exception
}
}
/* Our initial exponent is decPt, adjusted by the number of
* discarded zeros. Or, if there was no decPt,
* then its just nDigits adjusted by discarded trailing zeros.
*/
if ( decSeen ){
decExp = decPt - nLeadZero;
} else {
decExp = nDigits+nTrailZero;
}
/*
* Look for 'e' or 'E' and an optionally signed integer.
*/
if ( (i < l) && ((c = in.charAt(i) )=='e") || (c == 'E") ){
int expSign = 1;
int expVal = 0;
int reallyBig = Integer.MAX_VALUE / 10;
boolean expOverflow = false;
switch( in.charAt(++i) ){
case '-":
expSign = -1;
//FALLTHROUGH
case '+":
i++;
}
int expAt = i;
expLoop:
while ( i < l ){
if ( expVal >= reallyBig ){
// the next character will cause integer
// overflow.
expOverflow = true;
}
switch ( c = in.charAt(i++) ){
case '0":
case '1":
case '2":
case '3":
case '4":
case '5":
case '6":
case '7":
case '8":
case '9":
expVal = expVal*10 + ( (int)c - (int)'0" );
continue;
default:
i--; // back up.
break expLoop; // stop parsing exponent.
}
}
int expLimit = bigDecimalExponent+nDigits+nTrailZero;
if ( expOverflow || ( expVal > expLimit ) ){
//
// The intent here is to end up with
// infinity or zero, as appropriate.
// The reason for yielding such a small decExponent,
// rather than something intuitive such as
// expSign*Integer.MAX_VALUE, is that this value
// is subject to further manipulation in
// doubleValue() and floatValue(), and I don't want
// it to be able to cause overflow there!
// (The only way we can get into trouble here is for
// really outrageous nDigits+nTrailZero, such as 2 billion. )
//
decExp = expSign*expLimit;
} else {
// this should not overflow, since we tested
// for expVal > (MAX+N), where N >= abs(decExp)
decExp = decExp + expSign*expVal;
}
// if we saw something not a digit ( or end of string )
// after the [Ee][+-], without seeing any digits at all
// this is certainly an error. If we saw some digits,
// but then some trailing garbage, that might be ok.
// so we just fall through in that case.
// HUMBUG
if ( i == expAt )
break parseNumber; // certainly bad
}
/*
* We parsed everything we could.
* If there are leftovers, then this is not good input!
*/
if ( i < l &&
((i != l - 1) ||
(in.charAt(i) != 'f" &&
in.charAt(i) != 'F" &&
in.charAt(i) != 'd" &&
in.charAt(i) != 'D"))) {
break parseNumber; // go throw exception
}
return new FloatingDecimal( isNegative, decExp, digits, nDigits, false );
} catch ( StringIndexOutOfBoundsException e ){ }
throw new NumberFormatException( in );
| | private void | roundup()
int i;
int q = digits[ i = (nDigits-1)];
if ( q == '9" ){
while ( q == '9" && i > 0 ){
digits[i] = '0";
q = digits[--i];
}
if ( q == '9" ){
// carryout! High-order 1, rest 0s, larger exp.
decExponent += 1;
digits[0] = '1";
return;
}
// else fall through.
}
digits[i] = (char)(q+1);
| | float | stickyRound(double dval)
long lbits = Double.doubleToLongBits( dval );
long binexp = lbits & expMask;
if ( binexp == 0L || binexp == expMask ){
// what we have here is special.
// don't worry, the right thing will happen.
return (float) dval;
}
lbits += (long)roundDir; // hack-o-matic.
return (float)Double.longBitsToDouble( lbits );
| | public java.lang.String | toJavaFormatString()
char result[] = new char[ nDigits + 10 ];
int i = 0;
if ( isNegative ){ result[0] = '-"; i = 1; }
if ( isExceptional ){
System.arraycopy( digits, 0, result, i, nDigits );
i += nDigits;
} else {
if ( decExponent > 0 && decExponent < 8 ){
// print digits.digits.
int charLength = Math.min( nDigits, decExponent );
System.arraycopy( digits, 0, result, i, charLength );
i += charLength;
if ( charLength < decExponent ){
charLength = decExponent-charLength;
System.arraycopy( zero, 0, result, i, charLength );
i += charLength;
result[i++] = '.";
result[i++] = '0";
} else {
result[i++] = '.";
if ( charLength < nDigits ){
int t = nDigits - charLength;
System.arraycopy( digits, charLength, result, i, t );
i += t;
} else{
result[i++] = '0";
}
}
} else if ( decExponent <=0 && decExponent > -3 ){
result[i++] = '0";
result[i++] = '.";
if ( decExponent != 0 ){
System.arraycopy( zero, 0, result, i, -decExponent );
i -= decExponent;
}
System.arraycopy( digits, 0, result, i, nDigits );
i += nDigits;
} else {
result[i++] = digits[0];
result[i++] = '.";
if ( nDigits > 1 ){
System.arraycopy( digits, 1, result, i, nDigits-1 );
i += nDigits-1;
} else {
result[i++] = '0";
}
result[i++] = 'E";
int e;
if ( decExponent <= 0 ){
result[i++] = '-";
e = -decExponent+1;
} else {
e = decExponent-1;
}
// decExponent has 1, 2, or 3, digits
if ( e <= 9 ) {
result[i++] = (char)( e+'0" );
} else if ( e <= 99 ){
result[i++] = (char)( e/10 +'0" );
result[i++] = (char)( e%10 + '0" );
} else {
result[i++] = (char)(e/100+'0");
e %= 100;
result[i++] = (char)(e/10+'0");
result[i++] = (char)( e%10 + '0" );
}
}
}
return new String(result, 0, i);
| | public java.lang.String | toString()
// most brain-dead version
StringBuffer result = new StringBuffer( nDigits+8 );
if ( isNegative ){ result.append( '-" ); }
if ( isExceptional ){
result.append( digits, 0, nDigits );
} else {
result.append( "0.");
result.append( digits, 0, nDigits );
result.append('e");
result.append( decExponent );
}
return new String(result);
| | private static double | ulp(double dval, boolean subtracting)
long lbits = Double.doubleToLongBits( dval ) & ~signMask;
int binexp = (int)(lbits >>> expShift);
double ulpval;
if ( subtracting && ( binexp >= expShift ) && ((lbits&fractMask) == 0L) ){
// for subtraction from normalized, powers of 2,
// use next-smaller exponent
binexp -= 1;
}
if ( binexp > expShift ){
ulpval = Double.longBitsToDouble( ((long)(binexp-expShift))<<expShift );
} else if ( binexp == 0 ){
ulpval = Double.MIN_VALUE;
} else {
ulpval = Double.longBitsToDouble( 1L<<(binexp-1) );
}
if ( subtracting ) ulpval = - ulpval;
return ulpval;
|
|
