// Copyright 2018 The Abseil Authors. // // Licensed under the Apache License, Version 2.0 (the "License"); // you may not use this file except in compliance with the License. // You may obtain a copy of the License at // // https://www.apache.org/licenses/LICENSE-2.0 // // Unless required by applicable law or agreed to in writing, software // distributed under the License is distributed on an "AS IS" BASIS, // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. // See the License for the specific language governing permissions and // limitations under the License.
// The macro ABSL_BIT_PACK_FLOATS is defined on x86-64, where IEEE floating // point numbers have the same endianness in memory as a bitfield struct // containing the corresponding parts. // // When set, we replace calls to ldexp() with manual bit packing, which is // faster and is unaffected by floating point environment. #ifdef ABSL_BIT_PACK_FLOATS #error ABSL_BIT_PACK_FLOATS cannot be directly set #elifdefined(__x86_64__) || defined(_M_X64) #define ABSL_BIT_PACK_FLOATS 1 #endif
// A note about subnormals: // // The code below talks about "normals" and "subnormals". A normal IEEE float // has a fixed-width mantissa and power of two exponent. For example, a normal // `double` has a 53-bit mantissa. Because the high bit is always 1, it is not // stored in the representation. The implicit bit buys an extra bit of // resolution in the datatype. // // The downside of this scheme is that there is a large gap between DBL_MIN and // zero. (Large, at least, relative to the different between DBL_MIN and the // next representable number). This gap is softened by the "subnormal" numbers, // which have the same power-of-two exponent as DBL_MIN, but no implicit 53rd // bit. An all-bits-zero exponent in the encoding represents subnormals. (Zero // is represented as a subnormal with an all-bits-zero mantissa.) // // The code below, in calculations, represents the mantissa as a uint64_t. The // end result normally has the 53rd bit set. It represents subnormals by using // narrower mantissas.
namespace absl {
ABSL_NAMESPACE_BEGIN namespace {
template <typename FloatType> struct FloatTraits;
template <> struct FloatTraits<double> { using mantissa_t = uint64_t;
// The number of bits in the given float type. static constexpr int kTargetBits = 64;
// The number of exponent bits in the given float type. static constexpr int kTargetExponentBits = 11;
// The number of mantissa bits in the given float type. This includes the // implied high bit. static constexpr int kTargetMantissaBits = 53;
// The largest supported IEEE exponent, in our integral mantissa // representation. // // If `m` is the largest possible int kTargetMantissaBits bits wide, then // m * 2**kMaxExponent is exactly equal to DBL_MAX. static constexpr int kMaxExponent = 971;
// The smallest supported IEEE normal exponent, in our integral mantissa // representation. // // If `m` is the smallest possible int kTargetMantissaBits bits wide, then // m * 2**kMinNormalExponent is exactly equal to DBL_MIN. static constexpr int kMinNormalExponent = -1074;
// The IEEE exponent bias. It equals ((1 << (kTargetExponentBits - 1)) - 1). static constexpr int kExponentBias = 1023;
// The Eisel-Lemire "Shifting to 54/25 Bits" adjustment. It equals (63 - 1 - // kTargetMantissaBits). static constexpr int kEiselLemireShift = 9;
// The Eisel-Lemire high64_mask. It equals ((1 << kEiselLemireShift) - 1). static constexpr uint64_t kEiselLemireMask = uint64_t{0x1FF};
// The smallest negative integer N (smallest negative means furthest from // zero) such that parsing 9999999999999999999eN, with 19 nines, is still // positive. Parsing a smaller (more negative) N will produce zero. // // Adjusting the decimal point and exponent, without adjusting the value, // 9999999999999999999eN equals 9.999999999999999999eM where M = N + 18. // // 9999999999999999999, with 19 nines but no decimal point, is the largest // "repeated nines" integer that fits in a uint64_t. static constexpr int kEiselLemireMinInclusiveExp10 = -324 - 18;
// The smallest positive integer N such that parsing 1eN produces infinity. // Parsing a smaller N will produce something finite. static constexpr int kEiselLemireMaxExclusiveExp10 = 309;
staticdouble MakeNan(absl::Nonnull<constchar*> tagp) { #if ABSL_HAVE_BUILTIN(__builtin_nan) // Use __builtin_nan() if available since it has a fix for // https://bugs.llvm.org/show_bug.cgi?id=37778 // std::nan may use the glibc implementation. return __builtin_nan(tagp); #else // Support nan no matter which namespace it's in. Some platforms // incorrectly don't put it in namespace std. usingnamespace std; // NOLINT return nan(tagp); #endif
}
// Builds a nonzero floating point number out of the provided parts. // // This is intended to do the same operation as ldexp(mantissa, exponent), // but using purely integer math, to avoid -ffastmath and floating // point environment issues. Using type punning is also faster. We fall back // to ldexp on a per-platform basis for portability. // // `exponent` must be between kMinNormalExponent and kMaxExponent. // // `mantissa` must either be exactly kTargetMantissaBits wide, in which case // a normal value is made, or it must be less narrow than that, in which case // `exponent` must be exactly kMinNormalExponent, and a subnormal value is // made. staticdouble Make(mantissa_t mantissa, int exponent, bool sign) { #ifndef ABSL_BIT_PACK_FLOATS // Support ldexp no matter which namespace it's in. Some platforms // incorrectly don't put it in namespace std. usingnamespace std; // NOLINT return sign ? -ldexp(mantissa, exponent) : ldexp(mantissa, exponent); #else
constexpr uint64_t kMantissaMask =
(uint64_t{1} << (kTargetMantissaBits - 1)) - 1;
uint64_t dbl = static_cast<uint64_t>(sign) << 63; if (mantissa > kMantissaMask) { // Normal value. // Adjust by 1023 for the exponent representation bias, and an additional // 52 due to the implied decimal point in the IEEE mantissa // representation.
dbl += static_cast<uint64_t>(exponent + 1023 + kTargetMantissaBits - 1)
<< 52;
mantissa &= kMantissaMask;
} else { // subnormal value
assert(exponent == kMinNormalExponent);
}
dbl += mantissa; return absl::bit_cast<double>(dbl); #endif// ABSL_BIT_PACK_FLOATS
}
};
// Specialization of floating point traits for the `float` type. See the // FloatTraits<double> specialization above for meaning of each of the following // members and methods. template <> struct FloatTraits<float> { using mantissa_t = uint32_t;
static constexpr int kTargetBits = 32; static constexpr int kTargetExponentBits = 8; static constexpr int kTargetMantissaBits = 24; static constexpr int kMaxExponent = 104; static constexpr int kMinNormalExponent = -149; static constexpr int kExponentBias = 127; static constexpr int kEiselLemireShift = 38; static constexpr uint64_t kEiselLemireMask = uint64_t{0x3FFFFFFFFF}; static constexpr int kEiselLemireMinInclusiveExp10 = -46 - 18; static constexpr int kEiselLemireMaxExclusiveExp10 = 39;
staticfloat MakeNan(absl::Nonnull<constchar*> tagp) { #if ABSL_HAVE_BUILTIN(__builtin_nanf) // Use __builtin_nanf() if available since it has a fix for // https://bugs.llvm.org/show_bug.cgi?id=37778 // std::nanf may use the glibc implementation. return __builtin_nanf(tagp); #else // Support nanf no matter which namespace it's in. Some platforms // incorrectly don't put it in namespace std. usingnamespace std; // NOLINT return std::nanf(tagp); #endif
}
staticfloat Make(mantissa_t mantissa, int exponent, bool sign) { #ifndef ABSL_BIT_PACK_FLOATS // Support ldexpf no matter which namespace it's in. Some platforms // incorrectly don't put it in namespace std. usingnamespace std; // NOLINT return sign ? -ldexpf(mantissa, exponent) : ldexpf(mantissa, exponent); #else
constexpr uint32_t kMantissaMask =
(uint32_t{1} << (kTargetMantissaBits - 1)) - 1;
uint32_t flt = static_cast<uint32_t>(sign) << 31; if (mantissa > kMantissaMask) { // Normal value. // Adjust by 127 for the exponent representation bias, and an additional // 23 due to the implied decimal point in the IEEE mantissa // representation.
flt += static_cast<uint32_t>(exponent + 127 + kTargetMantissaBits - 1)
<< 23;
mantissa &= kMantissaMask;
} else { // subnormal value
assert(exponent == kMinNormalExponent);
}
flt += mantissa; return absl::bit_cast<float>(flt); #endif// ABSL_BIT_PACK_FLOATS
}
};
// Decimal-to-binary conversions require coercing powers of 10 into a mantissa // and a power of 2. The two helper functions Power10Mantissa(n) and // Power10Exponent(n) perform this task. Together, these represent a hand- // rolled floating point value which is equal to or just less than 10**n. // // The return values satisfy two range guarantees: // // Power10Mantissa(n) * 2**Power10Exponent(n) <= 10**n // < (Power10Mantissa(n) + 1) * 2**Power10Exponent(n) // // 2**63 <= Power10Mantissa(n) < 2**64. // // See the "Table of powers of 10" comment below for a "1e60" example. // // Lookups into the power-of-10 table must first check the Power10Overflow() and // Power10Underflow() functions, to avoid out-of-bounds table access. // // Indexes into these tables are biased by -kPower10TableMinInclusive. Valid // indexes range from kPower10TableMinInclusive to kPower10TableMaxExclusive. externconst uint64_t kPower10MantissaHighTable[]; // High 64 of 128 bits. externconst uint64_t kPower10MantissaLowTable[]; // Low 64 of 128 bits.
// The smallest (inclusive) allowed value for use with the Power10Mantissa() // and Power10Exponent() functions below. (If a smaller exponent is needed in // calculations, the end result is guaranteed to underflow.)
constexpr int kPower10TableMinInclusive = -342;
// The largest (exclusive) allowed value for use with the Power10Mantissa() and // Power10Exponent() functions below. (If a larger-or-equal exponent is needed // in calculations, the end result is guaranteed to overflow.)
constexpr int kPower10TableMaxExclusive = 309;
// Returns true if n is large enough that 10**n always results in an IEEE // overflow. bool Power10Overflow(int n) { return n >= kPower10TableMaxExclusive; }
// Returns true if n is small enough that 10**n times a ParsedFloat mantissa // always results in an IEEE underflow. bool Power10Underflow(int n) { return n < kPower10TableMinInclusive; }
// Returns true if Power10Mantissa(n) * 2**Power10Exponent(n) is exactly equal // to 10**n numerically. Put another way, this returns true if there is no // truncation error in Power10Mantissa(n). bool Power10Exact(int n) { return n >= 0 && n <= 27; }
// Sentinel exponent values for representing numbers too large or too close to // zero to represent in a double.
constexpr int kOverflow = 99999;
constexpr int kUnderflow = -99999;
// Struct representing the calculated conversion result of a positive (nonzero) // floating point number. // // The calculated number is mantissa * 2**exponent (mantissa is treated as an // integer.) `mantissa` is chosen to be the correct width for the IEEE float // representation being calculated. (`mantissa` will always have the same bit // width for normal values, and narrower bit widths for subnormals.) // // If the result of conversion was an underflow or overflow, exponent is set // to kUnderflow or kOverflow. struct CalculatedFloat {
uint64_t mantissa = 0; int exponent = 0;
};
// Returns the bit width of the given uint128. (Equivalently, returns 128 // minus the number of leading zero bits.) int BitWidth(uint128 value) { if (Uint128High64(value) == 0) { // This static_cast is only needed when using a std::bit_width() // implementation that does not have the fix for LWG 3656 applied. returnstatic_cast<int>(bit_width(Uint128Low64(value)));
} return128 - countl_zero(Uint128High64(value));
}
// Calculates how far to the right a mantissa needs to be shifted to create a // properly adjusted mantissa for an IEEE floating point number. // // `mantissa_width` is the bit width of the mantissa to be shifted, and // `binary_exponent` is the exponent of the number before the shift. // // This accounts for subnormal values, and will return a larger-than-normal // shift if binary_exponent would otherwise be too low. template <typename FloatType> int NormalizedShiftSize(int mantissa_width, int binary_exponent) { constint normal_shift =
mantissa_width - FloatTraits<FloatType>::kTargetMantissaBits; constint minimum_shift =
FloatTraits<FloatType>::kMinNormalExponent - binary_exponent; return std::max(normal_shift, minimum_shift);
}
// Right shifts a uint128 so that it has the requested bit width. (The // resulting value will have 128 - bit_width leading zeroes.) The initial // `value` must be wider than the requested bit width. // // Returns the number of bits shifted. int TruncateToBitWidth(int bit_width, absl::Nonnull<uint128*> value) { constint current_bit_width = BitWidth(*value); constint shift = current_bit_width - bit_width;
*value >>= shift; return shift;
}
// Checks if the given ParsedFloat represents one of the edge cases that are // not dependent on number base: zero, infinity, or NaN. If so, sets *value // the appropriate double, and returns true. template <typename FloatType> bool HandleEdgeCase(const strings_internal::ParsedFloat& input, bool negative,
absl::Nonnull<FloatType*> value) { if (input.type == strings_internal::FloatType::kNan) { // A bug in both clang < 7 and gcc would cause the compiler to optimize // away the buffer we are building below. Declaring the buffer volatile // avoids the issue, and has no measurable performance impact in // microbenchmarks. // // https://bugs.llvm.org/show_bug.cgi?id=37778 // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=86113
constexpr ptrdiff_t kNanBufferSize = 128; #if (defined(__GNUC__) && !defined(__clang__)) || \
(defined(__clang__) && __clang_major__ < 7) volatilechar n_char_sequence[kNanBufferSize]; #else char n_char_sequence[kNanBufferSize]; #endif if (input.subrange_begin == nullptr) {
n_char_sequence[0] = '\0';
} else {
ptrdiff_t nan_size = input.subrange_end - input.subrange_begin;
nan_size = std::min(nan_size, kNanBufferSize - 1);
std::copy_n(input.subrange_begin, nan_size, n_char_sequence);
n_char_sequence[nan_size] = '\0';
} char* nan_argument = const_cast<char*>(n_char_sequence);
*value = negative ? -FloatTraits<FloatType>::MakeNan(nan_argument)
: FloatTraits<FloatType>::MakeNan(nan_argument); returntrue;
} if (input.type == strings_internal::FloatType::kInfinity) {
*value = negative ? -std::numeric_limits<FloatType>::infinity()
: std::numeric_limits<FloatType>::infinity(); returntrue;
} if (input.mantissa == 0) {
*value = negative ? -0.0 : 0.0; returntrue;
} returnfalse;
}
// Given a CalculatedFloat result of a from_chars conversion, generate the // correct output values. // // CalculatedFloat can represent an underflow or overflow, in which case the // error code in *result is set. Otherwise, the calculated floating point // number is stored in *value. template <typename FloatType> void EncodeResult(const CalculatedFloat& calculated, bool negative,
absl::Nonnull<absl::from_chars_result*> result,
absl::Nonnull<FloatType*> value) { if (calculated.exponent == kOverflow) {
result->ec = std::errc::result_out_of_range;
*value = negative ? -std::numeric_limits<FloatType>::max()
: std::numeric_limits<FloatType>::max(); return;
} elseif (calculated.mantissa == 0 || calculated.exponent == kUnderflow) {
result->ec = std::errc::result_out_of_range;
*value = negative ? -0.0 : 0.0; return;
}
*value = FloatTraits<FloatType>::Make( static_cast<typename FloatTraits<FloatType>::mantissa_t>(
calculated.mantissa),
calculated.exponent, negative);
}
// Returns the given uint128 shifted to the right by `shift` bits, and rounds // the remaining bits using round_to_nearest logic. The value is returned as a // uint64_t, since this is the type used by this library for storing calculated // floating point mantissas. // // It is expected that the width of the input value shifted by `shift` will // be the correct bit-width for the target mantissa, which is strictly narrower // than a uint64_t. // // If `input_exact` is false, then a nonzero error epsilon is assumed. For // rounding purposes, the true value being rounded is strictly greater than the // input value. The error may represent a single lost carry bit. // // When input_exact, shifted bits of the form 1000000... represent a tie, which // is broken by rounding to even -- the rounding direction is chosen so the low // bit of the returned value is 0. // // When !input_exact, shifted bits of the form 10000000... represent a value // strictly greater than one half (due to the error epsilon), and so ties are // always broken by rounding up. // // When !input_exact, shifted bits of the form 01111111... are uncertain; // the true value may or may not be greater than 10000000..., due to the // possible lost carry bit. The correct rounding direction is unknown. In this // case, the result is rounded down, and `output_exact` is set to false. // // Zero and negative values of `shift` are accepted, in which case the word is // shifted left, as necessary.
uint64_t ShiftRightAndRound(uint128 value, int shift, bool input_exact,
absl::Nonnull<bool*> output_exact) { if (shift <= 0) {
*output_exact = input_exact; returnstatic_cast<uint64_t>(value << -shift);
} if (shift >= 128) { // Exponent is so small that we are shifting away all significant bits. // Answer will not be representable, even as a subnormal, so return a zero // mantissa (which represents underflow).
*output_exact = true; return0;
}
const uint128 shifted_bits = value & shift_mask;
value >>= shift; if (shifted_bits > halfway_point) { // Shifted bits greater than 10000... require rounding up. returnstatic_cast<uint64_t>(value + 1);
} if (shifted_bits == halfway_point) { // In exact mode, shifted bits of 10000... mean we're exactly halfway // between two numbers, and we must round to even. So only round up if // the low bit of `value` is set. // // In inexact mode, the nonzero error means the actual value is greater // than the halfway point and we must always round up. if ((value & 1) == 1 || !input_exact) {
++value;
} returnstatic_cast<uint64_t>(value);
} if (!input_exact && shifted_bits == halfway_point - 1) { // Rounding direction is unclear, due to error.
*output_exact = false;
} // Otherwise, round down. returnstatic_cast<uint64_t>(value);
}
// Checks if a floating point guess needs to be rounded up, using high precision // math. // // `guess_mantissa` and `guess_exponent` represent a candidate guess for the // number represented by `parsed_decimal`. // // The exact number represented by `parsed_decimal` must lie between the two // numbers: // A = `guess_mantissa * 2**guess_exponent` // B = `(guess_mantissa + 1) * 2**guess_exponent` // // This function returns false if `A` is the better guess, and true if `B` is // the better guess, with rounding ties broken by rounding to even. bool MustRoundUp(uint64_t guess_mantissa, int guess_exponent, const strings_internal::ParsedFloat& parsed_decimal) { // 768 is the number of digits needed in the worst case. We could determine a // better limit dynamically based on the value of parsed_decimal.exponent. // This would optimize pathological input cases only. (Sane inputs won't have // hundreds of digits of mantissa.)
absl::strings_internal::BigUnsigned<84> exact_mantissa; int exact_exponent = exact_mantissa.ReadFloatMantissa(parsed_decimal, 768);
// Adjust the `guess` arguments to be halfway between A and B.
guess_mantissa = guess_mantissa * 2 + 1;
guess_exponent -= 1;
// In our comparison: // lhs = exact = exact_mantissa * 10**exact_exponent // = exact_mantissa * 5**exact_exponent * 2**exact_exponent // rhs = guess = guess_mantissa * 2**guess_exponent // // Because we are doing integer math, we can't directly deal with negative // exponents. We instead move these to the other side of the inequality.
absl::strings_internal::BigUnsigned<84>& lhs = exact_mantissa; int comparison; if (exact_exponent >= 0) {
lhs.MultiplyByFiveToTheNth(exact_exponent);
absl::strings_internal::BigUnsigned<84> rhs(guess_mantissa); // There are powers of 2 on both sides of the inequality; reduce this to // a single bit-shift. if (exact_exponent > guess_exponent) {
lhs.ShiftLeft(exact_exponent - guess_exponent);
} else {
rhs.ShiftLeft(guess_exponent - exact_exponent);
}
comparison = Compare(lhs, rhs);
} else { // Move the power of 5 to the other side of the equation, giving us: // lhs = exact_mantissa * 2**exact_exponent // rhs = guess_mantissa * 5**(-exact_exponent) * 2**guess_exponent
absl::strings_internal::BigUnsigned<84> rhs =
absl::strings_internal::BigUnsigned<84>::FiveToTheNth(-exact_exponent);
rhs.MultiplyBy(guess_mantissa); if (exact_exponent > guess_exponent) {
lhs.ShiftLeft(exact_exponent - guess_exponent);
} else {
rhs.ShiftLeft(guess_exponent - exact_exponent);
}
comparison = Compare(lhs, rhs);
} if (comparison < 0) { returnfalse;
} elseif (comparison > 0) { returntrue;
} else { // When lhs == rhs, the decimal input is exactly between A and B. // Round towards even -- round up only if the low bit of the initial // `guess_mantissa` was a 1. We shifted guess_mantissa left 1 bit at // the beginning of this function, so test the 2nd bit here. return (guess_mantissa & 2) == 2;
}
}
// Constructs a CalculatedFloat from a given mantissa and exponent, but // with the following normalizations applied: // // If rounding has caused mantissa to increase just past the allowed bit // width, shift and adjust exponent. // // If exponent is too high, sets kOverflow. // // If mantissa is zero (representing a non-zero value not representable, even // as a subnormal), sets kUnderflow. template <typename FloatType>
CalculatedFloat CalculatedFloatFromRawValues(uint64_t mantissa, int exponent) {
CalculatedFloat result; if (mantissa == uint64_t{1} << FloatTraits<FloatType>::kTargetMantissaBits) {
mantissa >>= 1;
exponent += 1;
} if (exponent > FloatTraits<FloatType>::kMaxExponent) {
result.exponent = kOverflow;
} elseif (mantissa == 0) {
result.exponent = kUnderflow;
} else {
result.exponent = exponent;
result.mantissa = mantissa;
} return result;
}
template <typename FloatType>
CalculatedFloat CalculateFromParsedHexadecimal( const strings_internal::ParsedFloat& parsed_hex) {
uint64_t mantissa = parsed_hex.mantissa; int exponent = parsed_hex.exponent; // This static_cast is only needed when using a std::bit_width() // implementation that does not have the fix for LWG 3656 applied. int mantissa_width = static_cast<int>(bit_width(mantissa)); constint shift = NormalizedShiftSize<FloatType>(mantissa_width, exponent); bool result_exact;
exponent += shift;
mantissa = ShiftRightAndRound(mantissa, shift, /* input exact= */ true, &result_exact); // ParseFloat handles rounding in the hexadecimal case, so we don't have to // check `result_exact` here. return CalculatedFloatFromRawValues<FloatType>(mantissa, exponent);
}
// Large or small enough decimal exponents will always result in overflow // or underflow. if (Power10Underflow(parsed_decimal.exponent)) {
result.exponent = kUnderflow; return result;
} elseif (Power10Overflow(parsed_decimal.exponent)) {
result.exponent = kOverflow; return result;
}
// Otherwise convert our power of 10 into a power of 2 times an integer // mantissa, and multiply this by our parsed decimal mantissa.
uint128 wide_binary_mantissa = parsed_decimal.mantissa;
wide_binary_mantissa *= Power10Mantissa(parsed_decimal.exponent); int binary_exponent = Power10Exponent(parsed_decimal.exponent);
// Discard bits that are inaccurate due to truncation error. The magic // `mantissa_width` constants below are justified in // https://abseil.io/about/design/charconv. They represent the number of bits // in `wide_binary_mantissa` that are guaranteed to be unaffected by error // propagation. bool mantissa_exact; int mantissa_width; if (parsed_decimal.subrange_begin) { // Truncated mantissa
mantissa_width = 58;
mantissa_exact = false;
binary_exponent +=
TruncateToBitWidth(mantissa_width, &wide_binary_mantissa);
} elseif (!Power10Exact(parsed_decimal.exponent)) { // Exact mantissa, truncated power of ten
mantissa_width = 63;
mantissa_exact = false;
binary_exponent +=
TruncateToBitWidth(mantissa_width, &wide_binary_mantissa);
} else { // Product is exact
mantissa_width = BitWidth(wide_binary_mantissa);
mantissa_exact = true;
}
// Shift into an FloatType-sized mantissa, and round to nearest. constint shift =
NormalizedShiftSize<FloatType>(mantissa_width, binary_exponent); bool result_exact;
binary_exponent += shift;
uint64_t binary_mantissa = ShiftRightAndRound(wide_binary_mantissa, shift,
mantissa_exact, &result_exact); if (!result_exact) { // We could not determine the rounding direction using int128 math. Use // full resolution math instead. if (MustRoundUp(binary_mantissa, binary_exponent, parsed_decimal)) {
binary_mantissa += 1;
}
}
// As discussed in https://nigeltao.github.io/blog/2020/eisel-lemire.html the // primary goal of the Eisel-Lemire algorithm is speed, for 99+% of the cases, // not 100% coverage. As long as Eisel-Lemire doesn’t claim false positives, // the combined approach (falling back to an alternative implementation when // this function returns false) is both fast and correct. template <typename FloatType> bool EiselLemire(const strings_internal::ParsedFloat& input, bool negative,
absl::Nonnull<FloatType*> value,
absl::Nonnull<std::errc*> ec) {
uint64_t man = input.mantissa; int exp10 = input.exponent; if (exp10 < FloatTraits<FloatType>::kEiselLemireMinInclusiveExp10) {
*value = negative ? -0.0 : 0.0;
*ec = std::errc::result_out_of_range; returntrue;
} elseif (exp10 >= FloatTraits<FloatType>::kEiselLemireMaxExclusiveExp10) { // Return max (a finite value) consistent with from_chars and DR 3081. For // SimpleAtod and SimpleAtof, post-processing will return infinity.
*value = negative ? -std::numeric_limits<FloatType>::max()
: std::numeric_limits<FloatType>::max();
*ec = std::errc::result_out_of_range; returntrue;
}
// Assert kPower10TableMinInclusive <= exp10 < kPower10TableMaxExclusive. // Equivalently, !Power10Underflow(exp10) and !Power10Overflow(exp10).
static_assert(
FloatTraits<FloatType>::kEiselLemireMinInclusiveExp10 >=
kPower10TableMinInclusive, "(exp10-kPower10TableMinInclusive) in kPower10MantissaHighTable bounds");
static_assert(
FloatTraits<FloatType>::kEiselLemireMaxExclusiveExp10 <=
kPower10TableMaxExclusive, "(exp10-kPower10TableMinInclusive) in kPower10MantissaHighTable bounds");
// The terse (+) comments in this function body refer to sections of the // https://nigeltao.github.io/blog/2020/eisel-lemire.html blog post. // // That blog post discusses double precision (11 exponent bits with a -1023 // bias, 52 mantissa bits), but the same approach applies to single precision // (8 exponent bits with a -127 bias, 23 mantissa bits). Either way, the // computation here happens with 64-bit values (e.g. man) or 128-bit values // (e.g. x) before finally converting to 64- or 32-bit floating point. // // See also "Number Parsing at a Gigabyte per Second, Software: Practice and // Experience 51 (8), 2021" (https://arxiv.org/abs/2101.11408) for detail.
// (+) Normalization. int clz = countl_zero(man);
man <<= static_cast<unsignedint>(clz); // The 217706 etc magic numbers are from the Power10Exponent function.
uint64_t ret_exp2 = static_cast<uint64_t>((217706 * exp10 >> 16) + 64 +
FloatTraits<FloatType>::kExponentBias - clz);
// (+) Wider Approximation. static constexpr uint64_t high64_mask =
FloatTraits<FloatType>::kEiselLemireMask; if (((Uint128High64(x) & high64_mask) == high64_mask) &&
(man > (std::numeric_limits<uint64_t>::max() - Uint128Low64(x)))) {
uint128 y = static_cast<uint128>(man) * static_cast<uint128>(
kPower10MantissaLowTable[exp10 - kPower10TableMinInclusive]);
x += Uint128High64(y); // For example, parsing "4503599627370497.5" will take the if-true // branch here (for double precision), since: // - x = 0x8000000000000BFF_FFFFFFFFFFFFFFFF // - y = 0x8000000000000BFF_7FFFFFFFFFFFF400 // - man = 0xA000000000000F00 // Likewise, when parsing "0.0625" for single precision: // - x = 0x7FFFFFFFFFFFFFFF_FFFFFFFFFFFFFFFF // - y = 0x813FFFFFFFFFFFFF_8A00000000000000 // - man = 0x9C40000000000000 if (((Uint128High64(x) & high64_mask) == high64_mask) &&
((Uint128Low64(x) + 1) == 0) &&
(man > (std::numeric_limits<uint64_t>::max() - Uint128Low64(y)))) { returnfalse;
}
}
// (+) Shifting to 54 Bits (or for single precision, to 25 bits).
uint64_t msb = Uint128High64(x) >> 63;
uint64_t ret_man =
Uint128High64(x) >> (msb + FloatTraits<FloatType>::kEiselLemireShift);
ret_exp2 -= 1 ^ msb;
// (+) Half-way Ambiguity. // // For example, parsing "1e+23" will take the if-true branch here (for double // precision), since: // - x = 0x54B40B1F852BDA00_0000000000000000 // - ret_man = 0x002A5A058FC295ED // Likewise, when parsing "20040229.0" for single precision: // - x = 0x4C72894000000000_0000000000000000 // - ret_man = 0x000000000131CA25 if ((Uint128Low64(x) == 0) && ((Uint128High64(x) & high64_mask) == 0) &&
((ret_man & 3) == 1)) { returnfalse;
}
// (+) From 54 to 53 Bits (or for single precision, from 25 to 24 bits).
ret_man += ret_man & 1; // Line From54a.
ret_man >>= 1; // Line From54b. // Incrementing ret_man (at line From54a) may have overflowed 54 bits (53 // bits after the right shift by 1 at line From54b), so adjust for that. // // For example, parsing "9223372036854775807" will take the if-true branch // here (for double precision), since: // - ret_man = 0x0020000000000000 = (1 << 53) // Likewise, when parsing "2147483647.0" for single precision: // - ret_man = 0x0000000001000000 = (1 << 24) if ((ret_man >> FloatTraits<FloatType>::kTargetMantissaBits) > 0) {
ret_exp2 += 1; // Conceptually, we need a "ret_man >>= 1" in this if-block to balance // incrementing ret_exp2 in the line immediately above. However, we only // get here when line From54a overflowed (after adding a 1), so ret_man // here is (1 << 53). Its low 53 bits are therefore all zeroes. The only // remaining use of ret_man is to mask it with ((1 << 52) - 1), so only its // low 52 bits matter. A "ret_man >>= 1" would have no effect in practice. // // We omit the "ret_man >>= 1", even if it is cheap (and this if-branch is // rarely taken) and technically 'more correct', so that mutation tests // that would otherwise modify or omit that "ret_man >>= 1" don't complain // that such code mutations have no observable effect.
}
// ret_exp2 is a uint64_t. Zero or underflow means that we're in subnormal // space. max_exp2 (0x7FF for double precision, 0xFF for single precision) or // above means that we're in Inf/NaN space. // // The if block is equivalent to (but has fewer branches than): // if ((ret_exp2 <= 0) || (ret_exp2 >= max_exp2)) { etc } // // For example, parsing "4.9406564584124654e-324" will take the if-true // branch here, since ret_exp2 = -51. static constexpr uint64_t max_exp2 =
(1 << FloatTraits<FloatType>::kTargetExponentBits) - 1; if ((ret_exp2 - 1) >= (max_exp2 - 1)) { returnfalse;
}
bool negative = false; if (first != last && *first == '-') {
++first;
negative = true;
} // If the `hex` flag is *not* set, then we will accept a 0x prefix and try // to parse a hexadecimal float. if ((fmt_flags & chars_format::hex) == chars_format{} && last - first >= 2 &&
*first == '0' && (first[1] == 'x' || first[1] == 'X')) { constchar* hex_first = first + 2;
strings_internal::ParsedFloat hex_parse =
strings_internal::ParseFloat<16>(hex_first, last, fmt_flags); if (hex_parse.end == nullptr ||
hex_parse.type != strings_internal::FloatType::kNumber) { // Either we failed to parse a hex float after the "0x", or we read // "0xinf" or "0xnan" which we don't want to match. // // However, a string that begins with "0x" also begins with "0", which // is normally a valid match for the number zero. So we want these // strings to match zero unless fmt_flags is `scientific`. (This flag // means an exponent is required, which the string "0" does not have.) if (fmt_flags == chars_format::scientific) {
result.ec = std::errc::invalid_argument;
} else {
result.ptr = first + 1;
value = negative ? -0.0 : 0.0;
} return result;
} // We matched a value.
result.ptr = hex_parse.end; if (HandleEdgeCase(hex_parse, negative, &value)) { return result;
}
CalculatedFloat calculated =
CalculateFromParsedHexadecimal<FloatType>(hex_parse);
EncodeResult(calculated, negative, &result, &value); return result;
} // Otherwise, we choose the number base based on the flags. if ((fmt_flags & chars_format::hex) == chars_format::hex) {
strings_internal::ParsedFloat hex_parse =
strings_internal::ParseFloat<16>(first, last, fmt_flags); if (hex_parse.end == nullptr) {
result.ec = std::errc::invalid_argument; return result;
}
result.ptr = hex_parse.end; if (HandleEdgeCase(hex_parse, negative, &value)) { return result;
}
CalculatedFloat calculated =
CalculateFromParsedHexadecimal<FloatType>(hex_parse);
EncodeResult(calculated, negative, &result, &value); return result;
} else {
strings_internal::ParsedFloat decimal_parse =
strings_internal::ParseFloat<10>(first, last, fmt_flags); if (decimal_parse.end == nullptr) {
result.ec = std::errc::invalid_argument; return result;
}
result.ptr = decimal_parse.end; if (HandleEdgeCase(decimal_parse, negative, &value)) { return result;
} // A nullptr subrange_begin means that the decimal_parse.mantissa is exact // (not truncated), a precondition of the Eisel-Lemire algorithm. if ((decimal_parse.subrange_begin == nullptr) &&
EiselLemire<FloatType>(decimal_parse, negative, &value, &result.ec)) { return result;
}
CalculatedFloat calculated =
CalculateFromParsedDecimal<FloatType>(decimal_parse);
EncodeResult(calculated, negative, &result, &value); return result;
}
}
} // namespace
from_chars_result from_chars(absl::Nonnull<constchar*> first,
absl::Nonnull<constchar*> last, double& value,
chars_format fmt) { return FromCharsImpl(first, last, value, fmt);
}
from_chars_result from_chars(absl::Nonnull<constchar*> first,
absl::Nonnull<constchar*> last, float& value,
chars_format fmt) { return FromCharsImpl(first, last, value, fmt);
}
namespace {
// Table of powers of 10, from kPower10TableMinInclusive to // kPower10TableMaxExclusive. // // kPower10MantissaHighTable[i - kPower10TableMinInclusive] stores the 64-bit // mantissa. The high bit is always on. // // kPower10MantissaLowTable extends that 64-bit mantissa to 128 bits. // // Power10Exponent(i) calculates the power-of-two exponent. // // For a number i, this gives the unique mantissaHigh and exponent such that // (mantissaHigh * 2**exponent) <= 10**i < ((mantissaHigh + 1) * 2**exponent). // // For example, Python can confirm that the exact hexadecimal value of 1e60 is: // >>> a = 1000000000000000000000000000000000000000000000000000000000000 // >>> hex(a) // '0x9f4f2726179a224501d762422c946590d91000000000000000' // Adding underscores at every 8th hex digit shows 50 hex digits: // '0x9f4f2726_179a2245_01d76242_2c946590_d9100000_00000000_00'. // In this case, the high bit of the first hex digit, 9, is coincidentally set, // so we do not have to do further shifting to deduce the 128-bit mantissa: // - kPower10MantissaHighTable[60 - kP10TMI] = 0x9f4f2726179a2245U // - kPower10MantissaLowTable[ 60 - kP10TMI] = 0x01d762422c946590U // where kP10TMI is kPower10TableMinInclusive. The low 18 of those 50 hex // digits are truncated. // // 50 hex digits (with the high bit set) is 200 bits and mantissaHigh holds 64 // bits, so Power10Exponent(60) = 200 - 64 = 136. Again, Python can confirm: // >>> b = 0x9f4f2726179a2245 // >>> ((b+0)<<136) <= a // True // >>> ((b+1)<<136) <= a // False // // The tables were generated by // https://github.com/google/wuffs/blob/315b2e52625ebd7b02d8fac13e3cd85ea374fb80/script/print-mpb-powers-of-10.go // after re-formatting its output into two arrays of N uint64_t values (instead // of an N element array of uint64_t pairs).
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung und die Messung sind noch experimentell.
