in include/ylt/util/dragonbox.h [2170:2318]
JKJ_SAFEBUFFERS static ReturnType compute_nearest_normal(
carrier_uint const two_fc, int const exponent,
AdditionalArgs... additional_args) noexcept {
//////////////////////////////////////////////////////////////////////
// Step 1: Schubfach multiplier calculation
//////////////////////////////////////////////////////////////////////
ReturnType ret_value;
IntervalType interval_type{additional_args...};
// Compute k and beta.
int const minus_k = log::floor_log10_pow2(exponent) - kappa;
auto const cache = CachePolicy::template get_cache<format>(-minus_k);
int const beta = exponent + log::floor_log2_pow10(-minus_k);
// Compute zi and deltai.
// 10^kappa <= deltai < 10^(kappa + 1)
auto const deltai = compute_delta(cache, beta);
// For the case of binary32, the result of integer check is not correct for
// 29711844 * 2^-82
// = 6.1442653300000000008655037797566933477355632930994033813476... *
// 10^-18 and 29711844 * 2^-81
// = 1.2288530660000000001731007559513386695471126586198806762695... *
// 10^-17, and they are the unique counterexamples. However, since 29711844
// is even, this does not cause any problem for the endpoints calculations;
// it can only cause a problem when we need to perform integer check for the
// center. Fortunately, with these inputs, that branch is never executed, so
// we are fine.
auto const [zi, is_z_integer] = compute_mul((two_fc | 1) << beta, cache);
//////////////////////////////////////////////////////////////////////
// Step 2: Try larger divisor; remove trailing zeros if necessary
//////////////////////////////////////////////////////////////////////
constexpr auto big_divisor = compute_power<kappa + 1>(std::uint32_t(10));
constexpr auto small_divisor = compute_power<kappa>(std::uint32_t(10));
// Using an upper bound on zi, we might be able to optimize the division
// better than the compiler; we are computing zi / big_divisor here.
ret_value.significand = div::divide_by_pow10<
kappa + 1, carrier_uint,
(carrier_uint(1) << (significand_bits + 1)) * big_divisor - 1>(zi);
auto r = std::uint32_t(zi - big_divisor * ret_value.significand);
if (r < deltai) {
// Exclude the right endpoint if necessary.
if (r == 0 && (is_z_integer & !interval_type.include_right_endpoint())) {
if constexpr (BinaryToDecimalRoundingPolicy::tag ==
policy_impl::binary_to_decimal_rounding::tag_t::
do_not_care) {
ret_value.significand *= 10;
ret_value.exponent = minus_k + kappa;
--ret_value.significand;
TrailingZeroPolicy::template no_trailing_zeros<impl>(ret_value);
return ret_value;
}
else {
--ret_value.significand;
r = big_divisor;
goto small_divisor_case_label;
}
}
}
else if (r > deltai) {
goto small_divisor_case_label;
}
else {
// r == deltai; compare fractional parts.
auto const [xi_parity, x_is_integer] =
compute_mul_parity(two_fc - 1, cache, beta);
if (!(xi_parity |
(x_is_integer & interval_type.include_left_endpoint()))) {
goto small_divisor_case_label;
}
}
ret_value.exponent = minus_k + kappa + 1;
// We may need to remove trailing zeros.
TrailingZeroPolicy::template on_trailing_zeros<impl>(ret_value);
return ret_value;
//////////////////////////////////////////////////////////////////////
// Step 3: Find the significand with the smaller divisor
//////////////////////////////////////////////////////////////////////
small_divisor_case_label:
TrailingZeroPolicy::template no_trailing_zeros<impl>(ret_value);
ret_value.significand *= 10;
ret_value.exponent = minus_k + kappa;
if constexpr (BinaryToDecimalRoundingPolicy::tag ==
policy_impl::binary_to_decimal_rounding::tag_t::do_not_care) {
// Normally, we want to compute
// ret_value.significand += r / small_divisor
// and return, but we need to take care of the case that the resulting
// value is exactly the right endpoint, while that is not included in the
// interval.
if (!interval_type.include_right_endpoint()) {
// Is r divisible by 10^kappa?
if (is_z_integer &&
div::check_divisibility_and_divide_by_pow10<kappa>(r)) {
// This should be in the interval.
ret_value.significand += r - 1;
}
else {
ret_value.significand += r;
}
}
else {
ret_value.significand += div::small_division_by_pow10<kappa>(r);
}
}
else {
auto dist = r - (deltai / 2) + (small_divisor / 2);
bool const approx_y_parity = ((dist ^ (small_divisor / 2)) & 1) != 0;
// Is dist divisible by 10^kappa?
bool const divisible_by_small_divisor =
div::check_divisibility_and_divide_by_pow10<kappa>(dist);
// Add dist / 10^kappa to the significand.
ret_value.significand += dist;
if (divisible_by_small_divisor) {
// Check z^(f) >= epsilon^(f).
// We have either yi == zi - epsiloni or yi == (zi - epsiloni) - 1,
// where yi == zi - epsiloni if and only if z^(f) >= epsilon^(f).
// Since there are only 2 possibilities, we only need to care about the
// parity. Also, zi and r should have the same parity since the divisor
// is an even number.
auto const [yi_parity, is_y_integer] =
compute_mul_parity(two_fc, cache, beta);
if (yi_parity != approx_y_parity) {
--ret_value.significand;
}
else {
// If z^(f) >= epsilon^(f), we might have a tie
// when z^(f) == epsilon^(f), or equivalently, when y is an integer.
// For tie-to-up case, we can just choose the upper one.
if (BinaryToDecimalRoundingPolicy::prefer_round_down(ret_value) &
is_y_integer) {
--ret_value.significand;
}
}
}
}
return ret_value;
}