in absl/strings/numbers.cc [366:482]
static ExpDigits SplitToSix(const double value) {
ExpDigits exp_dig;
int exp = 5;
double d = value;
// First step: calculate a close approximation of the output, where the
// value d will be between 100,000 and 999,999, representing the digits
// in the output ASCII array, and exp is the base-10 exponent. It would be
// faster to use a table here, and to look up the base-2 exponent of value,
// however value is an IEEE-754 64-bit number, so the table would have 2,000
// entries, which is not cache-friendly.
if (d >= 999999.5) {
if (d >= 1e+261) exp += 256, d *= 1e-256;
if (d >= 1e+133) exp += 128, d *= 1e-128;
if (d >= 1e+69) exp += 64, d *= 1e-64;
if (d >= 1e+37) exp += 32, d *= 1e-32;
if (d >= 1e+21) exp += 16, d *= 1e-16;
if (d >= 1e+13) exp += 8, d *= 1e-8;
if (d >= 1e+9) exp += 4, d *= 1e-4;
if (d >= 1e+7) exp += 2, d *= 1e-2;
if (d >= 1e+6) exp += 1, d *= 1e-1;
} else {
if (d < 1e-250) exp -= 256, d *= 1e256;
if (d < 1e-122) exp -= 128, d *= 1e128;
if (d < 1e-58) exp -= 64, d *= 1e64;
if (d < 1e-26) exp -= 32, d *= 1e32;
if (d < 1e-10) exp -= 16, d *= 1e16;
if (d < 1e-2) exp -= 8, d *= 1e8;
if (d < 1e+2) exp -= 4, d *= 1e4;
if (d < 1e+4) exp -= 2, d *= 1e2;
if (d < 1e+5) exp -= 1, d *= 1e1;
}
// At this point, d is in the range [99999.5..999999.5) and exp is in the
// range [-324..308]. Since we need to round d up, we want to add a half
// and truncate.
// However, the technique above may have lost some precision, due to its
// repeated multiplication by constants that each may be off by half a bit
// of precision. This only matters if we're close to the edge though.
// Since we'd like to know if the fractional part of d is close to a half,
// we multiply it by 65536 and see if the fractional part is close to 32768.
// (The number doesn't have to be a power of two,but powers of two are faster)
uint64_t d64k = d * 65536;
int dddddd; // A 6-digit decimal integer.
if ((d64k % 65536) == 32767 || (d64k % 65536) == 32768) {
// OK, it's fairly likely that precision was lost above, which is
// not a surprise given only 52 mantissa bits are available. Therefore
// redo the calculation using 128-bit numbers. (64 bits are not enough).
// Start out with digits rounded down; maybe add one below.
dddddd = static_cast<int>(d64k / 65536);
// mantissa is a 64-bit integer representing M.mmm... * 2^63. The actual
// value we're representing, of course, is M.mmm... * 2^exp2.
int exp2;
double m = std::frexp(value, &exp2);
uint64_t mantissa = m * (32768.0 * 65536.0 * 65536.0 * 65536.0);
// std::frexp returns an m value in the range [0.5, 1.0), however we
// can't multiply it by 2^64 and convert to an integer because some FPUs
// throw an exception when converting an number higher than 2^63 into an
// integer - even an unsigned 64-bit integer! Fortunately it doesn't matter
// since m only has 52 significant bits anyway.
mantissa <<= 1;
exp2 -= 64; // not needed, but nice for debugging
// OK, we are here to compare:
// (dddddd + 0.5) * 10^(exp-5) vs. mantissa * 2^exp2
// so we can round up dddddd if appropriate. Those values span the full
// range of 600 orders of magnitude of IEE 64-bit floating-point.
// Fortunately, we already know they are very close, so we don't need to
// track the base-2 exponent of both sides. This greatly simplifies the
// the math since the 2^exp2 calculation is unnecessary and the power-of-10
// calculation can become a power-of-5 instead.
std::pair<uint64_t, uint64_t> edge, val;
if (exp >= 6) {
// Compare (dddddd + 0.5) * 5 ^ (exp - 5) to mantissa
// Since we're tossing powers of two, 2 * dddddd + 1 is the
// same as dddddd + 0.5
edge = PowFive(2 * dddddd + 1, exp - 5);
val.first = mantissa;
val.second = 0;
} else {
// We can't compare (dddddd + 0.5) * 5 ^ (exp - 5) to mantissa as we did
// above because (exp - 5) is negative. So we compare (dddddd + 0.5) to
// mantissa * 5 ^ (5 - exp)
edge = PowFive(2 * dddddd + 1, 0);
val = PowFive(mantissa, 5 - exp);
}
// printf("exp=%d %016lx %016lx vs %016lx %016lx\n", exp, val.first,
// val.second, edge.first, edge.second);
if (val > edge) {
dddddd++;
} else if (val == edge) {
dddddd += (dddddd & 1);
}
} else {
// Here, we are not close to the edge.
dddddd = static_cast<int>((d64k + 32768) / 65536);
}
if (dddddd == 1000000) {
dddddd = 100000;
exp += 1;
}
exp_dig.exponent = exp;
int two_digits = dddddd / 10000;
dddddd -= two_digits * 10000;
numbers_internal::PutTwoDigits(two_digits, &exp_dig.digits[0]);
two_digits = dddddd / 100;
dddddd -= two_digits * 100;
numbers_internal::PutTwoDigits(two_digits, &exp_dig.digits[2]);
numbers_internal::PutTwoDigits(dddddd, &exp_dig.digits[4]);
return exp_dig;
}