coffee/fast_float
1
0
Fork 0
mirror of https://github.com/fastfloat/fast_float.git synced 2026-10-28 18:07:16 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #408 from jadidbourbaki/float16-exactness

Browse Source 
Fix std::float16_t rounding of subnormal ties and fast-path overflow
This commit is contained in:
Daniel Lemire 2026-09-22 14:11:44 -04:00 committed by GitHub
commit cdc385d058
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
7 changed files with 667 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
include/fast_float/decimal_to_binary.h
View File

@ -156,10 +156,21 @@ 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 (binary::subnormal_ties_possible() && (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.

43
include/fast_float/float_common.h
View File

@ -588,6 +588,8 @@ 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 bool subnormal_ties_possible();
static constexpr int largest_power_of_ten();
static constexpr int smallest_power_of_ten();
static constexpr T exact_power_of_ten(int64_t power);
@ -619,6 +621,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 /
@ -634,9 +637,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
@ -853,7 +854,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 <>
@ -865,7 +868,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 <>
@ -889,7 +894,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 <>
@ -976,7 +982,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 <>
@ -1012,7 +1019,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 <>
@ -1021,6 +1029,25 @@ 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)());
}
// A subnormal needs w * 10^q < 2^(minimum_exponent() + 1), so q is at most
// (minimum_exponent() + 1) * log10(2), with 1233/4096 < log10(2). Only
// std::float16_t has such q in its round-to-even range. We compare
// 4096 * q with (minimum_exponent() + 1) * 1233 to avoid right-shifting a
// negative value (implementation-defined before C++20).
template <typename T>
inline constexpr bool binary_format<T>::subnormal_ties_possible() {
return min_exponent_round_to_even() * 4096 <= (minimum_exponent() + 1) * 1233;
}
template <>
inline constexpr uint64_t
binary_format<double>::max_mantissa_fast_path(int64_t power) {

6
include/fast_float/parse_number.h
View File

@ -241,6 +241,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;

416
script/format_parameters.py Normal file
View File

@ -0,0 +1,416 @@
#
# 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")
# subnormal_ties_possible() in float_common.h: an upper bound on the
# largest q that can give a subnormal, compared with the tie range
q_sub = ((fmt.minimum_exponent + 1) * 1233) >> 12
check(10**q_sub < 2**fmt.emin < 10 ** (q_sub + 2), "q <= %d bounds the subnormal exponents" % q_sub)
check((lo <= q_sub) == bool(sub), "subnormal_ties_possible() == %s" % bool(sub))
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())

2
tests/CMakeLists.txt
View File

@ -95,6 +95,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)

30
tests/basictest.cpp
View File

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

163
tests/exhaustive16_midpoint.cpp Normal file
View File

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