[Bug libstdc++/107468] New: std::from_chars doesn't always round to nearest

["jakub at gcc dot gnu.org" <gcc-bugzilla@gcc.gnu.org>]

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

            Bug ID: 107468
           Summary: std::from_chars doesn't always round to nearest
           Product: gcc
           Version: 13.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: libstdc++
          Assignee: unassigned at gcc dot gnu.org
          Reporter: jakub at gcc dot gnu.org
  Target Milestone: ---

#include <cfenv>
#include <charconv>
#include <iostream>
#include <string_view>

int
main()
{
  float f;
  char buf[] = "3.355447e+07";
  fesetround(FE_DOWNWARD);
  auto [ptr, ec] = std::from_chars(buf, buf + sizeof (buf) - 1, f,
std::chars_format::scientific);
  char buf2[128];
  auto [ptr2, ec2] = std::to_chars(buf2, buf2 + 128, f,
std::chars_format::fixed);
  std::cout << std::string_view(buf2, ptr2 - buf2) << std::endl;
}

33554470 isn't representable in IEEE single, only 33554468 and 33554472 are,
the former is 0x1.000012p+25 and the latter is 0x1.000014p+25 and the latter
has
0 as the least significant digit.  So, round to even should be 33554472 and
that is what happens without the fesetround(FE_DOWNWARD).

When fast_float isn't used, floating_from_chars.cc temporarily overrides the
rounding to nearest, but when it is used, it doesn't.
In most cases it is fine, because fast_float usually computes the mantissa and
exponent in integral code.  But it has two exceptions for this I think.
One is the fast path:
  // Next is Clinger's fast path.
  if (binary_format<T>::min_exponent_fast_path() <= pns.exponent &&
pns.exponent <= binary_format<T>::max_exponent_fast_path() && pns.mantissa
<=binary_format<T>::max_mantissa_fast_p
ath() && !pns.too_many_digits) {
    value = T(pns.mantissa);
    if (pns.exponent < 0) { value = value /
binary_format<T>::exact_power_of_ten(-pns.exponent); }
    else { value = value * binary_format<T>::exact_power_of_ten(pns.exponent);
}
    if (pns.negative) { value = -value; }
    return answer;
  }
I'm afraid we need to temporarily override rounding to nearest around the
multiplications in there.  And the other one is:
  T b;
  to_float(false, am_b, b);
  adjusted_mantissa theor = to_extended_halfway(b);
where I think we don't really need it but I don't understand the point of
jumping through the floating point format, it encodes am_b into an integer,
memcpys it to float/double (in to_float), then immediately memcpys it back to
an integer (in to_extended called from to_extended_halfway) and then unpacks it
with some adjustments.