diff --git a/docs/roadmap.md b/docs/roadmap.md index 829fb38..a8c729c 100644 --- a/docs/roadmap.md +++ b/docs/roadmap.md @@ -89,13 +89,18 @@ For completeness — these came up alongside the above but are *not* blocked by `bound` over a set of grids) rather than collapsing to the hull. Large surface (every operator/predicate would need a union story), so this stays a design sketch until a concrete use case demands it. -- **128-bit rounded store** — the cold assignment path forms `(rhs − Lower)/Notch` as - an exact 64-bit rational before rounding; a full-mantissa `double`-derived source - (denominator 2^52+) on a grid with large `|Lower|` can need 65+ bits *before* the - round even though the rounded slot index is tiny. Mixed-sign offsets are rescued in - 128-bit and the rest report `errc::overflow` today; computing the *rounded* index - directly in 128-bit (as `to_fixed` already does in `cmath.hpp`) would make those - stores exact instead of an error. +- **128-bit rounded store** — **done.** The cold assignment path forms + `(rhs − Lower)/Notch` as an exact 64-bit rational before rounding; a full-mantissa + `double`-derived source (denominator 2^52+) on a grid with large `|Lower|` can need + 65+ bits *before* the round even though the rounded slot index is tiny. When that + exact formation overflows, `wide_offset_quotient` (assignment.hpp) now computes the + slot directly — compile-time divisor factors gcd-reduced first, then one 128×64 + multiply and one 128÷64 divide (`mul128`/`divmod128` in rational.hpp, native + `__int128` or portable shift-subtract). Rounding in that regime is the offset rule + (`round_offset`), the same fallback `round_quotient` uses past 64 bits. Only + results beyond even the 128-bit envelope report `errc::overflow`. One caveat: at + *constant evaluation* the transient rational overflow still surfaces as the + intentional `constexpr_error` diagnostic before the fallback can engage. - **Modules / compile-time-footprint work** — parked for later; modules is C++20, not a blocker. diff --git a/include/bound/detail/assignment.hpp b/include/bound/detail/assignment.hpp index e66cd3b..d4235b8 100644 --- a/include/bound/detail/assignment.hpp +++ b/include/bound/detail/assignment.hpp @@ -258,6 +258,92 @@ namespace bnd::detail action.fn(lhs, q); } + // 128-bit rounded store — the offset slot of an in-range rhs computed + // directly in wide arithmetic when the exact 64-bit rational formation + // of (rhs − Lower)/Notch overflows (full-mantissa fp-derived sources on + // grids with large |Lower|): + // slot + remainder/divisor = (rhs − Lower)·d_n / (a_dr·a_dl·n_n) + // The compile-time divisor factors are gcd-reduced first, so the only + // wide operations are one 128×64 multiply and one 128÷64 divide. + // ok == false when the reduced divisor or dividend exceeds the 128-bit + // envelope (or rhs is out of range — callers check range first). + struct wide_quotient { umax slot; umax remainder; umax divisor; bool ok; }; + + static constexpr wide_quotient wide_offset_quotient(rational const& rv) + { + if constexpr (Notch == 0) + return {}; // continuous grids never index slots + else + { + constexpr umax n_l = Lower.Numerator; + constexpr umax a_dl = abs_den(Lower.Denominator); + constexpr bool low_neg = Lower.Denominator < 0; + constexpr umax n_n = Notch.Numerator; // Notch > 0: d_n > 0 + constexpr umax d_n = static_cast(Notch.Denominator); + + // Compile-time divisor part; a grid whose a_dl·n_n cannot fit umax + // is beyond the wide envelope entirely. + constexpr umax den_ct = []{ + umax p; + return mul_overflow(a_dl, n_n, &p) ? umax{0} : p; + }(); + if constexpr (den_ct == 0) + return {}; + else + { + constexpr umax g1 = std::gcd(d_n, den_ct); + constexpr umax d_n1 = d_n / g1; + constexpr umax den1 = den_ct / g1; + + const umax a_dr = abs_den(rv.Denominator); + const bool rhs_neg = rv.Denominator < 0; + + // Offset numerator over the common denominator a_dr·a_dl: + // s_r·n_r·a_dl − s_l·n_l·a_dr (≥ 0 for in-range rhs). + // Each product is < 2^127 (numerator < 2^64, denominator ≤ imax), + // so the same-sign sum below cannot carry out of 128 bits. + const u128 val = umul(rv.Numerator, a_dl); + const u128 low = umul(n_l, a_dr); + u128 offset; + if (!rhs_neg && low_neg) + offset = u128{val.hi + low.hi + (val.lo + low.lo < val.lo ? 1u : 0u), + val.lo + low.lo}; + else if (!rhs_neg && !low_neg) + { + if (cmp128(val, low) < 0) return {}; // rhs < Lower + offset = u128{val.hi - low.hi - (val.lo < low.lo ? 1u : 0u), + val.lo - low.lo}; + } + else if (rhs_neg && low_neg) + { + if (cmp128(low, val) < 0) return {}; // rhs < Lower + offset = u128{low.hi - val.hi - (low.lo < val.lo ? 1u : 0u), + low.lo - val.lo}; + } + else + return {}; // rhs < 0 ≤ Lower + + const umax g2 = std::gcd(d_n1, a_dr); + const umax d_n2 = d_n1 / g2; + const umax a_dr1 = a_dr / g2; + + umax divisor; + if (mul_overflow(a_dr1, den1, &divisor) + || divisor > static_cast(std::numeric_limits::max())) + return {}; + + const mul128_result dividend = mul128(offset, d_n2); + if (dividend.overflowed) + return {}; + + const divmod128_result qr = divmod128(dividend.value, divisor); + if (qr.quotient.hi != 0) + return {}; // slot beyond any 64-bit index space + return {qr.quotient.lo, qr.remainder, divisor, true}; + } + } + } + template static constexpr bool store_checked(L& lhs, R rhs, P&& policy, A&& action = {}) { @@ -346,16 +432,36 @@ namespace bnd::detail } // The exact quotient can overflow the 64-bit rational range (huge - // source denominator × fine notch). Report it as an overflow through - // the policy channel instead of dereferencing nullopt — a throw here - // would escape noexcept callers (the math engines) as terminate. + // source denominator × fine notch). Recompute the slot directly in + // 128-bit (wide_offset_quotient above); only a result beyond even + // that envelope reports errc::overflow — never a nullopt deref, + // which would escape noexcept callers (the math engines) as + // terminate. Rounding here is the offset rule (round_offset), the + // same semantics round_quotient falls back to past 64 bits. const auto quotient = (rhs - Lower)/Notch; if (!quotient.has_value()) [[unlikely]] { - if constexpr (error_action>) - { action.fn(lhs, errc::overflow, errc_message(errc::overflow)); return false; } - policy.report(errc::overflow); - return false; + const wide_quotient wide = wide_offset_quotient(rational{rhs}); + if (!wide.ok) + { + if constexpr (error_action>) + { action.fn(lhs, errc::overflow, errc_message(errc::overflow)); return false; } + policy.report(errc::overflow); + return false; + } + if (wide.remainder == 0) + { store_slot(wide.slot); return true; } + if constexpr (has_round_flag) + { store_slot(round_offset(wide.slot, wide.remainder, wide.divisor)); return true; } + if (policy.round_check()) [[unlikely]] + { + if constexpr (error_action>) + { action.fn(lhs, errc::rounding_error, errc_message(errc::rounding_error)); return false; } + policy.report(errc::rounding_error); + return false; + } + store_slot(round_offset(wide.slot, wide.remainder, wide.divisor)); + return true; } rational raw = *quotient; umax den = static_cast(raw.Denominator); diff --git a/include/bound/detail/rational.hpp b/include/bound/detail/rational.hpp index c2cf9c0..95d32a0 100644 --- a/include/bound/detail/rational.hpp +++ b/include/bound/detail/rational.hpp @@ -50,6 +50,44 @@ namespace bnd::detail constexpr std::strong_ordering cmp128(u128 a, u128 b) { return (a.hi != b.hi) ? (a.hi <=> b.hi) : (a.lo <=> b.lo); } + // 128×64 product with an overflow flag (result beyond 128 bits). + struct mul128_result { u128 value; bool overflowed; }; + constexpr mul128_result mul128(u128 a, umax b) + { + const u128 low = umul(a.lo, b); + const u128 high = umul(a.hi, b); + const umax hi_sum = high.lo + low.hi; + return {u128{hi_sum, low.lo}, high.hi != 0 || hi_sum < low.hi}; + } + + // Quotient/remainder of a 128-bit dividend by a 64-bit divisor. Requires + // 1 <= d <= imax_max (the rational-denominator domain) so the portable + // partial remainder can never overflow when shifted. + struct divmod128_result { u128 quotient; umax remainder; }; + constexpr divmod128_result divmod128(u128 n, umax d) + { +#if defined(__SIZEOF_INT128__) + using u128n = unsigned __int128; + const u128n wide = (static_cast(n.hi) << 64) | n.lo; + const u128n q = wide / d; + return {u128{static_cast(q >> 64), static_cast(q)}, + static_cast(wide % d)}; +#else + // Portable (MSVC): restoring shift-subtract divide — the same construction + // as cmath.hpp's to_fixed fallback. + u128 q{0, 0}; + umax r = 0; + for (int i = 127; i >= 0; --i) + { + r = (r << 1) | ((i >= 64 ? (n.hi >> (i - 64)) : (n.lo >> i)) & 1u); + q.hi = (q.hi << 1) | (q.lo >> 63); + q.lo <<= 1; + if (r >= d) { r -= d; q.lo |= 1; } + } + return {q, r}; +#endif + } + //--------------------------------------------------------------------------- // trim //--------------------------------------------------------------------------- diff --git a/include/bound/generic.hpp b/include/bound/generic.hpp index 5e0aca2..9582b6d 100644 --- a/include/bound/generic.hpp +++ b/include/bound/generic.hpp @@ -398,6 +398,25 @@ namespace bnd template inline constexpr bool HasPolicy = has_flag(BoundPolicy, F) || plain

::test(F); + // Rounds the split offset quotient q + r/den (r < den ≤ imax_max) per L's + // rounding policy — q/r form so no expression can overflow umax + // (num + den/2 could, for num near umax). Shared by round_quotient's + // offset rule and the 128-bit wide store (assignment.hpp). + template + [[nodiscard]] constexpr umax round_offset(umax q, umax r, umax den) noexcept + { + if constexpr (HasPolicy) return (r * 2 >= den) ? q + 1 : q; + else if constexpr (HasPolicy) return q; + else if constexpr (HasPolicy) return (r != 0) ? q + 1 : q; + else if constexpr (HasPolicy) + { + if (r * 2 < den) return q; + if (r * 2 > den) return q + 1; + return (q & 1) ? q + 1 : q; + } + else return q; + } + // Round the non-negative offset quotient num/den (den >= 1) to an integer // notch index per L's rounding policy. // @@ -420,27 +439,8 @@ namespace bnd : static_cast(zl.Numerator)) : imax{0}; - // Offset rule: rounds num/den (≥ 0) directly, in q/r form so no formula - // can overflow umax (num + den/2 could, for num near umax). Used for - // exotic Lower/Notch (no integral value index) and as the fallback when - // the signed value-index rebuild below cannot fit 64 bits. - const auto offset_rule = [num, den]() -> umax - { - const umax q = num / den, r = num % den; - if constexpr (HasPolicy) return (r * 2 >= den) ? q + 1 : q; - else if constexpr (HasPolicy) return q; - else if constexpr (HasPolicy) return (r != 0) ? q + 1 : q; - else if constexpr (HasPolicy) - { - if (r * 2 < den) return q; - if (r * 2 > den) return q + 1; - return (q & 1) ? q + 1 : q; - } - else return q; - }; - if constexpr (!vidx) - return offset_rule(); + return round_offset(num / den, num % den, den); else { // Round the signed value-index NUM/di exactly like detail::div_rounded. @@ -452,7 +452,7 @@ namespace bnd if (num > static_cast(std::numeric_limits::max()) || mul_overflow(m, di, &mdi) || add_overflow(mdi, static_cast(num), &NUM)) [[unlikely]] - return offset_rule(); + return round_offset(num / den, num % den, den); const imax t = NUM / di; // C++ truncation toward zero const imax rr = NUM % di; // sign of NUM, |rr| < di imax J; diff --git a/single_include/bound/bound.hpp b/single_include/bound/bound.hpp index 630c52a..bb59447 100644 --- a/single_include/bound/bound.hpp +++ b/single_include/bound/bound.hpp @@ -2183,6 +2183,44 @@ namespace bnd::detail constexpr std::strong_ordering cmp128(u128 a, u128 b) { return (a.hi != b.hi) ? (a.hi <=> b.hi) : (a.lo <=> b.lo); } + // 128×64 product with an overflow flag (result beyond 128 bits). + struct mul128_result { u128 value; bool overflowed; }; + constexpr mul128_result mul128(u128 a, umax b) + { + const u128 low = umul(a.lo, b); + const u128 high = umul(a.hi, b); + const umax hi_sum = high.lo + low.hi; + return {u128{hi_sum, low.lo}, high.hi != 0 || hi_sum < low.hi}; + } + + // Quotient/remainder of a 128-bit dividend by a 64-bit divisor. Requires + // 1 <= d <= imax_max (the rational-denominator domain) so the portable + // partial remainder can never overflow when shifted. + struct divmod128_result { u128 quotient; umax remainder; }; + constexpr divmod128_result divmod128(u128 n, umax d) + { +#if defined(__SIZEOF_INT128__) + using u128n = unsigned __int128; + const u128n wide = (static_cast(n.hi) << 64) | n.lo; + const u128n q = wide / d; + return {u128{static_cast(q >> 64), static_cast(q)}, + static_cast(wide % d)}; +#else + // Portable (MSVC): restoring shift-subtract divide — the same construction + // as cmath.hpp's to_fixed fallback. + u128 q{0, 0}; + umax r = 0; + for (int i = 127; i >= 0; --i) + { + r = (r << 1) | ((i >= 64 ? (n.hi >> (i - 64)) : (n.lo >> i)) & 1u); + q.hi = (q.hi << 1) | (q.lo >> 63); + q.lo <<= 1; + if (r >= d) { r -= d; q.lo |= 1; } + } + return {q, r}; +#endif + } + //--------------------------------------------------------------------------- // trim //--------------------------------------------------------------------------- @@ -4519,6 +4557,25 @@ namespace bnd template inline constexpr bool HasPolicy = has_flag(BoundPolicy, F) || plain

::test(F); + // Rounds the split offset quotient q + r/den (r < den ≤ imax_max) per L's + // rounding policy — q/r form so no expression can overflow umax + // (num + den/2 could, for num near umax). Shared by round_quotient's + // offset rule and the 128-bit wide store (assignment.hpp). + template + [[nodiscard]] constexpr umax round_offset(umax q, umax r, umax den) noexcept + { + if constexpr (HasPolicy) return (r * 2 >= den) ? q + 1 : q; + else if constexpr (HasPolicy) return q; + else if constexpr (HasPolicy) return (r != 0) ? q + 1 : q; + else if constexpr (HasPolicy) + { + if (r * 2 < den) return q; + if (r * 2 > den) return q + 1; + return (q & 1) ? q + 1 : q; + } + else return q; + } + // Round the non-negative offset quotient num/den (den >= 1) to an integer // notch index per L's rounding policy. // @@ -4541,27 +4598,8 @@ namespace bnd : static_cast(zl.Numerator)) : imax{0}; - // Offset rule: rounds num/den (≥ 0) directly, in q/r form so no formula - // can overflow umax (num + den/2 could, for num near umax). Used for - // exotic Lower/Notch (no integral value index) and as the fallback when - // the signed value-index rebuild below cannot fit 64 bits. - const auto offset_rule = [num, den]() -> umax - { - const umax q = num / den, r = num % den; - if constexpr (HasPolicy) return (r * 2 >= den) ? q + 1 : q; - else if constexpr (HasPolicy) return q; - else if constexpr (HasPolicy) return (r != 0) ? q + 1 : q; - else if constexpr (HasPolicy) - { - if (r * 2 < den) return q; - if (r * 2 > den) return q + 1; - return (q & 1) ? q + 1 : q; - } - else return q; - }; - if constexpr (!vidx) - return offset_rule(); + return round_offset(num / den, num % den, den); else { // Round the signed value-index NUM/di exactly like detail::div_rounded. @@ -4573,7 +4611,7 @@ namespace bnd if (num > static_cast(std::numeric_limits::max()) || mul_overflow(m, di, &mdi) || add_overflow(mdi, static_cast(num), &NUM)) [[unlikely]] - return offset_rule(); + return round_offset(num / den, num % den, den); const imax t = NUM / di; // C++ truncation toward zero const imax rr = NUM % di; // sign of NUM, |rr| < di imax J; @@ -4988,6 +5026,92 @@ namespace bnd::detail action.fn(lhs, q); } + // 128-bit rounded store — the offset slot of an in-range rhs computed + // directly in wide arithmetic when the exact 64-bit rational formation + // of (rhs − Lower)/Notch overflows (full-mantissa fp-derived sources on + // grids with large |Lower|): + // slot + remainder/divisor = (rhs − Lower)·d_n / (a_dr·a_dl·n_n) + // The compile-time divisor factors are gcd-reduced first, so the only + // wide operations are one 128×64 multiply and one 128÷64 divide. + // ok == false when the reduced divisor or dividend exceeds the 128-bit + // envelope (or rhs is out of range — callers check range first). + struct wide_quotient { umax slot; umax remainder; umax divisor; bool ok; }; + + static constexpr wide_quotient wide_offset_quotient(rational const& rv) + { + if constexpr (Notch == 0) + return {}; // continuous grids never index slots + else + { + constexpr umax n_l = Lower.Numerator; + constexpr umax a_dl = abs_den(Lower.Denominator); + constexpr bool low_neg = Lower.Denominator < 0; + constexpr umax n_n = Notch.Numerator; // Notch > 0: d_n > 0 + constexpr umax d_n = static_cast(Notch.Denominator); + + // Compile-time divisor part; a grid whose a_dl·n_n cannot fit umax + // is beyond the wide envelope entirely. + constexpr umax den_ct = []{ + umax p; + return mul_overflow(a_dl, n_n, &p) ? umax{0} : p; + }(); + if constexpr (den_ct == 0) + return {}; + else + { + constexpr umax g1 = std::gcd(d_n, den_ct); + constexpr umax d_n1 = d_n / g1; + constexpr umax den1 = den_ct / g1; + + const umax a_dr = abs_den(rv.Denominator); + const bool rhs_neg = rv.Denominator < 0; + + // Offset numerator over the common denominator a_dr·a_dl: + // s_r·n_r·a_dl − s_l·n_l·a_dr (≥ 0 for in-range rhs). + // Each product is < 2^127 (numerator < 2^64, denominator ≤ imax), + // so the same-sign sum below cannot carry out of 128 bits. + const u128 val = umul(rv.Numerator, a_dl); + const u128 low = umul(n_l, a_dr); + u128 offset; + if (!rhs_neg && low_neg) + offset = u128{val.hi + low.hi + (val.lo + low.lo < val.lo ? 1u : 0u), + val.lo + low.lo}; + else if (!rhs_neg && !low_neg) + { + if (cmp128(val, low) < 0) return {}; // rhs < Lower + offset = u128{val.hi - low.hi - (val.lo < low.lo ? 1u : 0u), + val.lo - low.lo}; + } + else if (rhs_neg && low_neg) + { + if (cmp128(low, val) < 0) return {}; // rhs < Lower + offset = u128{low.hi - val.hi - (low.lo < val.lo ? 1u : 0u), + low.lo - val.lo}; + } + else + return {}; // rhs < 0 ≤ Lower + + const umax g2 = std::gcd(d_n1, a_dr); + const umax d_n2 = d_n1 / g2; + const umax a_dr1 = a_dr / g2; + + umax divisor; + if (mul_overflow(a_dr1, den1, &divisor) + || divisor > static_cast(std::numeric_limits::max())) + return {}; + + const mul128_result dividend = mul128(offset, d_n2); + if (dividend.overflowed) + return {}; + + const divmod128_result qr = divmod128(dividend.value, divisor); + if (qr.quotient.hi != 0) + return {}; // slot beyond any 64-bit index space + return {qr.quotient.lo, qr.remainder, divisor, true}; + } + } + } + template static constexpr bool store_checked(L& lhs, R rhs, P&& policy, A&& action = {}) { @@ -5076,16 +5200,36 @@ namespace bnd::detail } // The exact quotient can overflow the 64-bit rational range (huge - // source denominator × fine notch). Report it as an overflow through - // the policy channel instead of dereferencing nullopt — a throw here - // would escape noexcept callers (the math engines) as terminate. + // source denominator × fine notch). Recompute the slot directly in + // 128-bit (wide_offset_quotient above); only a result beyond even + // that envelope reports errc::overflow — never a nullopt deref, + // which would escape noexcept callers (the math engines) as + // terminate. Rounding here is the offset rule (round_offset), the + // same semantics round_quotient falls back to past 64 bits. const auto quotient = (rhs - Lower)/Notch; if (!quotient.has_value()) [[unlikely]] { - if constexpr (error_action>) - { action.fn(lhs, errc::overflow, errc_message(errc::overflow)); return false; } - policy.report(errc::overflow); - return false; + const wide_quotient wide = wide_offset_quotient(rational{rhs}); + if (!wide.ok) + { + if constexpr (error_action>) + { action.fn(lhs, errc::overflow, errc_message(errc::overflow)); return false; } + policy.report(errc::overflow); + return false; + } + if (wide.remainder == 0) + { store_slot(wide.slot); return true; } + if constexpr (has_round_flag) + { store_slot(round_offset(wide.slot, wide.remainder, wide.divisor)); return true; } + if (policy.round_check()) [[unlikely]] + { + if constexpr (error_action>) + { action.fn(lhs, errc::rounding_error, errc_message(errc::rounding_error)); return false; } + policy.report(errc::rounding_error); + return false; + } + store_slot(round_offset(wide.slot, wide.remainder, wide.divisor)); + return true; } rational raw = *quotient; umax den = static_cast(raw.Denominator); diff --git a/tests/accuracy.cpp b/tests/accuracy.cpp index 75f9e2c..c42db4f 100644 --- a/tests/accuracy.cpp +++ b/tests/accuracy.cpp @@ -89,9 +89,9 @@ namespace // Input grids — dyadic 2^-14 lattices with `real` (f64) storage, the fp // engines' intended pairing: results snap through the fp store path. All // engines run the identical grids, so notch-unit errors compare 1:1. - // (Integer-index snap grids currently route fp-engine stores through the - // 64-bit rational cold path, which overflows for full-mantissa results on - // grids with |Lower| ≫ 1 — see the store_checked overflow reporting.) + // (Integer-index snap grids also work — full-mantissa results store via the + // 128-bit rounded path — but `real` keeps the published table's store + // semantics uniform across engines.) using angle_grid = bound<{{-8, 8}, notch<1, 16384>}, round_nearest | real>; using tan_grid = bound<{{-1.5, 1.5}, notch<1, 16384>}, round_nearest | real>; using unit_grid = bound<{{-1, 1}, notch<1, 16384>}, round_nearest | real>; diff --git a/tests/test_cmath_double.cpp b/tests/test_cmath_double.cpp index b50f5cf..2408885 100644 --- a/tests/test_cmath_double.cpp +++ b/tests/test_cmath_double.cpp @@ -274,4 +274,26 @@ TEST_CASE("dbl: transcendental tier on real bounds matches std::", "[dbl][real][ } } + +//--------------------------------------------------------------------------- +// 2026-07 regression: full-mantissa engine results must store onto plain +// integer-index snap grids (not just `real` ones). tan's auto output grid +// has |Lower| ~ 1024, so a full-mantissa double result once overflowed the +// exact 64-bit (rhs - Lower)/Notch store and terminated through the noexcept +// engine; the 128-bit rounded store now lands the correctly rounded slot. +//--------------------------------------------------------------------------- +TEST_CASE("dbl engine stores full-mantissa results onto integer-index snap grids", + "[cmath][dbl][storage][regression]") +{ + using Ang = bound<{{-8, 8}, notch<1, 16384>}, round_nearest>; // integer-index storage + constexpr double half_notch = 0.5 / 16384.0; + + for (double x : {1.0, 1.5, -1.5, 0.4636}) + { + auto t = math::dbl::tan(Ang{x}); + REQUIRE(t.has_value()); + REQUIRE(std::fabs(static_cast(*t) - std::tan(x)) <= half_notch * 1.01); + } +} + #endif // !BND_MATH_FIXED diff --git a/tests/test_storage_bugs.cpp b/tests/test_storage_bugs.cpp index 890d0aa..87be47f 100644 --- a/tests/test_storage_bugs.cpp +++ b/tests/test_storage_bugs.cpp @@ -167,16 +167,17 @@ TEST_CASE("Bug F: real div-by-zero is reported, not a silent inf", // 2026-07: fp-derived rational store on a snap grid with |Lower| ≫ 1. The // cold store path forms (rhs − Lower)/Notch exactly; with a full-mantissa // double source (den 2^54) that once dereferenced a nullopt (terminate -// through the noexcept math engines). Now: the fitting case stores the -// correctly rounded slot, the truly unrepresentable case reports -// errc::overflow through the policy channel. +// through the noexcept math engines). Now: offsets that fit 64 bits go +// through the rescued rational path, and offsets beyond it go through the +// 128-bit rounded store (wide_offset_quotient) — both land on the correctly +// rounded slot. //--------------------------------------------------------------------------- -TEST_CASE("fp-derived rational store on a wide snap grid rounds or reports", +TEST_CASE("fp-derived rational store on a wide snap grid uses the 128-bit path", "[storage][overflow][regression]") { using wide = bound<{{-1024, 1024}, notch<1, 16384>}, round_nearest>; - SECTION("negative value: offset fits after the 128-bit rescue") + SECTION("negative value: offset fits after the 128-bit add rescue") { wide slot{}; slot = rational{umax{9006646171630191}, imax{-18014398509481984}}; // ≈ -0.4999693 @@ -184,14 +185,29 @@ TEST_CASE("fp-derived rational store on a wide snap grid rounds or reports", REQUIRE(rational{slot} == rational{umax{8191}, imax{-16384}}); } - SECTION("positive value: exact offset needs > 64 bits → errc::overflow, no terminate") + SECTION("positive value: exact offset needs > 64 bits → 128-bit rounded store") { wide slot{}; + slot = rational{umax{9006646171630191}, imax{18014398509481984}}; // ≈ +0.4999693 + // (1024 + v)·16384 = 16785407.49692… → round_nearest → slot 16785407 + // → value 8191/16384 (0.49993896…, the nearest grid point) + REQUIRE(rational{slot} == rational{umax{8191}, imax{16384}}); + } + + // (No constexpr section: at constant evaluation the transient rational + // overflow surfaces as the intentional constexpr_error build diagnostic + // before the wide fallback can engage — the 128-bit path is runtime-only + // in practice, though itself constexpr-capable.) + + SECTION("strict policy off-notch in the wide regime → rounding_error") + { + using strict = bound<{{-1024, 1024}, notch<1, 16384>}>; // checked, no round flag + strict slot{}; try { - slot = rational{umax{9006646171630191}, imax{18014398509481984}}; // ≈ +0.4999693 + slot = rational{umax{9006646171630191}, imax{18014398509481984}}; FAIL("expected the default handler to throw"); } - catch (bound_error const& e) { REQUIRE(e.code == errc::overflow); } + catch (bound_error const& e) { REQUIRE(e.code == errc::rounding_error); } } }