Fix std::float16_t rounding of subnormal ties and fast-path overflow

include/fast_float/decimal_to_binary.h

@@ -156,10 +156,20 @@ compute_float(int64_t q, uint64_t w) noexcept {
       return answer;
     }
     // next line is safe because -answer.power2 + 1 < 64
-    answer.mantissa >>= -answer.power2 + 1;
-    // Thankfully, we can't have both "round-to-even" and subnormals because
-    // "round-to-even" only occurs for powers close to 0 in the 32-bit and
-    // and 64-bit case (with no more than 19 digits).
+    int const subnormal_shift = -answer.power2 + 1;
+    answer.mantissa >>= subnormal_shift;
+    // A subnormal result can also fall exactly between two floats. With at
+    // most 19 digits this never happens for float and double, but it does for
+    // std::float16_t (e.g., 2^-25 = 298023223876953125e-25), so we apply the
+    // same round-to-even test as in the normal case below.
+    // See script/format_parameters.py.
+    if ((product.low <= 1) && (q >= binary::min_exponent_round_to_even()) &&
+        (q <= binary::max_exponent_round_to_even()) &&
+        ((answer.mantissa & 3) == 1)) {
+      if (((answer.mantissa << subnormal_shift) << shift) == product.high) {
+        answer.mantissa &= ~uint64_t(1); // flip it so that we do not round up
+      }
+    }
     answer.mantissa += (answer.mantissa & 1); // round up
     answer.mantissa >>= 1;
     // There is a weird scenario where we don't have a subnormal but just.

include/fast_float/float_common.h

@@ -225,16 +225,12 @@ using parse_options = parse_options_t<char>;
 
 #ifndef FASTFLOAT_ASSERT
 #define FASTFLOAT_ASSERT(x)                                                    \
-  {                                                                            \
-    static_cast<void>(x);                                                      \
-  }
+  { static_cast<void>(x); }
 #endif
 
 #ifndef FASTFLOAT_DEBUG_ASSERT
 #define FASTFLOAT_DEBUG_ASSERT(x)                                              \
-  {                                                                            \
-    static_cast<void>(x);                                                      \
-  }
+  { static_cast<void>(x); }
 #endif
 
 // rust style `try!()` macro, or `?` operator
@@ -665,6 +661,7 @@ template <typename T> struct binary_format : binary_format_lookup_tables<T> {
   static constexpr uint64_t max_mantissa_fast_path(int64_t power);
   static constexpr uint64_t
   max_mantissa_fast_path(); // used when fegetround() == FE_TONEAREST
+  static constexpr bool fast_path_can_overflow();
   static constexpr int largest_power_of_ten();
   static constexpr int smallest_power_of_ten();
   static constexpr T exact_power_of_ten(int64_t power);
@@ -696,6 +693,7 @@ template <typename U> struct binary_format_lookup_tables<double, U> {
       0x20000000000000 / (constant_55555 * constant_55555 * 5),
       0x20000000000000 / (constant_55555 * constant_55555 * 5 * 5),
       0x20000000000000 / (constant_55555 * constant_55555 * 5 * 5 * 5),
+      0x20000000000000 / (constant_55555 * constant_55555 * 5 * 5 * 5 * 5),
       0x20000000000000 / (constant_55555 * constant_55555 * constant_55555),
       0x20000000000000 / (constant_55555 * constant_55555 * constant_55555 * 5),
       0x20000000000000 /
@@ -711,9 +709,7 @@ template <typename U> struct binary_format_lookup_tables<double, U> {
       0x20000000000000 / (constant_55555 * constant_55555 * constant_55555 *
                           constant_55555 * 5 * 5),
       0x20000000000000 / (constant_55555 * constant_55555 * constant_55555 *
-                          constant_55555 * 5 * 5 * 5),
-      0x20000000000000 / (constant_55555 * constant_55555 * constant_55555 *
-                          constant_55555 * 5 * 5 * 5 * 5)};
+                          constant_55555 * 5 * 5 * 5)};
 };
 
 #if FASTFLOAT_DETAIL_MUST_DEFINE_CONSTEXPR_VARIABLE
@@ -930,7 +926,9 @@ binary_format<std::float16_t>::max_mantissa_fast_path(int64_t power) {
 
 template <>
 inline constexpr int binary_format<std::float16_t>::min_exponent_fast_path() {
-  return 0;
+  // w / 10^k with w <= 2^11 and k <= 4 rounds correctly even when evaluated
+  // in float or double first (checked in script/format_parameters.py).
+  return -4;
 }
 
 template <>
@@ -942,7 +940,9 @@ binary_format<std::float16_t>::max_exponent_round_to_even() {
 template <>
 inline constexpr int
 binary_format<std::float16_t>::min_exponent_round_to_even() {
-  return -22;
+  // -22 covers the normal ties; subnormal ties such as
+  // 2^-25 = 298023223876953125e-25 need q = -25 and q = -26.
+  return -26;
 }
 
 template <>
@@ -966,7 +966,8 @@ inline constexpr int binary_format<std::float16_t>::largest_power_of_ten() {
 
 template <>
 inline constexpr int binary_format<std::float16_t>::smallest_power_of_ten() {
-  return -27;
+  // (10^19 - 1) * 10^-27 < 2^-25, so any q < -26 rounds to zero.
+  return -26;
 }
 
 template <>
@@ -1053,7 +1054,8 @@ binary_format<std::bfloat16_t>::max_mantissa_fast_path(int64_t power) {
 
 template <>
 inline constexpr int binary_format<std::bfloat16_t>::min_exponent_fast_path() {
-  return 0;
+  // Same argument as for std::float16_t (w <= 2^8, k <= 3).
+  return -3;
 }
 
 template <>
@@ -1089,7 +1091,8 @@ inline constexpr int binary_format<std::bfloat16_t>::largest_power_of_ten() {
 
 template <>
 inline constexpr int binary_format<std::bfloat16_t>::smallest_power_of_ten() {
-  return -60;
+  // (10^19 - 1) * 10^-60 < 2^-134, so any q < -59 rounds to zero.
+  return -59;
 }
 
 template <>
@@ -1098,6 +1101,15 @@ inline constexpr size_t binary_format<std::bfloat16_t>::max_digits() {
 }
 #endif // __STDCPP_BFLOAT16_T__
 
+// Whether Clinger's fast path can overflow: only for std::float16_t, where
+// 2^11 * 10^4 > 65504.
+template <typename T>
+inline constexpr bool binary_format<T>::fast_path_can_overflow() {
+  return double(max_mantissa_fast_path()) *
+             double(exact_power_of_ten(max_exponent_fast_path())) >
+         double((std::numeric_limits<T>::max)());
+}
+
 template <>
 inline constexpr uint64_t
 binary_format<double>::max_mantissa_fast_path(int64_t power) {

include/fast_float/parse_number.h

@@ -221,6 +221,12 @@ clinger_fast_path_impl(uint64_t mantissa, int64_t exponent, bool is_negative,
         value = value / binary_format<T>::exact_power_of_ten(-exponent);
       } else {
         value = value * binary_format<T>::exact_power_of_ten(exponent);
+        // Only std::float16_t can overflow here (e.g., "656e2"); let the
+        // slow path report result_out_of_range.
+        if (binary_format<T>::fast_path_can_overflow() &&
+            value > (std::numeric_limits<T>::max)()) {
+          return false;
+        }
       }
       if (is_negative) {
         value = -value;

script/format_parameters.py

@@ -0,0 +1,411 @@
+#
+# Checks the binary_format<T> constants in include/fast_float/float_common.h
+# with exact integer arithmetic, for double, float, std::float16_t and
+# std::bfloat16_t. Every string is parsed as w * 10^q with w < 10^19, so each
+# constant has a condition on (w, q) that it must satisfy.
+#
+# References:
+# Daniel Lemire, Number Parsing at a Gigabyte per Second,
+# Software: Practice and Experience 51 (8), 2021 https://arxiv.org/abs/2101.11408
+# Noble Mushtak and Daniel Lemire, Fast Number Parsing Without Fallback,
+# Software: Practice and Experience 53 (6), 2023 https://arxiv.org/abs/2212.06644
+#
+# Usage: python3 script/format_parameters.py [--types float16,bfloat16]
+#        [--sweep-width 3] [--jobs N]
+#
+
+import argparse
+import multiprocessing
+import os
+import re
+import sys
+
+HEADER = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..",
+                      "include", "fast_float", "float_common.h")
+MASK64 = (1 << 64) - 1
+MAX_W = 10**19
+CXX_TYPES = {"double": "double", "float": "float",
+             "float16": "std::float16_t", "bfloat16": "std::bfloat16_t"}
+
+
+class Format:
+    def __init__(self, name):
+        self.name = name
+        cxx = CXX_TYPES[name]
+        with open(HEADER) as f:
+            text = f.read()
+        pattern = r"binary_format<" + re.escape(cxx) + r">::(\w+)\(\)\s*\{(.*?)\n\}"
+        for fn, body in re.findall(pattern, text, re.S):
+            returns = re.findall(r"return\s+([^;]+);", body)
+            # the last return is the FLT_EVAL_METHOD 0/1 branch
+            if returns and re.fullmatch(r"-?(0x[0-9A-Fa-f]+|\d+)", returns[-1].strip()):
+                setattr(self, fn, int(returns[-1], 0))
+        tables = re.search(r"binary_format_lookup_tables<" + re.escape(cxx) +
+                           r", U>\s*\{(.*?)\n\};", text, re.S).group(1)
+        powers = re.search(r"powers_of_ten\[\]\s*=\s*\{(.*?)\}", tables, re.S).group(1)
+        self.powers_of_ten = [float(x) for x in
+                              re.sub(r"(bf16|f16|f)\b", "", powers).split(",") if x.strip()]
+        mant = re.search(r"max_mantissa\[\]\s*=\s*\{(.*?)\}", tables, re.S).group(1)
+        mant = re.sub(r"//[^\n]*", "", mant).replace("/", "//").replace("\n", " ")
+        self.max_mantissa = [eval(x.strip(), {"constant_55555": 5**5})
+                             for x in mant.split(",") if x.strip()]
+        self.p = self.mantissa_explicit_bits
+        self.bias = self.p - self.minimum_exponent
+        self.denorm_min_exp = 1 - self.bias  # smallest subnormal is 2^denorm_min_exp
+        self.emin = self.minimum_exponent + 1  # leading bit of the smallest normal
+        self.emax = self.infinite_power - 1 + self.minimum_exponent
+
+    def grid(self):
+        """Every finite value and every midpoint as (num, exp2, kind) with
+        kind in value, tie_down (even lower neighbour), tie_up, overflow."""
+        out = []
+        two_p = 1 << self.p
+        d = self.denorm_min_exp
+        for m in range(1, two_p):
+            out.append((m, d, "value"))
+        for m in range(two_p):
+            out.append((2 * m + 1, d - 1, "tie_down" if m % 2 == 0 else "tie_up"))
+        for e in range(self.emin, self.emax + 1):
+            for f in range(two_p):
+                out.append((two_p + f, e - self.p, "value"))
+                kind = "tie_down" if f % 2 == 0 else "tie_up"
+                if e == self.emax and f == two_p - 1:
+                    kind = "overflow"
+                out.append((2 * (two_p + f) + 1, e - self.p - 1, kind))
+        return out
+
+
+def ratio(num, exp2):
+    return (num << exp2, 1) if exp2 >= 0 else (num, 1 << -exp2)
+
+
+def ratio10(w, q):
+    return (w * 10**q, 1) if q >= 0 else (w, 10**-q)
+
+
+def le(n1, d1, n2, d2):
+    return n1 * d2 <= n2 * d1
+
+
+def floor_log2(n, d):
+    e = n.bit_length() - d.bit_length()
+    if (d << e if e >= 0 else d) > (n if e >= 0 else n << -e):
+        e -= 1
+    return e
+
+
+def round_ties_even(n, d):
+    m, r = divmod(n, d)
+    if 2 * r > d or (2 * r == d and m & 1):
+        m += 1
+    return m
+
+
+def reference_bits(fmt, n, d):
+    """Correctly rounded bit pattern of n/d >= 0."""
+    if n == 0:
+        return 0
+    e = floor_log2(n, d)
+    s = fmt.denorm_min_exp if e < fmt.emin else e - fmt.p
+    m = round_ties_even(n, d << s) if s >= 0 else round_ties_even(n << -s, d)
+    if e < fmt.emin:
+        return m
+    if m == 2 << fmt.p:
+        m >>= 1
+        e += 1
+    if e > fmt.emax:
+        return fmt.infinite_power << fmt.p
+    return ((e - fmt.minimum_exponent) << fmt.p) | (m - (1 << fmt.p))
+
+
+def ceil_log2(x):
+    z = 0
+    while (1 << z) < x:
+        z += 1
+    return z
+
+
+def table_entry(q):
+    """T[q] as stored in fast_table.h (see table_generation.py)."""
+    if q < 0:
+        power5 = 5**-q
+        z = ceil_log2(power5)
+        b = z + 127 if q >= -27 else 2 * z + 128
+        c = 2**b // power5 + 1
+        while c >= 1 << 128:
+            c //= 2
+        return c
+    power5 = 5**q
+    while power5 < 1 << 127:
+        power5 *= 2
+    while power5 >= 1 << 128:
+        power5 //= 2
+    return power5
+
+
+TABLE = {q: table_entry(q) for q in range(-342, 309)}
+
+
+def compute_product_approximation(q, w, bit_precision):
+    t = TABLE[q]
+    first = w * (t >> 64)
+    high, low = first >> 64, first & MASK64
+    if (high & (MASK64 >> bit_precision)) == MASK64 >> bit_precision:
+        second_high = (w * (t & MASK64)) >> 64
+        low = (low + second_high) & MASK64
+        if second_high > low:
+            high += 1
+    return high, low
+
+
+def compute_float(fmt, q, w, subnormal_tie_test=True):
+    """Model of compute_float in decimal_to_binary.h. Returns
+    (mantissa, power2, truncated mantissa, k) with truncated == floor(w 10^q 2^k)."""
+    if w == 0 or q < fmt.smallest_power_of_ten:
+        return 0, 0, None, None
+    if q > fmt.largest_power_of_ten:
+        return 0, fmt.infinite_power, None, None
+    lz = 64 - w.bit_length()
+    w <<= lz
+    high, low = compute_product_approximation(q, w, fmt.p + 3)
+    upperbit = high >> 63
+    shift = upperbit + 64 - fmt.p - 3
+    mantissa = high >> shift
+    power2 = ((152170 + 65536) * q >> 16) + 63 + upperbit - lz - fmt.minimum_exponent
+    tie_q = fmt.min_exponent_round_to_even <= q <= fmt.max_exponent_round_to_even
+    if power2 <= 0:
+        if -power2 + 1 >= 64:
+            return 0, 0, None, None
+        extra = -power2 + 1
+        mantissa >>= extra
+        truncated, k = mantissa, fmt.bias
+        if (subnormal_tie_test and low <= 1 and tie_q and (mantissa & 3) == 1
+                and ((mantissa << extra) << shift) == high):
+            mantissa &= ~1
+        mantissa += mantissa & 1
+        mantissa >>= 1
+        return mantissa, int(mantissa >= 1 << fmt.p), truncated, k
+    truncated, k = mantissa, 1 + fmt.bias - power2
+    if low <= 1 and tie_q and (mantissa & 3) == 1 and (mantissa << shift) == high:
+        mantissa &= ~1
+    mantissa += mantissa & 1
+    mantissa >>= 1
+    if mantissa >= 2 << fmt.p:
+        mantissa = 1 << fmt.p
+        power2 += 1
+    mantissa &= ~(1 << fmt.p)
+    if power2 >= fmt.infinite_power:
+        return 0, fmt.infinite_power, truncated, k
+    return mantissa, power2, truncated, k
+
+
+def parsed_bits(fmt, w, q, subnormal_tie_test=True):
+    """from_chars on the digits w * 10^q. Clinger's fast path is one correctly
+    rounded operation on exact operands, so it equals the reference."""
+    if fmt.min_exponent_fast_path <= q <= fmt.max_exponent_fast_path and w <= 2 << fmt.p:
+        return reference_bits(fmt, *ratio10(w, q)), None, None
+    mantissa, power2, truncated, k = compute_float(fmt, q, w, subnormal_tie_test)
+    return mantissa | (power2 << fmt.p), truncated, k
+
+
+failures = 0
+
+
+def check(ok, text):
+    global failures
+    print(("  ok    " if ok else "  FAIL  ") + text)
+    failures += not ok
+
+
+def note(text):
+    print("  note  " + text)
+
+
+def continued_fraction(numer, denom):
+    cf = []
+    while denom:
+        quot, rem = divmod(numer, denom)
+        cf.append(quot)
+        numer, denom = denom, rem
+    return cf
+
+
+def convergents(cf):
+    p2, q2, p1, q1 = 0, 1, 1, 0
+    for a in cf:
+        p1, p2, q1, q2 = a * p1 + p2, p1, a * q1 + q2, q1
+        yield p1, q1
+
+
+def check_constants(fmt):
+    half_denorm = ratio(1, fmt.denorm_min_exp - 1)
+    q = fmt.smallest_power_of_ten
+    check(le(*ratio10(MAX_W - 1, q - 1), *half_denorm),
+          "smallest_power_of_ten=%d: q < %d rounds to zero" % (q, q))
+    check(not le(*ratio10(MAX_W - 1, q), *half_denorm),
+          "smallest_power_of_ten is tight")
+    threshold = ratio((4 << fmt.p) - 1, fmt.emax - fmt.p - 1)  # max + ulp/2
+    q = fmt.largest_power_of_ten
+    check(le(*threshold, *ratio10(1, q + 1)),
+          "largest_power_of_ten=%d: q > %d overflows" % (q, q))
+    check(not le(*threshold, *ratio10(1, q)), "largest_power_of_ten is tight")
+
+    two_p1 = 2 << fmt.p
+    q = fmt.max_exponent_fast_path
+    check(5**q <= two_p1 < 5 ** (q + 1),
+          "max_exponent_fast_path=%d: 10^q exact in T iff q <= %d" % (q, q))
+    check(fmt.powers_of_ten[: q + 1] == [10.0**i for i in range(q + 1)], "powers_of_ten")
+    check(fmt.max_mantissa[: q + 1] == [two_p1 // 5**i for i in range(q + 1)],
+          "max_mantissa[q] == floor(2^%d / 5^q)" % (fmt.p + 1))
+    overflow = not le(two_p1 * 10**q, 1, *ratio(two_p1 - 1, fmt.emax - fmt.p))
+    note("fast path product 2^%d * 10^%d %s the largest finite value"
+         % (fmt.p + 1, q, "exceeds" if overflow else "is below"))
+    k = -fmt.min_exponent_fast_path
+    check(0 <= k <= fmt.max_exponent_fast_path, "min_exponent_fast_path=%d" % -k)
+    if fmt.p <= 10:
+        # double rounding: w / 10^k rounded to `wide` bits, then to T
+        for wide in (24, 53, 64):
+            bad = 0
+            for kk in range(1, fmt.max_exponent_fast_path + 1):
+                for w in range(1, two_p1 + 1):
+                    e = floor_log2(w, 10**kk)
+                    s = e - wide + 1
+                    m = round_ties_even(w, 10**kk << s) if s >= 0 else round_ties_even(w << -s, 10**kk)
+                    if reference_bits(fmt, *ratio(m, s)) != reference_bits(fmt, w, 10**kk):
+                        bad += 1
+            check(bad == 0, "w / 10^k, k <= %d, via %d-bit intermediate never double-rounds"
+                  % (fmt.max_exponent_fast_path, wide))
+
+
+def check_ties(fmt):
+    """A tie is m * 2^e with m odd. It is w * 10^q, w < 10^19, iff q <= e and
+    w = m 5^-q 2^(e-q) (q <= 0) or 5^q | m (q > 0). Ties whose lower neighbour
+    is even (m == 1 mod 4) must fall in the round-to-even range."""
+    p = fmt.p
+    lo, hi = fmt.min_exponent_round_to_even, fmt.max_exponent_round_to_even
+    q = 0
+    while ((1 << (p + 1)) + 1) * 5 ** (q + 1) < MAX_W:
+        q += 1
+    normal_lo = -q
+    q = 0
+    while 5 ** (q + 1) < 1 << (p + 2):
+        q += 1
+    while True:
+        c = 1
+        while 5**q * c <= 1 << (p + 1):
+            c += 4
+        if 5**q * c < 1 << (p + 2) and q <= fmt.emax - p - 1:
+            break
+        q -= 1
+    normal_hi = q
+    e = fmt.denorm_min_exp - 1
+    sub = []
+    q = e
+    while 5**-q << (e - q) < MAX_W:
+        sub.append(q)
+        q -= 1
+    check(lo <= normal_lo and hi >= normal_hi,
+          "round_to_even range [%d, %d] covers normal ties [%d, %d]" % (lo, hi, normal_lo, normal_hi))
+    if sub:
+        check(lo <= min(sub), "and the subnormal ties at q in [%d, %d]" % (min(sub), max(sub)))
+    else:
+        note("no subnormal tie has fewer than 20 digits")
+    longest = len(str(((1 << (p + 2)) - 1) * 5**fmt.bias))
+    check(fmt.max_digits >= longest, "max_digits=%d >= %d, the longest tie" % (fmt.max_digits, longest))
+
+
+def check_table(fmt):
+    """Facts behind the tie test (product.low <= 1 and dropped bits zero)."""
+    q_lo, q_hi = max(fmt.smallest_power_of_ten, -27), -1
+    check(fmt.min_exponent_round_to_even >= -27, "tie test stays where 5^-q < 2^64")
+    check(all(TABLE[q] & MASK64 for q in range(q_lo, q_hi + 1)),
+          "T[q] low word nonzero for %d <= q <= %d" % (q_lo, q_hi))
+    fake = False
+    for q in range(q_lo, q_hi + 1):
+        t_hi = TABLE[q] >> 64
+        v2 = (t_hi & -t_hi).bit_length() - 1
+        z = ceil_log2(5**-q)
+        for shift in (61 - fmt.p, 62 - fmt.p):
+            # skipped second product: w * t_hi == 0 or 1 mod 2^(shift+64)
+            if shift + 64 - v2 <= 63:
+                fake = True
+            if v2 == 0 and (1 << 63) <= pow(t_hi, -1, 1 << (shift + 64)) < 1 << 64:
+                fake = True
+            # computed second product: a non-tie within 2^-62 of a rational
+            # with denominator 2^(z-1-shift) is impossible
+            if z - 1 - shift > 62:
+                fake = True
+    check(not fake, "no false tie for %d <= q <= %d" % (q_lo, q_hi))
+    bad = False
+    for q in range(fmt.smallest_power_of_ten, fmt.largest_power_of_ten + 1):
+        for _, w in convergents(continued_fraction(TABLE[q], 2**137)):
+            if w >= 2**64:
+                break
+            if (TABLE[q] * w) % 2**137 > 2**137 - 2**64:
+                bad = True
+    check(not bad, "Mushtak-Lemire condition for %d <= q <= %d"
+          % (fmt.smallest_power_of_ten, fmt.largest_power_of_ten))
+
+
+def sweep_q(args):
+    fmt, grid, q, width = args
+    cands = set()
+    for num, exp2, _ in grid:
+        n, d = ratio(num, exp2)
+        if q >= 0:
+            d *= 10**q
+        else:
+            n *= 10**-q
+        cands.update(w for w in range(n // d - width, n // d + width + 2) if 1 <= w < MAX_W)
+    wrong, wrong_without, bad_trunc = [], 0, 0
+    for w in cands:
+        n, d = ratio10(w, q)
+        want = reference_bits(fmt, n, d)
+        got, truncated, k = parsed_bits(fmt, w, q)
+        if got != want:
+            wrong.append((w, q, got, want))
+        if truncated is not None:
+            exact = (n << k) // d if k >= 0 else n // (d << -k)
+            bad_trunc += truncated != exact
+        if parsed_bits(fmt, w, q, False)[0] != want:
+            wrong_without += 1
+    return len(cands), wrong, wrong_without, bad_trunc
+
+
+def sweep(fmt, width, jobs):
+    """Compare the model with the reference on every w within `width` of
+    G / 10^q for every value, tie and the overflow threshold G."""
+    grid = fmt.grid()
+    qs = range(fmt.smallest_power_of_ten - 1, fmt.largest_power_of_ten + 2)
+    tasks = [(fmt, grid, q, width) for q in qs]
+    with multiprocessing.Pool(jobs) as pool:
+        results = pool.map(sweep_q, tasks)
+    wrong = [x for r in results for x in r[1]]
+    note("swept %d pairs (w, q)" % sum(r[0] for r in results))
+    check(sum(r[3] for r in results) == 0, "truncated mantissa is exact")
+    check(not wrong, "compute_float matches the reference")
+    for w, q, got, want in sorted(wrong)[:20]:
+        print("        w=%d q=%d got 0x%04x want 0x%04x" % (w, q, got, want))
+    note("without the subnormal tie test: %d wrong" % sum(r[2] for r in results))
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--types", default="double,float,float16,bfloat16")
+    ap.add_argument("--sweep-width", type=int, default=3)
+    ap.add_argument("--jobs", type=int, default=os.cpu_count() or 1)
+    args = ap.parse_args()
+    for name in args.types.split(","):
+        fmt = Format(name.strip())
+        print("== %s" % fmt.name)
+        check_constants(fmt)
+        check_table(fmt)
+        check_ties(fmt)
+        if fmt.p <= 10 and args.sweep_width > 0:
+            sweep(fmt, args.sweep_width, args.jobs)
+    print("%d failure(s)" % failures)
+    return 1 if failures else 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

tests/CMakeLists.txt

@@ -94,6 +94,8 @@ endif()
 if (FASTFLOAT_FIXEDWIDTH_TESTS)
   fast_float_add_cpp_test(fixedwidthtest)
   target_compile_features(fixedwidthtest PUBLIC cxx_std_23)
+  fast_float_add_cpp_test(exhaustive16_midpoint)
+  target_compile_features(exhaustive16_midpoint PUBLIC cxx_std_23)
 endif()
 
 option(FASTFLOAT_EXHAUSTIVE "Exhaustive tests" OFF)

tests/basictest.cpp

@@ -1832,6 +1832,20 @@ TEST_CASE("float16.inf") {
          std::errc::result_out_of_range);
   verify("3.5028234666e38", std::numeric_limits<std::float16_t>::infinity(),
          std::errc::result_out_of_range);
+  // overflow in the fast path
+  verify("656e2", std::numeric_limits<std::float16_t>::infinity(),
+         std::errc::result_out_of_range);
+  verify("-656e2", -std::numeric_limits<std::float16_t>::infinity(),
+         std::errc::result_out_of_range);
+  verify("7e4", std::numeric_limits<std::float16_t>::infinity(),
+         std::errc::result_out_of_range);
+  verify("2048e4", std::numeric_limits<std::float16_t>::infinity(),
+         std::errc::result_out_of_range);
+  verify("655e2", 0x1.ffcp+15f16); // 65500 rounds to max
+  // max + ulp/2 rounds to even
+  verify("6552e1", std::numeric_limits<std::float16_t>::infinity(),
+         std::errc::result_out_of_range);
+  verify("65519", 0x1.ffcp+15f16);
 }
 
 TEST_CASE("float16.general") {
@@ -1852,6 +1866,22 @@ TEST_CASE("float16.general") {
   verify("-0.000000059604644775390625", -0x1p-24f16);
   verify("-5.9604644775390625e-8", -0x1p-24f16);
 
+  // subnormal ties round to even
+  verify("2.98023223876953125e-8", 0.0f16, std::errc::result_out_of_range);
+  verify("0.0000000298023223876953125", 0.0f16, std::errc::result_out_of_range);
+  verify("-2.98023223876953125e-8", -0.0f16, std::errc::result_out_of_range);
+  verify("8.94069671630859375e-8", 0x1p-23f16);
+  verify("1.490116119384765625e-7", 0x1p-23f16);
+  verify("9.834766387939453125e-7", 0x1p-20f16);
+  // just above or below a tie
+  verify("2.980232238769531251e-8", 0x1p-24f16);
+  verify("2.980232238769531249e-8", 0.0f16, std::errc::result_out_of_range);
+  verify("1.4901161193847656251e-7", 0x1.8p-23f16);
+  verify("1.4901161193847656249e-7", 0x1p-23f16);
+  // more than 19 digits
+  verify("2.9802322387695312500e-8", 0.0f16, std::errc::result_out_of_range);
+  verify("1.4901161193847656250000e-7", 0x1p-23f16);
+
   verify("-1e-999", -0.0f16, std::errc::result_out_of_range);
   verify("6.0975551605224609375", 0x1.864p+2f16);
   verify_runtime(append_zeros("6.0975551605224609375", 655), 0x1.864p+2f16);

tests/exhaustive16_midpoint.cpp

@@ -0,0 +1,163 @@
+// For std::float16_t and std::bfloat16_t, parse the exact decimal expansion
+// of every finite value and of every midpoint between adjacent values, plus
+// strings just above and below each midpoint and midpoints padded past 19
+// digits. The expected bits follow from the bit pattern that generated the
+// string, so no C library reference is needed.
+#include "fast_float/fast_float.h"
+
+#include <cstdint>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <stdfloat>
+#include <string>
+#include <system_error>
+
+namespace {
+
+// Decimal digit strings, most significant digit first, no leading zeros.
+std::string times(std::string const &digits, unsigned factor) {
+  std::string out;
+  unsigned carry = 0;
+  for (size_t i = digits.size(); i-- > 0;) {
+    unsigned v = static_cast<unsigned>(digits[i] - '0') * factor + carry;
+    out.push_back(static_cast<char>('0' + v % 10));
+    carry = v / 10;
+  }
+  while (carry != 0) {
+    out.push_back(static_cast<char>('0' + carry % 10));
+    carry /= 10;
+  }
+  return std::string(out.rbegin(), out.rend());
+}
+
+std::string minus_one(std::string digits) {
+  size_t i = digits.size();
+  while (i-- > 0) {
+    if (digits[i] != '0') {
+      digits[i] = static_cast<char>(digits[i] - 1);
+      break;
+    }
+    digits[i] = '9';
+  }
+  return digits;
+}
+
+// num * 2^exp2 written as "<digits>e<exponent>" with an exact expansion.
+struct Decimal {
+  std::string digits;
+  int exponent;
+
+  std::string str() const { return digits + "e" + std::to_string(exponent); }
+
+  Decimal above() const { return {digits + "1", exponent - 1}; }
+
+  Decimal below() const { return {minus_one(digits) + "9", exponent - 1}; }
+
+  Decimal padded(size_t count) const {
+    return {digits + std::string(count, '0'),
+            exponent - static_cast<int>(count)};
+  }
+};
+
+Decimal exact_decimal(uint64_t num, int exp2) {
+  std::string digits = std::to_string(num);
+  if (exp2 >= 0) {
+    for (int i = 0; i < exp2; i++) {
+      digits = times(digits, 2);
+    }
+    return {digits, 0};
+  }
+  for (int i = 0; i < -exp2; i++) {
+    digits = times(digits, 5);
+  }
+  return {digits, exp2};
+}
+
+template <typename T> struct Layout;
+
+template <> struct Layout<std::float16_t> {
+  static constexpr char const *name = "float16";
+  static constexpr int p = 10;
+  static constexpr int bias = 25; // value = mantissa * 2^(field - bias)
+  static constexpr uint16_t infinity = 0x7C00;
+};
+
+template <> struct Layout<std::bfloat16_t> {
+  static constexpr char const *name = "bfloat16";
+  static constexpr int p = 7;
+  static constexpr int bias = 134;
+  static constexpr uint16_t infinity = 0x7F80;
+};
+
+struct Failure {
+  long count = 0;
+  long checked = 0;
+};
+
+template <typename T>
+void expect(Failure &f, std::string const &s, uint16_t want) {
+  T value{};
+  auto r = fast_float::from_chars(s.data(), s.data() + s.size(), value);
+  uint16_t got = 0;
+  std::memcpy(&got, &value, sizeof(got));
+  bool zero_or_inf = want == 0 || want == Layout<T>::infinity;
+  std::errc want_ec =
+      zero_or_inf ? std::errc::result_out_of_range : std::errc();
+  f.checked++;
+  if (got != want || r.ec != want_ec || r.ptr != s.data() + s.size()) {
+    f.count++;
+    if (f.count <= 20) {
+      std::printf(
+          "%s: \"%s\" parsed to 0x%04x (ec %d), expected 0x%04x (ec %d)\n",
+          Layout<T>::name, s.c_str(), got, static_cast<int>(r.ec), want,
+          static_cast<int>(want_ec));
+    }
+  }
+}
+
+template <typename T> Failure sweep() {
+  Failure f;
+  constexpr int p = Layout<T>::p;
+  constexpr uint16_t mantissa_mask = static_cast<uint16_t>((1u << p) - 1);
+  for (uint32_t b = 0; b < Layout<T>::infinity; b++) {
+    uint16_t bits = static_cast<uint16_t>(b);
+    int field = bits >> p;
+    uint64_t num = bits & mantissa_mask;
+    int exp2 = (field == 0 ? 1 : field) - Layout<T>::bias;
+    if (field != 0) {
+      num |= uint64_t(1) << p;
+    }
+    if (num != 0) {
+      Decimal d = exact_decimal(num, exp2);
+      expect<T>(f, d.str(), bits);
+      expect<T>(f, d.padded(21).str(), bits);
+    }
+    // The tie between bits and bits + 1 is (2 num + 1) * 2^(exp2 - 1); the
+    // last one (bits + 1 == infinity) is the overflow threshold.
+    Decimal tie = exact_decimal(2 * num + 1, exp2 - 1);
+    uint16_t lower = bits;
+    uint16_t upper = static_cast<uint16_t>(bits + 1);
+    uint16_t even = (bits & 1) ? upper : lower;
+    expect<T>(f, tie.str(), even);
+    expect<T>(f, tie.padded(21).str(), even);
+    expect<T>(f, tie.above().str(), upper);
+    expect<T>(f, tie.below().str(), lower);
+  }
+  return f;
+}
+
+} // namespace
+
+int main() {
+  Failure f16 = sweep<std::float16_t>();
+  std::printf("float16: %ld strings, %ld failures\n", f16.checked, f16.count);
+  Failure bf16 = sweep<std::bfloat16_t>();
+  std::printf("bfloat16: %ld strings, %ld failures\n", bf16.checked,
+              bf16.count);
+  if (f16.count != 0 || bf16.count != 0) {
+    return EXIT_FAILURE;
+  }
+  std::printf("all ok\n");
+  return EXIT_SUCCESS;
+}