fast_float_strtod: fix ±1 ULP rounding mismatch in widened fast path#15111
Conversation
Imports the precomputed extended-precision powers of five from upstream fast_float (fast_table.h, MIT-licensed). Covers the full range needed to handle any 64-bit decimal mantissa producing a finite double: 5^-342 ... 5^308, 651 entries, 1302 uint64_t values total. This is data-only; the next commit wires it into compute_float() to replace the current __SIZEOF_INT128__ widened branch which has known 1-ULP rounding mismatches vs strtod() on 17-19 digit mantissas (e.g. 9007199255094284e-19 and 2489830482329185244e1). Source: https://github.com/fastfloat/fast_float include/fast_float/fast_table.h
Adds the algorithmic body of the Eisel-Lemire decimal-to-binary
conversion as a double-only specialisation. Components:
- value128 / adjusted_mantissa structs (port of upstream's typedefs)
- leading_zeroes_u64 (__builtin_clzll wrapper)
- full_multiplication (64x64 -> 128, __uint128_t fast path with a
32-bit-target fallback so this remains correct on platforms
without 128-bit ints)
- eisel_lemire_power, compute_product_approximation_d
- compute_float_d — the main algorithm, returns (mantissa, power2)
- am_to_double — packs the result into IEEE-754 binary64 layout
Mathematical guarantees: for inputs with ≤19 mantissa digits and the
exponent ranges produced by parse_number_string, the result is provably
bit-exact with strtod (Mushtak & Lemire, 'Fast Number Parsing Without
Fallback'). compute_float_d's negative-power2 branch is therefore
unreachable from our parser, but is preserved in the API so the next
commit can wire a defensive strtod fallback for any future caller.
Still dead code at this stage; the next commit replaces the buggy
__SIZEOF_INT128__ widened branch in parse_number_string with a call
into compute_float_d + am_to_double.
Source: https://github.com/fastfloat/fast_float
include/fast_float/{decimal_to_binary,float_common}.h
The previous __SIZEOF_INT128__ widened branch used a homebrew '(double)hi * 2^64 + (double)lo' conversion that double-rounded the 128-bit product, producing ±1 ULP mismatches vs strtod() on inputs where the high half exceeded 2^53 or where the integer division truncated the remainder. Two reported reproducers: 9007199255094284e-19 -> -1 ULP 2489830482329185244e1 -> +1 ULP Replaces that branch with a call into the Eisel-Lemire core (compute_float_d + am_to_double) added in the previous commit. For inputs the Clinger fast path can handle (mantissa <= 2^53, |exp| <= 22) we keep Clinger — it's a single double multiply and stays the cheapest path. Everything else now goes through Eisel-Lemire, which is correctly-rounded by construction (Mushtak & Lemire, 'Fast Number Parsing Without Fallback') and covers the full double exponent range, not just |exp| <= 19. Existing fastfloat unit tests pass (98/98). Both reproducer inputs above now match strtod() bit-exactly.
…cases
Adds 22 new ff_testcase entries (subnormal boundaries, round-half-to-
even tie cases, near-DBL_MAX overflow, wide-exponent inputs now reachable
via Eisel-Lemire) plus an explicit differential block that bit-compares
fast_float_strtod() against libc strtod() on 22 hard inputs covering:
- 2^53 boundary mantissas (9007199254740992..996)
- DBL_MAX / smallest-normal / smallest-subnormal limits
- The two reproducer inputs the previous widened branch missed
(9007199255094284e-19 and 2489830482329185244e1)
- 19-digit Mushtak-Lemire stress mantissas
- Common scientific constants
Test count: 98 -> 140. All pass.
Reported-by: vitahlin
Suggested-by: vitahlin (rounding-boundary cases on PR redis#14661)
Whitespace-only change. The 651-entry table previously occupied 651 lines (one entry per line, matching upstream fast_float's per-entry formatting). For a flat C uint64_t array this is wasteful; pack 8 entries (16 hex literals, ~280 chars) per line instead. Cuts the table body from 651 to 82 lines. No behaviour change. All 140 fastfloat tests still pass.
🤖 Augment PR SummarySummary: This PR fixes a ±1 ULP mismatch between Redis’ Changes:
Tests: Expands the REDIS_TEST suite with the reported reproducers, subnormal/tie-boundary cases, and a differential bitwise comparison loop against libc 🤖 Was this summary useful? React with 👍 or 👎 |
| value128 firstproduct = full_multiplication(w, power_of_five_128[index]); | ||
| /* For double, bit_precision = mantissa_explicit_bits (52) + 3 = 55. */ | ||
| const uint64_t precision_mask = | ||
| (uint64_t)0xFFFFFFFFFFFFFFFFULL >> 55; |
There was a problem hiding this comment.
src/fast_float_strtod.c:926 — precision_mask looks computed with the wrong shift: UINT64_MAX >> 55 leaves only 9 one-bits, but the comment indicates bit_precision is 55 so this mask likely needs 55 one-bits. If this mask is wrong, the “rounding-boundary fixup” path in compute_product_approximation_d() may rarely trigger (or trigger incorrectly), reintroducing hard-to-reproduce rounding mismatches vs strtod().
Severity: high
🤖 Was this useful? React with 👍 or 👎, or 🚀 if it prevented an incident/outage.
There was a problem hiding this comment.
When these 9 bits are all 1 (i.e., firstproduct.high & 0x1FF == 0x1FF), it indicates that the value is very close to the round-half-to-even boundary, so the extension table (power_of_five_128[index+1]) needs to be triggered for a second-pass correction.
so this is correct.
| }; | ||
| run_ff_tests(decimal_ok, COUNTOF(decimal_ok), 0); | ||
|
|
||
| /* Differential cross-check: every accepted input must produce the |
There was a problem hiding this comment.
src/fast_float_strtod.c:1410 — The differential block compares fast_float_strtod() vs libc strtod(), but it would still pass if these inputs silently fall back to libc (e.g., if fast_float_try_fast() rejects them). Since this PR is specifically about correctness/perf of the widened fast path, it may be worth asserting that the key reproducers (and a couple of representative hard cases) are actually accepted by fast_float_try_fast() so a future regression doesn’t get masked.
Severity: low
🤖 Was this useful? React with 👍 or 👎, or 🚀 if it prevented an incident/outage.
Whitespace-only change. The previous 8-entries-per-line layout produced ~280-character lines, which hurt diff/code-review readability. 4 entries (8 hex literals) per line keeps lines around 140 characters. No behaviour change. All 140 fastfloat tests still pass.
|
I just realized that libc also use Eisel-Lemire algorithm. In fact, the main expense at present lies in the need to copy the memory. But I'm thinking that if we know for sure that the input is terminated with '\0', then we can directly fall back to using strtod(). [EDIT] ignore the above comment, tested, and it's not good enough. |
Reformatting only, no behavior change: - compute_float_d: collapse power2 normal-path computation onto a single line. - compute_float_d: collapse the subnormal→normal post-rounding power2 ternary onto a single line. - parse_number_string: keep the EISEL_LEMIRE bounds check condition on one line and place `return 0;` on its own indented line. vitahlin's two reproducers still match libc strtod() bit-exactly (0 ULP delta), and unit/scripting passes.
|
Hi @sundb — addressed your review in
Reformatting only — no behavior change. @vitahlin's two reproducers still match libc |
Summary
Fixes a ±1 ULP rounding mismatch between
fast_float_strtod()and libcstrtod()in the widened (mantissa > 2^53) fast path introduced by #15061. Reported by @vitahlin in #14661 (comment) with two minimal reproducers:Redis treats
strtod()as the fallback forfast_float_strtod(), so every fast-path-accepted input is contractually expected to be bit-exact withstrtod(). The two cases above are accepted by the widened branch but produce a different IEEE-754 representation, breaking the contract.Root cause
The widened branch added in #15061 used a homebrew shortcut to convert a 128-bit integer product to a double:
This is not a single-rounding operation —
(double)hirounds whenhi > 2^53, and the subsequent+ (double)lorounds again. For inputs near the round-half-to-even boundary, the two roundings can compose into the wrong direction.The negative-exponent branch is doubly affected: integer division
scaled / divisortruncates the remainder before the conversion, so even a hypothetical correcthi*2^64+lostep would round down on inputs that should round up.Fix
Replace the homebrew widened branch with the Eisel-Lemire algorithm from upstream
fast_float. This is the same algorithm thatfast_float's own widened path uses; bit-exact-with-strtod for inputs with ≤19-digit mantissa is proved in:Pieces ported (MIT-licensed, from
fast_float/include/fast_float/fast_table.h,decimal_to_binary.h,float_common.h):5^-342 ... 5^308, 651 entries) — pure data.compute_product_approximation— 128-bit multiply with high-half rounding-boundary fixup.compute_float— the main algorithm; returns(mantissa, power2)ready to be packed.am_to_double— IEEE-754 binary64 bit-pack.__builtin_clzlland__uint128_twrappers (with a 32-bit fallback for portability).Clinger's strict fast path (
mantissa ≤ 2^53and|exp| ≤ 22) is kept unchanged — it's a single double multiply/divide and is faster than Eisel-Lemire on its domain. Only the buggy widened branch is replaced.The port stays minimal:
float, nolong double)digit_compare unreachable fromparse_number_string's ≤19-digit mantissa per the Mushtak-Lemire proof, but a defensiveam.power2 < 0check is preserved that falls back to libcstrtod()if any future caller widens the input domain.Why not just revert #15061?
Considered. Reverting restores correctness at the cost of the +73-84 % zset listpack-load wins #15061 measured on 17-19 digit double scores. Eisel-Lemire is the algorithm that gives both correctness and the wider mantissa range — preserving #15061's wins while fixing the rounding regression.
A "tightened admission filter" (only accept widened-path inputs where the conversion happens to be single-rounding) was also considered. The math shows the filter conditions are essentially unsatisfied for typical inputs (
lo == 0requires the 128-bit product be divisible by 2^64; only ~1 in 10^13 random inputs qualify), making it equivalent to a revert with extra dead code. Eisel-Lemire is the only widened-path solution that preserves perf on the typical case.Testing
./src/redis-server test fastfloat— 140 assertions, all pass.decimal_okcasesdecimal_okcases5e-324, smallest/largest subnormals), round-half-to-even ties, near-DBL_MAXoverflow, repunit mantissas, wide-exponent cases reachable via Eisel-Lemire9007199255094284e-19and2489830482329185244e1) — both now matchstrtod()exactlymemcmpbetweenfast_float_strtod()andstrtod()on hard cases (2^53 boundaries,DBL_MAX, smallest-subnormal, Mushtak-Lemire stress range, scientific constants)Performance — preliminary, no regression detected
ZSET-heavy benchmarks on
redis/redis:unstablebaseline vs this branch,oss-standalonetopology, AWS bare-metal runners (Intel SPRm7i.metal-24xland ARM Graviton4m8g.metal-24xl), n=3 dps per cell on the PR side:1Mkeys-load-zset-listpack-with-100-elements-double-score(#15061's headline)1key-zset-1K-elements-zrange-all-elements1key-zset-600K-elements-zrangestore-1K-elements1key-zset-1M-elements-zremrangebyscore-pipeline-10No regression beyond fleet baseline noise (~1.5 %) on any test on either arch. PR #15061's +73-84 % win on
load-zset-listpack-with-100-elements-double-scoreis preserved (Intel flat at −2.45 %, well within fleet drift; ARM improves +7.59 %, likely from Eisel-Lemire's wider exponent coverage hitting the fast path on more input shapes than the previous branch). Final cell + same-day baseline pinning will be added in a follow-up comment once the runs complete.Commits
The branch is structured as five atomic commits to ease review:
5e16ef626— Add the 128-bit power-of-five table (data only, no behaviour change)56b981261— Portcompute_float+ helpers (dead code at this point, verifies build)31fb41cbd— Wire Eisel-Lemire intoparse_number_string, remove the buggy widened branchadabcc0ec— ExtendREDIS_TESTblock with reproducers + boundary cases + differential vsstrtod()53c0aaaa8— Pack the table at 8 entries per line (whitespace-only)Each commit builds cleanly and the existing test suite passes after each one.
Credits
Test plan
make testpassesmake REDIS_CFLAGS=-DREDIS_TEST -j && ./src/redis-server test fastfloat— 140/140 passstrtod()unstablebaseline (full table in PR body; same-day baseline pinning to follow)Note
Medium Risk
Changes core decimal-to-double conversion logic and adds a large precomputed table; mistakes could cause subtle numeric rounding differences or regress parsing behavior/perf across platforms.
Overview
Fixes
fast_float_strtod()correctness for cases where mantissas exceed2^53by replacing the prior 128-bit “widened” conversion shortcut with a port of the Eisel–Lemire decimal-to-binary algorithm, including IEEE-754 packing (am_to_double) and required helpers.Adds the upstream
5^q128-bit lookup table (coveringqin[-342, 308]) and updatesparse_number_string()to choose between the existing Clinger fast path and the new Eisel–Lemire path, falling back to libcstrtod()when out of range.Expands the
REDIS_TESTsuite with reproducers for the reported ±1 ULP issues plus additional boundary and differential tests to ensure bit-for-bit parity withstrtod()across subnormals, ties-to-even, near-overflow, and wide exponents.Reviewed by Cursor Bugbot for commit aeec6af. Bugbot is set up for automated code reviews on this repo. Configure here.