case COMPATIBLE: break; case DECIMAL_FLOAT: exp = checkExponent(decExponent + precision); digits = applyPrecision(decExponent + precision); break; case SCIENTIFIC: exp = checkExponent(precision + 1); digits = applyPrecision(precision + 1); break; case GENERAL: exp = checkExponent(precision); digits = applyPrecision(precision); // adjust precision to be the number of digits to right of decimal // the real exponent to be output is actually exp - 1, not exp if (exp - 1 < -4 || exp - 1 >= precision) { form = Form.SCIENTIFIC; precision--; } else { form = Form.DECIMAL_FLOAT; precision = precision - exp; } break; default: assert false; } decExponentRounded = exp; if (exp > 0 && ((form == Form.COMPATIBLE && (exp < 8)) || (form == Form.DECIMAL_FLOAT))) { // print digits.digits. int charLength = Math.min(nDigits, exp); System.arraycopy(digits, 0, result, i, charLength); i += charLength; if (charLength < exp) { charLength = exp-charLength; for (int nz = 0; nz < charLength; nz++) result[i++] = '0'; // Do not append ".0" for formatted floats since the user // may request that it be omitted. It is added as necessary // by the Formatter. if (form == Form.COMPATIBLE) { result[i++] = '.'; result[i++] = '0'; } } else { // Do not append ".0" for formatted floats since the user // may request that it be omitted. It is added as necessary // by the Formatter. if (form == Form.COMPATIBLE) { result[i++] = '.'; if (charLength < nDigits) { int t = Math.min(nDigits - charLength, precision); System.arraycopy(digits, charLength, result, i, t); i += t; } else { result[i++] = '0'; } } else { int t = Math.min(nDigits - charLength, precision); if (t > 0) { result[i++] = '.'; System.arraycopy(digits, charLength, result, i, t); i += t; } } } } else if (exp <= 0 && ((form == Form.COMPATIBLE && exp > -3) || (form == Form.DECIMAL_FLOAT))) { // print 0.0* digits result[i++] = '0'; if (exp != 0) { // write '0' s before the significant digits int t = Math.min(-exp, precision); if (t > 0) { result[i++] = '.'; for (int nz = 0; nz < t; nz++) result[i++] = '0'; } } int t = Math.min(digits.length, precision + exp); if (t > 0) { if (i == 1) result[i++] = '.'; // copy only when significant digits are within the precision System.arraycopy(digits, 0, result, i, t); i += t; } } else { result[i++] = digits[0]; if (form == Form.COMPATIBLE) { result[i++] = '.'; if (nDigits > 1) { System.arraycopy(digits, 1, result, i, nDigits-1); i += nDigits-1; } else { result[i++] = '0'; } result[i++] = 'E'; } else { if (nDigits > 1) { int t = Math.min(nDigits -1, precision); if (t > 0) { result[i++] = '.'; System.arraycopy(digits, 1, result, i, t); i += t; } } result[i++] = 'e'; } int e; if (exp <= 0) { result[i++] = '-'; e = -exp+1; } else { if (form != Form.COMPATIBLE) result[i++] = '+'; e = exp-1; } // decExponent has 1, 2, or 3, digits if (e <= 9) { if (form != Form.COMPATIBLE) result[i++] = '0'; 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 i; } // Per-thread buffer for string/stringbuffer conversion private static ThreadLocal perThreadBuffer = new ThreadLocal() { protected synchronized Object initialValue() { return new char[26]; } }; /* * Take a FormattedFloatingDecimal, which we presumably just scanned in, * and find out what its value is, as a double. * * AS A SIDE EFFECT, SET roundDir TO INDICATE PREFERRED * ROUNDING DIRECTION in case the result is really destined * for a single-precision float. */ public strictfp double doubleValue(){ int kDigits = Math.min( nDigits, maxDecimalDigits+1 ); long lValue; double dValue; double rValue, tValue; // First, check for NaN and Infinity values if(digits == infinity || digits == notANumber) { if(digits == notANumber) return Double.NaN; else return (isNegative?Double.NEGATIVE_INFINITY:Double.POSITIVE_INFINITY); } else { if (mustSetRoundDir) { 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 multipy 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 fequency break correctionLoop; } FDBigInt halfUlp = constructPow52( B5, Ulp2 ); if ( (cmpResult = diff.cmp( halfUlp ) ) < 0 ){ // difference is small. // this is close enough if (mustSetRoundDir) { 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?? if (mustSetRoundDir) { 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.Lines 1001–1500 of 1767. Previous lines Load more
case COMPATIBLE: break; case DECIMAL_FLOAT: exp = checkExponent(decExponent + precision); digits = applyPrecision(decExponent + precision); break; case SCIENTIFIC: exp = checkExponent(precision + 1); digits = applyPrecision(precision + 1); break; case GENERAL: exp = checkExponent(precision); digits = applyPrecision(precision); // adjust precision to be the number of digits to right of decimal // the real exponent to be output is actually exp - 1, not exp if (exp - 1 < -4 || exp - 1 >= precision) { form = Form.SCIENTIFIC; precision--; } else { form = Form.DECIMAL_FLOAT; precision = precision - exp; } break; default: assert false; } decExponentRounded = exp; if (exp > 0 && ((form == Form.COMPATIBLE && (exp < 8)) || (form == Form.DECIMAL_FLOAT))) { // print digits.digits. int charLength = Math.min(nDigits, exp); System.arraycopy(digits, 0, result, i, charLength); i += charLength; if (charLength < exp) { charLength = exp-charLength; for (int nz = 0; nz < charLength; nz++) result[i++] = '0'; // Do not append ".0" for formatted floats since the user // may request that it be omitted. It is added as necessary // by the Formatter. if (form == Form.COMPATIBLE) { result[i++] = '.'; result[i++] = '0'; } } else { // Do not append ".0" for formatted floats since the user // may request that it be omitted. It is added as necessary // by the Formatter. if (form == Form.COMPATIBLE) { result[i++] = '.'; if (charLength < nDigits) { int t = Math.min(nDigits - charLength, precision); System.arraycopy(digits, charLength, result, i, t); i += t; } else { result[i++] = '0'; } } else { int t = Math.min(nDigits - charLength, precision); if (t > 0) { result[i++] = '.'; System.arraycopy(digits, charLength, result, i, t); i += t; } } } } else if (exp <= 0 && ((form == Form.COMPATIBLE && exp > -3) || (form == Form.DECIMAL_FLOAT))) { // print 0.0* digits result[i++] = '0'; if (exp != 0) { // write '0' s before the significant digits int t = Math.min(-exp, precision); if (t > 0) { result[i++] = '.'; for (int nz = 0; nz < t; nz++) result[i++] = '0'; } } int t = Math.min(digits.length, precision + exp); if (t > 0) { if (i == 1) result[i++] = '.'; // copy only when significant digits are within the precision System.arraycopy(digits, 0, result, i, t); i += t; } } else { result[i++] = digits[0]; if (form == Form.COMPATIBLE) { result[i++] = '.'; if (nDigits > 1) { System.arraycopy(digits, 1, result, i, nDigits-1); i += nDigits-1; } else { result[i++] = '0'; } result[i++] = 'E'; } else { if (nDigits > 1) { int t = Math.min(nDigits -1, precision); if (t > 0) { result[i++] = '.'; System.arraycopy(digits, 1, result, i, t); i += t; } } result[i++] = 'e'; } int e; if (exp <= 0) { result[i++] = '-'; e = -exp+1; } else { if (form != Form.COMPATIBLE) result[i++] = '+'; e = exp-1; } // decExponent has 1, 2, or 3, digits if (e <= 9) { if (form != Form.COMPATIBLE) result[i++] = '0'; 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 i; } // Per-thread buffer for string/stringbuffer conversion private static ThreadLocal perThreadBuffer = new ThreadLocal() { protected synchronized Object initialValue() { return new char[26]; } }; /* * Take a FormattedFloatingDecimal, which we presumably just scanned in, * and find out what its value is, as a double. * * AS A SIDE EFFECT, SET roundDir TO INDICATE PREFERRED * ROUNDING DIRECTION in case the result is really destined * for a single-precision float. */ public strictfp double doubleValue(){ int kDigits = Math.min( nDigits, maxDecimalDigits+1 ); long lValue; double dValue; double rValue, tValue; // First, check for NaN and Infinity values if(digits == infinity || digits == notANumber) { if(digits == notANumber) return Double.NaN; else return (isNegative?Double.NEGATIVE_INFINITY:Double.POSITIVE_INFINITY); } else { if (mustSetRoundDir) { 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 multipy 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 fequency break correctionLoop; } FDBigInt halfUlp = constructPow52( B5, Ulp2 ); if ( (cmpResult = diff.cmp( halfUlp ) ) < 0 ){ // difference is small. // this is close enough if (mustSetRoundDir) { 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?? if (mustSetRoundDir) { 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.