// 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.

#include "absl/strings/charconv.h"

#include <algorithm>

#include <cassert>

#include <cstddef>

#include <cstdint>

#include <limits>

#include <system_error>  // NOLINT(build/c++11)

#include "absl/base/casts.h"

#include "absl/base/config.h"

#include "absl/base/nullability.h"

#include "absl/numeric/bits.h"

#include "absl/numeric/int128.h"

#include "absl/strings/internal/charconv_bigint.h"

#include "absl/strings/internal/charconv_parse.h"

// 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

#elif defined(__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;

  static double MakeNan(const char* absl_nonnull 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.

    using namespace 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.

  static double 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.

    using namespace 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;

  static float MakeNan(const char* absl_nonnull 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.

    using namespace std;  // NOLINT

    return std::nanf(tagp);

#endif

  }

  static float 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.

    using namespace 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.

extern const uint64_t kPower10MantissaHighTable[];  // High 64 of 128 bits.

extern const 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;

uint64_t Power10Mantissa(int n) {

  return kPower10MantissaHighTable[n - kPower10TableMinInclusive];

}

int Power10Exponent(int n) {

  // The 217706 etc magic numbers encode the results as a formula instead of a

  // table. Their equivalence (over the kPower10TableMinInclusive ..

  // kPower10TableMaxExclusive range) is confirmed by

  // https://github.com/google/wuffs/blob/315b2e52625ebd7b02d8fac13e3cd85ea374fb80/script/print-mpb-powers-of-10.go

  return (217706 * n >> 16) - 63;

}

// 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.

    return static_cast<int>(bit_width(Uint128Low64(value)));

  }

  return 128 - 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) {

  const int normal_shift =

      mantissa_width - FloatTraits<FloatType>::kTargetMantissaBits;

  const int minimum_shift =

      FloatTraits<FloatType>::kMinNormalExponent - binary_exponent;

  return std::max(normal_shift, minimum_shift);

}

charconv.cc:int absl::(anonymous namespace)::NormalizedShiftSize<double>(int, int)

Line

Count

Source

335

552k

int NormalizedShiftSize(int mantissa_width, int binary_exponent) {

336

552k

  const int normal_shift =

337

552k

      mantissa_width - FloatTraits<FloatType>::kTargetMantissaBits;

338

552k

  const int minimum_shift =

339

552k

      FloatTraits<FloatType>::kMinNormalExponent - binary_exponent;

340

552k

  return std::max(normal_shift, minimum_shift);

341

552k

}

Unexecuted instantiation: charconv.cc:int absl::(anonymous namespace)::NormalizedShiftSize<float>(int, int)

// 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, uint128* absl_nonnull value) {

  const int current_bit_width = BitWidth(*value);

  const int 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,

                    FloatType* absl_nonnull value) {

  if (input.type == strings_internal::FloatType::kNan) {

    // A bug in 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://gcc.gnu.org/bugzilla/show_bug.cgi?id=86113

    constexpr ptrdiff_t kNanBufferSize = 128;

#if (defined(__GNUC__) && !defined(__clang__))

    volatile char 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);

    return true;

  }

  if (input.type == strings_internal::FloatType::kInfinity) {

    *value = negative ? -std::numeric_limits<FloatType>::infinity()

                      : std::numeric_limits<FloatType>::infinity();

    return true;

  }

  if (input.mantissa == 0) {

    *value = negative ? -0.0f : 0.0f;

    return true;

  }

  return false;

}

charconv.cc:bool absl::(anonymous namespace)::HandleEdgeCase<double>(absl::strings_internal::ParsedFloat const&, bool, double*)

Line

Count

Source

360

3.97M

                    FloatType* absl_nonnull value) {

361

3.97M

  if (input.type == strings_internal::FloatType::kNan) {

362

    // A bug in gcc would cause the compiler to optimize away the buffer we are

363

    // building below.  Declaring the buffer volatile avoids the issue, and has

364

    // no measurable performance impact in microbenchmarks.

365

    //

366

    // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=86113

367

1.36k

    constexpr ptrdiff_t kNanBufferSize = 128;

368

#if (defined(__GNUC__) && !defined(__clang__))

369

    volatile char n_char_sequence[kNanBufferSize];

370

#else

371

1.36k

    char n_char_sequence[kNanBufferSize];

372

1.36k

#endif

373

1.36k

    if (input.subrange_begin == nullptr) {

374

476

      n_char_sequence[0] = '\0';

375

889

    } else {

376

889

      ptrdiff_t nan_size = input.subrange_end - input.subrange_begin;

377

889

      nan_size = std::min(nan_size, kNanBufferSize - 1);

378

889

      std::copy_n(input.subrange_begin, nan_size, n_char_sequence);

379

889

      n_char_sequence[nan_size] = '\0';

380

889

    }

381

1.36k

    char* nan_argument = const_cast<char*>(n_char_sequence);

382

1.36k

    *value = negative ? -FloatTraits<FloatType>::MakeNan(nan_argument)

383

1.36k

                      : FloatTraits<FloatType>::MakeNan(nan_argument);

384

1.36k

    return true;

385

1.36k

  }

386

3.97M

  if (input.type == strings_internal::FloatType::kInfinity) {

387

145

    *value = negative ? -std::numeric_limits<FloatType>::infinity()

388

145

                      : std::numeric_limits<FloatType>::infinity();

389

145

    return true;

390

145

  }

391

3.97M

  if (input.mantissa == 0) {

392

177k

    *value = negative ? -0.0f : 0.0f;

393

177k

    return true;

394

177k

  }

395

3.80M

  return false;

396

3.97M

}

Unexecuted instantiation: charconv.cc:bool absl::(anonymous namespace)::HandleEdgeCase<float>(absl::strings_internal::ParsedFloat const&, bool, float*)

// 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::from_chars_result* absl_nonnull result,

                  FloatType* absl_nonnull 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;

  } else if (calculated.mantissa == 0 || calculated.exponent == kUnderflow) {

    result->ec = std::errc::result_out_of_range;

    *value = negative ? -0.0f : 0.0f;

    return;

  }

  *value = FloatTraits<FloatType>::Make(

      static_cast<typename FloatTraits<FloatType>::mantissa_t>(

          calculated.mantissa),

      calculated.exponent, negative);

}

charconv.cc:void absl::(anonymous namespace)::EncodeResult<double>(absl::(anonymous namespace)::CalculatedFloat const&, bool, absl::from_chars_result*, double*)

Line

Count

Source

407

553k

                  FloatType* absl_nonnull value) {

408

553k

  if (calculated.exponent == kOverflow) {

409

678

    result->ec = std::errc::result_out_of_range;

410

678

    *value = negative ? -std::numeric_limits<FloatType>::max()

411

678

                      : std::numeric_limits<FloatType>::max();

412

678

    return;

413

553k

  } else if (calculated.mantissa == 0 || calculated.exponent == kUnderflow) {

414

1.78k

    result->ec = std::errc::result_out_of_range;

415

1.78k

    *value = negative ? -0.0f : 0.0f;

416

1.78k

    return;

417

1.78k

  }

418

551k

  *value = FloatTraits<FloatType>::Make(

419

551k

      static_cast<typename FloatTraits<FloatType>::mantissa_t>(

420

551k

          calculated.mantissa),

421

551k

      calculated.exponent, negative);

422

551k

}

Unexecuted instantiation: charconv.cc:void absl::(anonymous namespace)::EncodeResult<float>(absl::(anonymous namespace)::CalculatedFloat const&, bool, absl::from_chars_result*, float*)

// 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,

                            bool* absl_nonnull output_exact) {

  if (shift <= 0) {

    *output_exact = input_exact;

    return static_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;

    return 0;

  }

  *output_exact = true;

  const uint128 shift_mask = (uint128(1) << shift) - 1;

  const uint128 halfway_point = uint128(1) << (shift - 1);

  const uint128 shifted_bits = value & shift_mask;

  value >>= shift;

  if (shifted_bits > halfway_point) {

    // Shifted bits greater than 10000... require rounding up.

    return static_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;

    }

    return static_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.

  return static_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) {

    return false;

  } else if (comparison > 0) {

    return true;

  } 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;

  } else if (mantissa == 0) {

    result.exponent = kUnderflow;

  } else {

    result.exponent = exponent;

    result.mantissa = mantissa;

  }

  return result;

}

charconv.cc:absl::(anonymous namespace)::CalculatedFloat absl::(anonymous namespace)::CalculatedFloatFromRawValues<double>(unsigned long, int)

Line

Count

Source

580

552k

CalculatedFloat CalculatedFloatFromRawValues(uint64_t mantissa, int exponent) {

581

552k

  CalculatedFloat result;

582

552k

  if (mantissa == uint64_t{1} << FloatTraits<FloatType>::kTargetMantissaBits) {

583

1.28k

    mantissa >>= 1;

584

1.28k

    exponent += 1;

585

1.28k

  }

586

552k

  if (exponent > FloatTraits<FloatType>::kMaxExponent) {

587

632

    result.exponent = kOverflow;

588

552k

  } else if (mantissa == 0) {

589

932

    result.exponent = kUnderflow;

590

551k

  } else {

591

551k

    result.exponent = exponent;

592

551k

    result.mantissa = mantissa;

593

551k

  }

594

552k

  return result;

595

552k

}

Unexecuted instantiation: charconv.cc:absl::(anonymous namespace)::CalculatedFloat absl::(anonymous namespace)::CalculatedFloatFromRawValues<float>(unsigned long, int)

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));

  const int 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);

}

charconv.cc:absl::(anonymous namespace)::CalculatedFloat absl::(anonymous namespace)::CalculateFromParsedHexadecimal<double>(absl::strings_internal::ParsedFloat const&)

Line

Count

Source

599

13.9k

    const strings_internal::ParsedFloat& parsed_hex) {

600

13.9k

  uint64_t mantissa = parsed_hex.mantissa;

601

13.9k

  int exponent = parsed_hex.exponent;

602

  // This static_cast is only needed when using a std::bit_width()

603

  // implementation that does not have the fix for LWG 3656 applied.

604

13.9k

  int mantissa_width = static_cast<int>(bit_width(mantissa));

605

13.9k

  const int shift = NormalizedShiftSize<FloatType>(mantissa_width, exponent);

606

13.9k

  bool result_exact;

607

13.9k

  exponent += shift;

608

13.9k

  mantissa = ShiftRightAndRound(mantissa, shift,

609

13.9k

                                /* input exact= */ true, &result_exact);

610

  // ParseFloat handles rounding in the hexadecimal case, so we don't have to

611

  // check `result_exact` here.

612

13.9k

  return CalculatedFloatFromRawValues<FloatType>(mantissa, exponent);

613

13.9k

}

Unexecuted instantiation: charconv.cc:absl::(anonymous namespace)::CalculatedFloat absl::(anonymous namespace)::CalculateFromParsedHexadecimal<float>(absl::strings_internal::ParsedFloat const&)

template <typename FloatType>

CalculatedFloat CalculateFromParsedDecimal(

    const strings_internal::ParsedFloat& parsed_decimal) {

  CalculatedFloat result;

  // Large or small enough decimal exponents will always result in overflow

  // or underflow.

  if (Power10Underflow(parsed_decimal.exponent)) {

    result.exponent = kUnderflow;

    return result;

  } else if (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);

  } else if (!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.

  const int 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;

    }

  }

  return CalculatedFloatFromRawValues<FloatType>(binary_mantissa,

                                                 binary_exponent);

}

charconv.cc:absl::(anonymous namespace)::CalculatedFloat absl::(anonymous namespace)::CalculateFromParsedDecimal<double>(absl::strings_internal::ParsedFloat const&)

Line

Count

Source

617

539k

    const strings_internal::ParsedFloat& parsed_decimal) {

618

539k

  CalculatedFloat result;

619

620

  // Large or small enough decimal exponents will always result in overflow

621

  // or underflow.

622

539k

  if (Power10Underflow(parsed_decimal.exponent)) {

623

854

    result.exponent = kUnderflow;

624

854

    return result;

625

539k

  } else if (Power10Overflow(parsed_decimal.exponent)) {

626

46

    result.exponent = kOverflow;

627

46

    return result;

628

46

  }

629

630

  // Otherwise convert our power of 10 into a power of 2 times an integer

631

  // mantissa, and multiply this by our parsed decimal mantissa.

632

539k

  uint128 wide_binary_mantissa = parsed_decimal.mantissa;

633

539k

  wide_binary_mantissa *= Power10Mantissa(parsed_decimal.exponent);

634

539k

  int binary_exponent = Power10Exponent(parsed_decimal.exponent);

635

636

  // Discard bits that are inaccurate due to truncation error.  The magic

637

  // `mantissa_width` constants below are justified in

638

  // https://abseil.io/about/design/charconv. They represent the number of bits

639

  // in `wide_binary_mantissa` that are guaranteed to be unaffected by error

640

  // propagation.

641

539k

  bool mantissa_exact;

642

539k

  int mantissa_width;

643

539k

  if (parsed_decimal.subrange_begin) {

644

    // Truncated mantissa

645

31.6k

    mantissa_width = 58;

646

31.6k

    mantissa_exact = false;

647

31.6k

    binary_exponent +=

648

31.6k

        TruncateToBitWidth(mantissa_width, &wide_binary_mantissa);

649

507k

  } else if (!Power10Exact(parsed_decimal.exponent)) {

650

    // Exact mantissa, truncated power of ten

651

4.96k

    mantissa_width = 63;

652

4.96k

    mantissa_exact = false;

653

4.96k

    binary_exponent +=

654

4.96k

        TruncateToBitWidth(mantissa_width, &wide_binary_mantissa);

655

502k

  } else {

656

    // Product is exact

657

502k

    mantissa_width = BitWidth(wide_binary_mantissa);

658

502k

    mantissa_exact = true;

659

502k

  }

660

661

  // Shift into an FloatType-sized mantissa, and round to nearest.

662

539k

  const int shift =

663

539k

      NormalizedShiftSize<FloatType>(mantissa_width, binary_exponent);

664

539k

  bool result_exact;

665

539k

  binary_exponent += shift;

666

539k

  uint64_t binary_mantissa = ShiftRightAndRound(wide_binary_mantissa, shift,

667

539k

                                                mantissa_exact, &result_exact);

668

539k

  if (!result_exact) {

669

    // We could not determine the rounding direction using int128 math.  Use

670

    // full resolution math instead.

671

11.5k

    if (MustRoundUp(binary_mantissa, binary_exponent, parsed_decimal)) {

672

5.60k

      binary_mantissa += 1;

673

5.60k

    }

674

11.5k

  }

675

676

539k

  return CalculatedFloatFromRawValues<FloatType>(binary_mantissa,

677

539k

                                                 binary_exponent);

678

539k

}

Unexecuted instantiation: charconv.cc:absl::(anonymous namespace)::CalculatedFloat absl::(anonymous namespace)::CalculateFromParsedDecimal<float>(absl::strings_internal::ParsedFloat const&)

// 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,

                 FloatType* absl_nonnull value, std::errc* absl_nonnull ec) {

  uint64_t man = input.mantissa;

  int exp10 = input.exponent;

  if (exp10 < FloatTraits<FloatType>::kEiselLemireMinInclusiveExp10) {

    *value = negative ? -0.0f : 0.0f;

    *ec = std::errc::result_out_of_range;

    return true;

  } else if (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;

    return true;

  }

  // 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<unsigned int>(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);

  // (+) Multiplication.

  uint128 x = static_cast<uint128>(man) *

              static_cast<uint128>(

                  kPower10MantissaHighTable[exp10 - kPower10TableMinInclusive]);

  // (+) 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)))) {

      return false;

    }

  }

  // (+) 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)) {

    return false;

  }

  // (+) 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)) {

    return false;

  }

#ifndef ABSL_BIT_PACK_FLOATS

  if (FloatTraits<FloatType>::kTargetBits == 64) {

    *value = FloatTraits<FloatType>::Make(

        (ret_man & 0x000FFFFFFFFFFFFFu) | 0x0010000000000000u,

        static_cast<int>(ret_exp2) - 1023 - 52, negative);

    return true;

  } else if (FloatTraits<FloatType>::kTargetBits == 32) {

    *value = FloatTraits<FloatType>::Make(

        (static_cast<uint32_t>(ret_man) & 0x007FFFFFu) | 0x00800000u,

        static_cast<int>(ret_exp2) - 127 - 23, negative);

    return true;

  }

#else

  if (FloatTraits<FloatType>::kTargetBits == 64) {

    uint64_t ret_bits = (ret_exp2 << 52) | (ret_man & 0x000FFFFFFFFFFFFFu);

    if (negative) {

      ret_bits |= 0x8000000000000000u;

    }

    *value = static_cast<FloatType>(absl::bit_cast<double>(ret_bits));

    return true;

  } else if (FloatTraits<FloatType>::kTargetBits == 32) {

    uint32_t ret_bits = (static_cast<uint32_t>(ret_exp2) << 23) |

                        (static_cast<uint32_t>(ret_man) & 0x007FFFFFu);

    if (negative) {

      ret_bits |= 0x80000000u;

    }

    *value = static_cast<FloatType>(absl::bit_cast<float>(ret_bits));

    return true;

  }

#endif  // ABSL_BIT_PACK_FLOATS

  return false;

}

charconv.cc:bool absl::(anonymous namespace)::EiselLemire<double>(absl::strings_internal::ParsedFloat const&, bool, double*, std::__1::errc*)

Line

Count

Source

687

3.75M

                 FloatType* absl_nonnull value, std::errc* absl_nonnull ec) {

688

3.75M

  uint64_t man = input.mantissa;

689

3.75M

  int exp10 = input.exponent;

690

3.75M

  if (exp10 < FloatTraits<FloatType>::kEiselLemireMinInclusiveExp10) {

691

377

    *value = negative ? -0.0f : 0.0f;

692

377

    *ec = std::errc::result_out_of_range;

693

377

    return true;

694

3.75M

  } else if (exp10 >= FloatTraits<FloatType>::kEiselLemireMaxExclusiveExp10) {

695

    // Return max (a finite value) consistent with from_chars and DR 3081. For

696

    // SimpleAtod and SimpleAtof, post-processing will return infinity.

697

504

    *value = negative ? -std::numeric_limits<FloatType>::max()

698

504

                      : std::numeric_limits<FloatType>::max();

699

504

    *ec = std::errc::result_out_of_range;

700

504

    return true;

701

504

  }

702

703

  // Assert kPower10TableMinInclusive <= exp10 < kPower10TableMaxExclusive.

704

  // Equivalently, !Power10Underflow(exp10) and !Power10Overflow(exp10).

705

3.75M

  static_assert(

706

3.75M

      FloatTraits<FloatType>::kEiselLemireMinInclusiveExp10 >=

707

3.75M

          kPower10TableMinInclusive,

708

3.75M

      "(exp10-kPower10TableMinInclusive) in kPower10MantissaHighTable bounds");

709

3.75M

  static_assert(

710

3.75M

      FloatTraits<FloatType>::kEiselLemireMaxExclusiveExp10 <=

711

3.75M

          kPower10TableMaxExclusive,

712

3.75M

      "(exp10-kPower10TableMinInclusive) in kPower10MantissaHighTable bounds");

713

714

  // The terse (+) comments in this function body refer to sections of the

715

  // https://nigeltao.github.io/blog/2020/eisel-lemire.html blog post.

716

  //

717

  // That blog post discusses double precision (11 exponent bits with a -1023

718

  // bias, 52 mantissa bits), but the same approach applies to single precision

719

  // (8 exponent bits with a -127 bias, 23 mantissa bits). Either way, the

720

  // computation here happens with 64-bit values (e.g. man) or 128-bit values

721

  // (e.g. x) before finally converting to 64- or 32-bit floating point.

722

  //

723

  // See also "Number Parsing at a Gigabyte per Second, Software: Practice and

724

  // Experience 51 (8), 2021" (https://arxiv.org/abs/2101.11408) for detail.

725

726

  // (+) Normalization.

727

3.75M

  int clz = countl_zero(man);