From c830b2cb2c06b2b3300a59dd4b39f9a1dc01475b Mon Sep 17 00:00:00 2001 From: Peter Neiss Date: Sat, 18 Jul 2026 18:51:22 +0200 Subject: [PATCH] Adopt prior-art implementation patterns; fix fp-store overflow terminate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Comparison pass over the resources.md libraries (bounded::integer, CNL, Boost.SafeNumerics, fpm, SafeInt, Intel safe-arithmetic, type_safe, PSsst, google/integers), adopting the patterns bound was missing: - tests/fail/: compile-fail diagnostics suite (bounded::integer/CNL) — each TU names its expected static_assert message in an EXPECT line; a control TU proves the harness builds valid code - tests/accuracy.cpp + docs/accuracy.md: fpm-style per-function accuracy tables for all three math engines against a long-double reference, regenerated by the accuracy_report target - detail/rep.hpp: fp representation propagation deduped out of addition/multiplication/division into one fp_rep trait (CNL single-concern style); codegen-neutral - grid hull + std::common_type + mixed-grid min/max (bounded::integer): hull(grid, grid) = interval hull + notch gcd, common_bound_t, and a std::common_type specialisation so mixed-grid bounds interoperate - docs: resources.md "implementation patterns compared" section (adopted/rejected with reasons), roadmap notes for type-level interval unions and the 128-bit rounded store The accuracy harness immediately found a real bug: storing a full-mantissa double-derived rational onto a snap grid with large |Lower| (e.g. dbl::tan's auto output grid) overflowed the exact 64-bit (rhs - Lower)/Notch rational and terminated through the noexcept engine as bad_optional_access. Fixed: - rational::add_impl rescues mixed-sign cross-product overflow in 128-bit when the difference fits umax - store_checked reports errc::overflow instead of dereferencing nullopt - round_quotient is now total: overflow-safe q/r offset rule replaces the historical formulas and serves as fallback when the signed value-index rebuild cannot fit 64 bits Same-sign exact offsets needing >64 bits still report errc::overflow; the full fix (compute the rounded index in 128-bit like to_fixed) is on the roadmap. Regression tests in test_rational.cpp and test_storage_bugs.cpp. Verified: 428/428 ctest (incl. compile-fail, amalgamate, codegen guard), 418/418 under BND_MATH_FIXED, fuzz seed 42 x 20k iters clean, clang-tidy clean under the CI header filter. Co-Authored-By: Claude Fable 5 --- CMakeLists.txt | 18 +- docs/accuracy.md | 61 ++++++ docs/math.md | 4 + docs/resources.md | 42 ++++ docs/roadmap.md | 15 ++ include/bound/arithmetic.hpp | 41 ++++ include/bound/detail/addition.hpp | 21 +- include/bound/detail/assignment.hpp | 14 +- include/bound/detail/division.hpp | 21 +- include/bound/detail/multiplication.hpp | 22 +- include/bound/detail/rational.hpp | 31 ++- include/bound/detail/rep.hpp | 46 +++++ include/bound/generic.hpp | 33 ++- include/bound/grid.hpp | 20 ++ include/bound/numeric_limits.hpp | 14 +- single_include/bound/bound.hpp | 255 ++++++++++++++++++------ tests/accuracy.cpp | 156 +++++++++++++++ tests/fail/CMakeLists.txt | 44 ++++ tests/fail/control_compiles.cpp | 12 ++ tests/fail/disjoint_interval_assign.cpp | 11 + tests/fail/exclusive_policies.cpp | 10 + tests/fail/math_without_snap.cpp | 12 ++ tests/fail/modulo_non_integer.cpp | 14 ++ tests/fail/notch_mismatch_assign.cpp | 14 ++ tests/fail/scalar_operand.cpp | 11 + tests/test_common_type.cpp | 122 ++++++++++++ tests/test_rational.cpp | 22 ++ tests/test_storage_bugs.cpp | 33 +++ 28 files changed, 988 insertions(+), 131 deletions(-) create mode 100644 docs/accuracy.md create mode 100644 include/bound/detail/rep.hpp create mode 100644 tests/accuracy.cpp create mode 100644 tests/fail/CMakeLists.txt create mode 100644 tests/fail/control_compiles.cpp create mode 100644 tests/fail/disjoint_interval_assign.cpp create mode 100644 tests/fail/exclusive_policies.cpp create mode 100644 tests/fail/math_without_snap.cpp create mode 100644 tests/fail/modulo_non_integer.cpp create mode 100644 tests/fail/notch_mismatch_assign.cpp create mode 100644 tests/fail/scalar_operand.cpp create mode 100644 tests/test_common_type.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 15191aa..fd5b608 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -158,6 +158,10 @@ target_include_directories(bound_tests SYSTEM PRIVATE ${CATCH2_INCLUDE_DIRS}) # PRE_TEST defers discovery to test invocation, under the configured emulator. catch_discover_tests(bound_tests DISCOVERY_MODE PRE_TEST) +# Compile-fail diagnostics suite — regression-tests that ill-formed uses stay +# ill-formed and keep their intended messages (see tests/fail/CMakeLists.txt). +add_subdirectory(tests/fail) + #--------------------------------------------------------------------------- # ctrack (benchmarking) #--------------------------------------------------------------------------- @@ -183,6 +187,18 @@ target_include_directories(bench SYSTEM PUBLIC ${CTRACK_INCLUDE_DIRS}) add_executable(bound_fuzz EXCLUDE_FROM_ALL tests/fuzz.cpp) target_link_libraries(bound_fuzz PUBLIC bound::bound) +# Accuracy characterization (fpm-style): sweeps every math engine against a +# long-double reference and emits a markdown error table. On demand: +# `--target bound_accuracy` to build, `--target accuracy_report` to regenerate +# the committed docs/accuracy.md. +add_executable(bound_accuracy EXCLUDE_FROM_ALL tests/accuracy.cpp) +target_link_libraries(bound_accuracy PUBLIC bound::bound) + +add_custom_target(accuracy_report + COMMAND bound_accuracy ${CMAKE_CURRENT_SOURCE_DIR}/docs/accuracy.md + COMMENT "Regenerating docs/accuracy.md") +add_dependencies(accuracy_report bound_accuracy) + #--------------------------------------------------------------------------- # MinGW link workaround #--------------------------------------------------------------------------- @@ -192,7 +208,7 @@ target_link_libraries(bound_fuzz PUBLIC bound::bound) # "multiple definition of std::type_info::operator==". Tolerate it on MinGW # (the two definitions are identical); harmless elsewhere, so it stays scoped. if(WIN32 AND CMAKE_CXX_COMPILER_ID STREQUAL "GNU") - foreach(tgt bound_tests bench bound_fuzz) + foreach(tgt bound_tests bench bound_fuzz bound_accuracy) target_link_options(${tgt} PRIVATE -Wl,--allow-multiple-definition) endforeach() endif() diff --git a/docs/accuracy.md b/docs/accuracy.md new file mode 100644 index 0000000..2a40002 --- /dev/null +++ b/docs/accuracy.md @@ -0,0 +1,61 @@ +# Engine accuracy + +Generated by `bound_accuracy` (tests/accuracy.cpp) — do not edit by hand. +Regenerate with `cmake --build --target accuracy_report`. + +Error is measured against a long-double `` reference, in units of +the **output grid notch**: max < 0.5 means every sampled result is the +correctly rounded grid point. Engines: `cordic` (integer/CORDIC, FPU-free, +bit-exact everywhere), `dbl` (reproducible binary64), `flt` (binary32). +See [math.md](math.md) for the engine model. + +| fn | engine | domain | max (notch) | mean (notch) | samples | +|---------|--------|----------------|-------------|--------------|---------| +| sin | cordic | [-8, 8] | 0.7516 | 0.263837 | 20165 | +| sin | dbl | [-8, 8] | 0.5000 | 0.249295 | 20165 | +| sin | flt | [-8, 8] | 0.5004 | 0.249295 | 20165 | +| cos | cordic | [-8, 8] | 0.7423 | 0.260149 | 20165 | +| cos | dbl | [-8, 8] | 0.5000 | 0.248520 | 20165 | +| cos | flt | [-8, 8] | 0.5005 | 0.248520 | 20165 | +| tan | cordic | [-1.5, 1.5] | 0.5006 | 0.247359 | 24577 | +| tan | dbl | [-1.5, 1.5] | 0.5000 | 0.247359 | 24577 | +| tan | flt | [-1.5, 1.5] | 0.5178 | 0.247367 | 24577 | +| asin | cordic | [-1, 1] | 0.5001 | 0.243364 | 32769 | +| asin | dbl | [-1, 1] | 0.5000 | 0.243364 | 32769 | +| asin | flt | [-1, 1] | 0.5015 | 0.243365 | 32769 | +| acos | cordic | [-1, 1] | 0.5000 | 0.244882 | 32769 | +| acos | dbl | [-1, 1] | 0.5000 | 0.244882 | 32769 | +| acos | flt | [-1, 1] | 0.5035 | 0.244887 | 32769 | +| atan | cordic | [-8, 8] | 0.5218 | 0.248050 | 20165 | +| atan | dbl | [-8, 8] | 0.5000 | 0.247945 | 20165 | +| atan | flt | [-8, 8] | 0.5018 | 0.247946 | 20165 | +| sinh | cordic | [-8, 8] | 0.8942 | 0.248963 | 20165 | +| sinh | dbl | [-8, 8] | 0.5000 | 0.247560 | 20165 | +| sinh | flt | [-8, 8] | 3.1212 | 0.298317 | 20165 | +| cosh | cordic | [-8, 8] | 0.8530 | 0.253277 | 20165 | +| cosh | dbl | [-8, 8] | 0.5000 | 0.251411 | 20165 | +| cosh | flt | [-8, 8] | 2.7390 | 0.297414 | 20165 | +| tanh | cordic | [-8, 8] | 0.5000 | 0.204307 | 20165 | +| tanh | dbl | [-8, 8] | 0.5000 | 0.204307 | 20165 | +| tanh | flt | [-8, 8] | 0.5007 | 0.204308 | 20165 | +| exp | cordic | [-8, 8] | 0.8889 | 0.250463 | 20165 | +| exp | dbl | [-8, 8] | 0.5000 | 0.249378 | 20165 | +| exp | flt | [-8, 8] | 2.7675 | 0.284219 | 20165 | +| exp2 | cordic | [-8, 8] | 0.6270 | 0.250352 | 20165 | +| exp2 | dbl | [-8, 8] | 0.5000 | 0.250156 | 20165 | +| exp2 | flt | [-8, 8] | 1.5950 | 0.265997 | 20165 | +| log | cordic | [1, 1000] | 0.5234 | 0.249774 | 20010 | +| log | dbl | [1, 1000] | 0.5000 | 0.249657 | 20010 | +| log | flt | [1, 1000] | 0.5044 | 0.249671 | 20010 | +| log2 | cordic | [1, 1000] | 0.5186 | 0.248581 | 20010 | +| log2 | dbl | [1, 1000] | 0.5000 | 0.248472 | 20010 | +| log2 | flt | [1, 1000] | 0.5106 | 0.248516 | 20010 | +| log10 | cordic | [1, 1000] | 0.5000 | 0.252197 | 20010 | +| log10 | dbl | [1, 1000] | 0.5000 | 0.252197 | 20010 | +| log10 | flt | [1, 1000] | 0.5045 | 0.252205 | 20010 | +| sqrt | cordic | [0, 16] | 0.5087 | 0.248846 | 20165 | +| sqrt | dbl | [0, 16] | 0.5000 | 0.248835 | 20165 | +| sqrt | flt | [0, 16] | 0.5019 | 0.248838 | 20165 | +| cbrt | cordic | [-8, 8] | 0.5004 | 0.250647 | 20165 | +| cbrt | dbl | [-8, 8] | 0.5000 | 0.250647 | 20165 | +| cbrt | flt | [-8, 8] | 0.5016 | 0.250649 | 20165 | diff --git a/docs/math.md b/docs/math.md index c109431..28f159b 100644 --- a/docs/math.md +++ b/docs/math.md @@ -117,6 +117,10 @@ other representation flags. output grid. - Algebraically-exact results (e.g. `cbrt(8)`, `hypot(3,4)`, `pow(2,10)`) land exactly under both engines. + - Measured per-function error tables for all three engines (max/mean in + output-notch units, against a long-double reference) are in + [accuracy.md](accuracy.md), regenerated by the `accuracy_report` build + target. - **Input-range limits** below are engine-shared `static_assert` envelopes, kept identical for both engines so the same programs compile everywhere. The trig/root/atan limits (±2^20) come from the integer engine's working diff --git a/docs/resources.md b/docs/resources.md index fe3f234..a842550 100644 --- a/docs/resources.md +++ b/docs/resources.md @@ -87,6 +87,48 @@ column. Pick the tool that fits: if you need overflow-safe plain integers, `bounded::integer` or Boost.SafeNumerics are proven; for pure fixed-point DSP, `fpm` or CNL; `bound` is for when the *domain* of a value is part of its type. +## Implementation patterns compared + +A 2026-07 pass over these libraries' *implementation* patterns (not just their +feature sets), asking what `bound` should adopt. Kept, with the source named: + +- **Compile-fail diagnostics tests** (bounded::integer, CNL) — `tests/fail/` + regression-tests that ill-formed uses stay ill-formed *and* keep their + intended `static_assert` message (each TU names its expected diagnostic in a + `// EXPECT:` line; a control TU proves the harness builds valid code). +- **Published accuracy characterization** (fpm) — fpm ships per-function + accuracy tables for its math; `bound` now does too: + [accuracy.md](accuracy.md), generated by `tests/accuracy.cpp` sweeping all + three engines against a long-double reference. (Its first run found a real + bug: a 64-bit overflow terminate in the fp→snap-grid store path.) +- **`std::common_type` interop** (bounded::integer) — mixed-grid bounds now + have a common type (the grid hull: interval hull + notch gcd), plus + mixed-grid `bnd::min`/`max` returning it. +- **Single-concern helpers** (CNL's component style) — the per-operator fp + representation-propagation logic, once copy-pasted across + addition/multiplication/division, is one trait (`detail/rep.hpp`). + +Considered and **not** adopted, with reasons: + +- **CRTP operator mix-ins** (PSsst, and Boost's operator-generation style) — + `bound`'s operators dispatch on storage shape with `if constexpr`, which a + generic mix-in base can't see; hand-written operators stay. +- **Selectable promotion policies** (Boost.SafeNumerics `native`/`automatic`) — + `bound` always auto-widens the result grid; a "native" mode would reintroduce + exactly the silent narrowing the library exists to prevent. +- **Constraint/verifier splitting** (type_safe) — the policy bitmask plus + tagged `on_*` actions already separate *what* is enforced from *what happens* + on violation, without a second customization axis. +- **Type-per-concern wrapper stacking** (CNL) — composing + `overflow>>`-style layers trades one flat policy set for + N interacting wrapper semantics; `bound`'s single `policy_flag` set keeps one + place to reason about flag interaction (mutual exclusions are + `static_assert`ed in one spot). +- **Disjoint interval unions** (Intel safe-arithmetic) — genuinely attractive + (zero-excluding divisors make division total); parked as a + [roadmap sketch](roadmap.md) because every operator and predicate would need + a union story. + ## Related libraries - **[bounded::integer](https://github.com/davidstone/bounded-integer)** (David diff --git a/docs/roadmap.md b/docs/roadmap.md index 151f6c9..829fb38 100644 --- a/docs/roadmap.md +++ b/docs/roadmap.md @@ -81,6 +81,21 @@ For completeness — these came up alongside the above but are *not* blocked by with no language dependency. See [freestanding.md](freestanding.md). - **`dyn_bound`** (runtime-valued bounds) — evaluated and **declined** on design grounds (no space/time-efficient implementation), not deferred. +- **Type-level interval unions** — Intel's [safe-arithmetic](resources.md) models + *disjoint* interval unions in the type (`ival<-1000,-1> || ival<1,1000>`), so e.g. a + zero-excluding divisor makes division provably total at compile time. `bound`'s + `grid::operator/` already computes the two zero-free halves internally + (`grid.hpp`) — the missing piece is expressing the union in the *type* (a + `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. - **Modules / compile-time-footprint work** — parked for later; modules is C++20, not a blocker. diff --git a/include/bound/arithmetic.hpp b/include/bound/arithmetic.hpp index a876fab..d783221 100644 --- a/include/bound/arithmetic.hpp +++ b/include/bound/arithmetic.hpp @@ -6,6 +6,7 @@ #include "bound/core.hpp" #include "bound/casts.hpp" +#include "bound/detail/rep.hpp" #include #include @@ -245,6 +246,36 @@ namespace bnd boundable auto t) { return a + (b - a) * t; } + //--------------------------------------------------------------------------- + // common_bound — the "hull" type able to hold every value of L and R exactly: + // interval hull + notch gcd (grid `hull`), representation propagated by the + // same widest-wins rule as arithmetic results (detail::fp_rep). Backs the + // std::common_type specialisation (numeric_limits.hpp) and mixed-grid + // min/max below. The primary has no `type` when the hull grid is + // unrepresentable, so common_type_t SFINAEs away instead of erroring. + //--------------------------------------------------------------------------- + namespace detail + { + template + struct common_bound {}; + + // Same type stays itself (policy included) — mirrors std::common_type. + template + struct common_bound { using type = Same; }; + + template + requires (!std::same_as) && (hull(Grid, Grid).has_value()) + struct common_bound + { + static constexpr grid hull_grid = *hull(Grid, Grid); + using type = bound::result_policy>; + }; + } + + template + using common_bound_t = typename detail::common_bound::type; + //--------------------------------------------------------------------------- // std-vocabulary helpers — ADL-found `min` / `max` / `midpoint` for generic // code. min/max mirror std; midpoint returns the *exact* average on a refined @@ -257,6 +288,16 @@ namespace bnd template [[nodiscard]] constexpr T max(T a, T b) { return (a < b) ? b : a; } + // Mixed-grid forms return the common hull type (both operands convert + // losslessly — the hull is assignable from each by construction). + template requires (!std::same_as) + [[nodiscard]] constexpr auto min(Lhs a, Rhs b) -> common_bound_t + { common_bound_t ca{a}, cb{b}; return (cb < ca) ? cb : ca; } + + template requires (!std::same_as) + [[nodiscard]] constexpr auto max(Lhs a, Rhs b) -> common_bound_t + { common_bound_t ca{a}, cb{b}; return (ca < cb) ? cb : ca; } + template [[nodiscard]] constexpr auto midpoint(T a, T b) { return (a + b) * just>; } diff --git a/include/bound/detail/addition.hpp b/include/bound/detail/addition.hpp index 82e51e7..aff3c41 100644 --- a/include/bound/detail/addition.hpp +++ b/include/bound/detail/addition.hpp @@ -4,6 +4,7 @@ #ifndef BNDadditionHPP #define BNDadditionHPP +#include "bound/detail/rep.hpp" #include "bound/generic.hpp" #include "bound/grid.hpp" #include "bound/policy.hpp" @@ -24,23 +25,9 @@ namespace bnd::detail "addition: result grid's notch/interval exceeds the representable rational " "range — coarsen the operand grids"); static constexpr grid result_grid = (Grid + Grid).value(); - // Propagate fp storage only when the result grid stays exactly representable - // in the chosen width; otherwise demote (f32→f64) or drop it so storage_pick - // deduces an exact representation (the fp sum would diverge from the exact sum - // — see grid::double_exact / float_exact). Widest-wins: prefer f32 only when - // both operands are f32-only and the result fits float; an f64 operand or a - // too-fine-for-float result widens to f64; too fine for double → exact. - static constexpr bool any_f64 = - (BoundPolicy & bnd::real) == bnd::real || (BoundPolicy & bnd::real) == bnd::real; - static constexpr bool any_f32 = - (BoundPolicy & bnd::f32) == bnd::f32 || (BoundPolicy & bnd::f32) == bnd::f32; - static constexpr bool keep_f32 = any_f32 && !any_f64 && float_exact; - static constexpr bool keep_f64 = !keep_f32 && (any_f64 || any_f32) && double_exact; - // Carry both operands' representation flags (widest-wins at storage selection). - static constexpr policy_flag rep = - ((BoundPolicy | BoundPolicy) & (bnd::exact | bnd::direct | bnd::indexed)) - | (keep_f64 ? bnd::real : none) | (keep_f32 ? bnd::f32 : none); - using result = bound; + // fp / representation propagation — shared rule in detail/rep.hpp. + using rep_t = fp_rep; + using result = bound; template static constexpr bool needs_overflow_check = diff --git a/include/bound/detail/assignment.hpp b/include/bound/detail/assignment.hpp index bbd6fae..e66cd3b 100644 --- a/include/bound/detail/assignment.hpp +++ b/include/bound/detail/assignment.hpp @@ -345,7 +345,19 @@ namespace bnd::detail } } - rational raw = ((rhs - Lower)/Notch).value(); + // 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. + 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; + } + rational raw = *quotient; umax den = static_cast(raw.Denominator); if (den == 1) { store_slot(raw.Numerator); return true; } diff --git a/include/bound/detail/division.hpp b/include/bound/detail/division.hpp index 64f3571..7230762 100644 --- a/include/bound/detail/division.hpp +++ b/include/bound/detail/division.hpp @@ -4,6 +4,7 @@ #ifndef BNDdivisionHPP #define BNDdivisionHPP +#include "bound/detail/rep.hpp" #include "bound/generic.hpp" #include "bound/grid.hpp" #include "bound/policy.hpp" @@ -154,22 +155,10 @@ namespace bnd::detail ? grid{interval{rational{0}, (Upper / Notch).value()}, Notch} : *(Grid / Grid); - static constexpr bool any_f64 = - (BoundPolicy & bnd::real) == bnd::real || (BoundPolicy & bnd::real) == bnd::real; - static constexpr bool any_f32 = - (BoundPolicy & bnd::f32) == bnd::f32 || (BoundPolicy & bnd::f32) == bnd::f32; - // Keep fp for a continuous result (Notch 0: the raw stores the quotient - // verbatim) or an fp-exact dyadic result; otherwise demote f32→f64 / drop f64 - // (the fp quotient would not land on the result grid). Widest-wins as in mul. - static constexpr bool keep_f32 = - any_f32 && !any_f64 && (result_grid.Notch == 0 || float_exact); - static constexpr bool keep_f64 = - !keep_f32 && (any_f64 || any_f32) && (result_grid.Notch == 0 || double_exact); - // Carry both operands' representation flags (widest-wins at storage selection). - static constexpr policy_flag rep = - ((BoundPolicy | BoundPolicy) & (bnd::exact | bnd::direct | bnd::indexed)) - | (keep_f64 ? bnd::real : none) | (keep_f32 ? bnd::f32 : none); - using result = bound; + // fp / representation propagation — shared rule in detail/rep.hpp. + // AllowContinuous: a continuous quotient (Notch 0) keeps fp verbatim. + using rep_t = fp_rep; + using result = bound; template static constexpr bool needs_overflow_check = diff --git a/include/bound/detail/multiplication.hpp b/include/bound/detail/multiplication.hpp index d0a2a9c..57681fb 100644 --- a/include/bound/detail/multiplication.hpp +++ b/include/bound/detail/multiplication.hpp @@ -4,6 +4,7 @@ #ifndef BNDmultiplicationHPP #define BNDmultiplicationHPP +#include "bound/detail/rep.hpp" #include "bound/generic.hpp" #include "bound/grid.hpp" #include "bound/policy.hpp" @@ -24,22 +25,11 @@ namespace bnd::detail "multiplication: result grid's notch/interval exceeds the representable " "rational range — coarsen the operand grids"); static constexpr grid result_grid = (Grid * Grid).value(); - // fp storage propagation (see addition.hpp): the product grid (notch = N_L·N_R) - // is finer, so demote f32→f64 / drop f64 when it outgrows the width — the fp - // product would round below the result notch. Widest-wins: f32 only if both - // operands f32-only and product fits float; else widen to f64; else exact. - static constexpr bool any_f64 = - (BoundPolicy & bnd::real) == bnd::real || (BoundPolicy & bnd::real) == bnd::real; - static constexpr bool any_f32 = - (BoundPolicy & bnd::f32) == bnd::f32 || (BoundPolicy & bnd::f32) == bnd::f32; - static constexpr bool keep_f32 = any_f32 && !any_f64 && float_exact; - static constexpr bool keep_f64 = !keep_f32 && (any_f64 || any_f32) && double_exact; - static constexpr bool dropped_fp = (any_f64 || any_f32) && !keep_f64 && !keep_f32; - // Carry both operands' representation flags (widest-wins at storage selection). - static constexpr policy_flag rep = - ((BoundPolicy | BoundPolicy) & (bnd::exact | bnd::direct | bnd::indexed)) - | (keep_f64 ? bnd::real : none) | (keep_f32 ? bnd::f32 : none); - using result = bound; + // fp / representation propagation — shared rule in detail/rep.hpp. The product + // grid (notch = N_L·N_R) is finer, so demotion/dropping is the common case. + using rep_t = fp_rep; + static constexpr bool dropped_fp = rep_t::dropped_fp; + using result = bound; // The dropped-fp case lands on a rational result when the product grid outgrows // uint index space; its product numerator can exceed `umax`, so check it (the diff --git a/include/bound/detail/rational.hpp b/include/bound/detail/rational.hpp index 626faad..c2cf9c0 100644 --- a/include/bound/detail/rational.hpp +++ b/include/bound/detail/rational.hpp @@ -585,10 +585,37 @@ namespace bnd::detail if constexpr (Checked) { if (mul_overflow(a_ad, b_ad_r, &denominator) || // = lcm(a_ad, b_ad) - mul_overflow(a.Numerator, b_ad_r, &A) || - mul_overflow(b.Numerator, a_ad_r, &B) || denominator > static_cast(std::numeric_limits::max())) { + if (std::is_constant_evaluated()) { constexpr_error<"rational +: denominator overflow">(); } + return ret_t{slim::nullopt}; + } + if (mul_overflow(a.Numerator, b_ad_r, &A) || + mul_overflow(b.Numerator, a_ad_r, &B)) + { + // Mixed signs: |A − B| can fit umax even when a cross-product alone + // does not (e.g. 1024 − m/2^54 forms 1024·2^54 == 2^64 before the + // subtraction brings it back in range — the dbl-engine store path hit + // exactly this). Retry the difference in 128-bit before giving up. + if (a_neg != b_neg) + { + const u128 A128 = umul(a.Numerator, b_ad_r); + const u128 B128 = umul(b.Numerator, a_ad_r); + const bool a_bigger = cmp128(A128, B128) > 0; + const u128 big = a_bigger ? A128 : B128; + const u128 small = a_bigger ? B128 : A128; + const u128 diff{big.hi - small.hi - (big.lo < small.lo ? 1u : 0u), + big.lo - small.lo}; + if (diff.hi == 0) + { + rational r; + r.Numerator = diff.lo; + r.Denominator = (a_neg ? a_bigger : !a_bigger) ? -denominator + : denominator; + trim(r.Numerator, r.Denominator); + return ret_t{r}; + } + } if (std::is_constant_evaluated()) { constexpr_error<"rational +: cross-multiplication overflow">(); } return ret_t{slim::nullopt}; } diff --git a/include/bound/detail/rep.hpp b/include/bound/detail/rep.hpp new file mode 100644 index 0000000..a9ea9e7 --- /dev/null +++ b/include/bound/detail/rep.hpp @@ -0,0 +1,46 @@ +//--------------------------------------------------------------------------- +// Copyright (C) 2026 Peter Neiss +//--------------------------------------------------------------------------- +#ifndef BNDrepHPP +#define BNDrepHPP + +#include "bound/generic.hpp" +#include "bound/grid.hpp" +#include "bound/policy_flag.hpp" + +//--------------------------------------------------------------------------- +// rep — representation-flag propagation for binary arithmetic results, shared +// by addition/multiplication/division. Propagate fp storage only when the +// result grid stays exactly representable in the chosen width; otherwise +// demote (f32→f64) or drop it so storage_pick deduces an exact representation +// (the fp result would diverge from the exact result — see grid::double_exact +// / float_exact). Widest-wins: prefer f32 only when both operands are +// f32-only and the result fits float; an f64 operand or a too-fine-for-float +// result widens to f64; too fine for double → exact. Division sets +// AllowContinuous: a continuous result (Notch 0) keeps fp regardless — the +// raw stores the quotient verbatim, so there is no grid to land on. +//--------------------------------------------------------------------------- +namespace bnd::detail +{ + template + struct fp_rep + { + static constexpr bool any_f64 = + (BoundPolicy & bnd::real) == bnd::real || (BoundPolicy & bnd::real) == bnd::real; + static constexpr bool any_f32 = + (BoundPolicy & bnd::f32) == bnd::f32 || (BoundPolicy & bnd::f32) == bnd::f32; + static constexpr bool continuous_ok = AllowContinuous && ResultGrid.Notch == 0; + static constexpr bool keep_f32 = + any_f32 && !any_f64 && (continuous_ok || float_exact); + static constexpr bool keep_f64 = + !keep_f32 && (any_f64 || any_f32) && (continuous_ok || double_exact); + static constexpr bool dropped_fp = (any_f64 || any_f32) && !keep_f64 && !keep_f32; + // Carry both operands' representation flags (widest-wins at storage selection). + static constexpr policy_flag rep = + ((BoundPolicy | BoundPolicy) & (bnd::exact | bnd::direct | bnd::indexed)) + | (keep_f64 ? bnd::real : none) | (keep_f32 ? bnd::f32 : none); + // The result bound's policy: the propagated representation, or plain checked. + static constexpr policy_flag result_policy = rep != none ? rep : checked; + }; +} +#endif diff --git a/include/bound/generic.hpp b/include/bound/generic.hpp index d6a1c77..5e0aca2 100644 --- a/include/bound/generic.hpp +++ b/include/bound/generic.hpp @@ -420,26 +420,39 @@ namespace bnd : static_cast(zl.Numerator)) : imax{0}; - if constexpr (!vidx) + // 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 { - // Exotic Lower/Notch (no integral value index) — historical offset rule. - if constexpr (HasPolicy) return (num + den / 2) / den; - else if constexpr (HasPolicy) return num / den; - else if constexpr (HasPolicy) return (num + den - 1) / den; + 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) { - umax q = num / den, r = num % den; if (r * 2 < den) return q; if (r * 2 > den) return q + 1; return (q & 1) ? q + 1 : q; } - else return num / den; - } + else return q; + }; + + if constexpr (!vidx) + return offset_rule(); else { // Round the signed value-index NUM/di exactly like detail::div_rounded. - const imax di = static_cast(den); - const imax NUM = m * di + static_cast(num); + // A numerator or m·di beyond imax (fp-derived sources on grids with + // large |Lower·count|) cannot rebuild the signed index — fall back to + // the offset rule, which differs only at exact ties on negative values. + const imax di = static_cast(den); + imax mdi, NUM; + if (num > static_cast(std::numeric_limits::max()) + || mul_overflow(m, di, &mdi) + || add_overflow(mdi, static_cast(num), &NUM)) [[unlikely]] + return offset_rule(); const imax t = NUM / di; // C++ truncation toward zero const imax rr = NUM % di; // sign of NUM, |rr| < di imax J; diff --git a/include/bound/grid.hpp b/include/bound/grid.hpp index 32f62ac..c455ba5 100644 --- a/include/bound/grid.hpp +++ b/include/bound/grid.hpp @@ -463,6 +463,26 @@ namespace bnd lhs.Interval / interval{rhs.Interval.Lower, -step}); } } + + //--------------------------------------------------------------------------- + // hull + //--------------------------------------------------------------------------- + // The smallest grid that represents every value of both operands exactly: + // interval hull + notch gcd. A valid grid anchors Lower on a multiple of its + // notch, so both lattices are sub-lattices of the gcd lattice — no offset + // term is needed, and the hull is a valid grid by construction. A continuous + // operand (Notch 0) makes the hull continuous. nullopt when the notch gcd's + // combined denominator exceeds the representable rational range. + //--------------------------------------------------------------------------- + inline constexpr slim::optional hull(const grid& lhs, const grid& rhs) + { + const interval iv{std::min(lhs.Interval.Lower, rhs.Interval.Lower), + std::max(lhs.Interval.Upper, rhs.Interval.Upper)}; + if (lhs.Notch == 0 || rhs.Notch == 0) + return grid{iv, detail::rational{0}}; + return lift([iv](detail::rational g){ return grid{iv, g}; }, + detail::gcd(lhs.Notch, rhs.Notch)); + } } // namespace bnd namespace slim diff --git a/include/bound/numeric_limits.hpp b/include/bound/numeric_limits.hpp index 5d61368..bfe3bec 100644 --- a/include/bound/numeric_limits.hpp +++ b/include/bound/numeric_limits.hpp @@ -8,13 +8,21 @@ #include #include +#include // std::common_type //--------------------------------------------------------------------------- -// numeric_limits / hash — std:: specialisations for bound. numeric_limits -// reports the *grid* bounds (Lower/Upper), not the raw type's limits. std::hash -// hashes the Raw member (rational raw: Numerator+Denominator, boost-style combine). +// numeric_limits / hash / common_type — std:: specialisations for bound. +// numeric_limits reports the *grid* bounds (Lower/Upper), not the raw type's +// limits. std::hash hashes the Raw member (rational raw: Numerator+Denominator, +// boost-style combine). std::common_type is the grid hull (see bnd::common_bound +// in arithmetic.hpp) so mixed-grid bounds interoperate with generic code. //--------------------------------------------------------------------------- +template +struct std::common_type, bnd::bound> + : bnd::detail::common_bound, bnd::bound> {}; + + template struct std::numeric_limits> { diff --git a/single_include/bound/bound.hpp b/single_include/bound/bound.hpp index 049d17e..630c52a 100644 --- a/single_include/bound/bound.hpp +++ b/single_include/bound/bound.hpp @@ -2718,10 +2718,37 @@ namespace bnd::detail if constexpr (Checked) { if (mul_overflow(a_ad, b_ad_r, &denominator) || // = lcm(a_ad, b_ad) - mul_overflow(a.Numerator, b_ad_r, &A) || - mul_overflow(b.Numerator, a_ad_r, &B) || denominator > static_cast(std::numeric_limits::max())) { + if (std::is_constant_evaluated()) { constexpr_error<"rational +: denominator overflow">(); } + return ret_t{slim::nullopt}; + } + if (mul_overflow(a.Numerator, b_ad_r, &A) || + mul_overflow(b.Numerator, a_ad_r, &B)) + { + // Mixed signs: |A − B| can fit umax even when a cross-product alone + // does not (e.g. 1024 − m/2^54 forms 1024·2^54 == 2^64 before the + // subtraction brings it back in range — the dbl-engine store path hit + // exactly this). Retry the difference in 128-bit before giving up. + if (a_neg != b_neg) + { + const u128 A128 = umul(a.Numerator, b_ad_r); + const u128 B128 = umul(b.Numerator, a_ad_r); + const bool a_bigger = cmp128(A128, B128) > 0; + const u128 big = a_bigger ? A128 : B128; + const u128 small = a_bigger ? B128 : A128; + const u128 diff{big.hi - small.hi - (big.lo < small.lo ? 1u : 0u), + big.lo - small.lo}; + if (diff.hi == 0) + { + rational r; + r.Numerator = diff.lo; + r.Denominator = (a_neg ? a_bigger : !a_bigger) ? -denominator + : denominator; + trim(r.Numerator, r.Denominator); + return ret_t{r}; + } + } if (std::is_constant_evaluated()) { constexpr_error<"rational +: cross-multiplication overflow">(); } return ret_t{slim::nullopt}; } @@ -4055,6 +4082,26 @@ namespace bnd lhs.Interval / interval{rhs.Interval.Lower, -step}); } } + + //--------------------------------------------------------------------------- + // hull + //--------------------------------------------------------------------------- + // The smallest grid that represents every value of both operands exactly: + // interval hull + notch gcd. A valid grid anchors Lower on a multiple of its + // notch, so both lattices are sub-lattices of the gcd lattice — no offset + // term is needed, and the hull is a valid grid by construction. A continuous + // operand (Notch 0) makes the hull continuous. nullopt when the notch gcd's + // combined denominator exceeds the representable rational range. + //--------------------------------------------------------------------------- + inline constexpr slim::optional hull(const grid& lhs, const grid& rhs) + { + const interval iv{std::min(lhs.Interval.Lower, rhs.Interval.Lower), + std::max(lhs.Interval.Upper, rhs.Interval.Upper)}; + if (lhs.Notch == 0 || rhs.Notch == 0) + return grid{iv, detail::rational{0}}; + return lift([iv](detail::rational g){ return grid{iv, g}; }, + detail::gcd(lhs.Notch, rhs.Notch)); + } } // namespace bnd namespace slim @@ -4494,26 +4541,39 @@ namespace bnd : static_cast(zl.Numerator)) : imax{0}; - if constexpr (!vidx) + // 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 { - // Exotic Lower/Notch (no integral value index) — historical offset rule. - if constexpr (HasPolicy) return (num + den / 2) / den; - else if constexpr (HasPolicy) return num / den; - else if constexpr (HasPolicy) return (num + den - 1) / den; + 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) { - umax q = num / den, r = num % den; if (r * 2 < den) return q; if (r * 2 > den) return q + 1; return (q & 1) ? q + 1 : q; } - else return num / den; - } + else return q; + }; + + if constexpr (!vidx) + return offset_rule(); else { // Round the signed value-index NUM/di exactly like detail::div_rounded. - const imax di = static_cast(den); - const imax NUM = m * di + static_cast(num); + // A numerator or m·di beyond imax (fp-derived sources on grids with + // large |Lower·count|) cannot rebuild the signed index — fall back to + // the offset rule, which differs only at exact ties on negative values. + const imax di = static_cast(den); + imax mdi, NUM; + if (num > static_cast(std::numeric_limits::max()) + || mul_overflow(m, di, &mdi) + || add_overflow(mdi, static_cast(num), &NUM)) [[unlikely]] + return offset_rule(); const imax t = NUM / di; // C++ truncation toward zero const imax rr = NUM % di; // sign of NUM, |rr| < di imax J; @@ -5015,7 +5075,19 @@ namespace bnd::detail } } - rational raw = ((rhs - Lower)/Notch).value(); + // 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. + 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; + } + rational raw = *quotient; umax den = static_cast(raw.Denominator); if (den == 1) { store_slot(raw.Numerator); return true; } @@ -5656,6 +5728,49 @@ namespace bnd //--------------------------------------------------------------------------- +// ====================================================================== +// bound/detail/rep.hpp +// ====================================================================== +//--------------------------------------------------------------------------- +//--------------------------------------------------------------------------- + + +//--------------------------------------------------------------------------- +// rep — representation-flag propagation for binary arithmetic results, shared +// by addition/multiplication/division. Propagate fp storage only when the +// result grid stays exactly representable in the chosen width; otherwise +// demote (f32→f64) or drop it so storage_pick deduces an exact representation +// (the fp result would diverge from the exact result — see grid::double_exact +// / float_exact). Widest-wins: prefer f32 only when both operands are +// f32-only and the result fits float; an f64 operand or a too-fine-for-float +// result widens to f64; too fine for double → exact. Division sets +// AllowContinuous: a continuous result (Notch 0) keeps fp regardless — the +// raw stores the quotient verbatim, so there is no grid to land on. +//--------------------------------------------------------------------------- +namespace bnd::detail +{ + template + struct fp_rep + { + static constexpr bool any_f64 = + (BoundPolicy & bnd::real) == bnd::real || (BoundPolicy & bnd::real) == bnd::real; + static constexpr bool any_f32 = + (BoundPolicy & bnd::f32) == bnd::f32 || (BoundPolicy & bnd::f32) == bnd::f32; + static constexpr bool continuous_ok = AllowContinuous && ResultGrid.Notch == 0; + static constexpr bool keep_f32 = + any_f32 && !any_f64 && (continuous_ok || float_exact); + static constexpr bool keep_f64 = + !keep_f32 && (any_f64 || any_f32) && (continuous_ok || double_exact); + static constexpr bool dropped_fp = (any_f64 || any_f32) && !keep_f64 && !keep_f32; + // Carry both operands' representation flags (widest-wins at storage selection). + static constexpr policy_flag rep = + ((BoundPolicy | BoundPolicy) & (bnd::exact | bnd::direct | bnd::indexed)) + | (keep_f64 ? bnd::real : none) | (keep_f32 ? bnd::f32 : none); + // The result bound's policy: the propagated representation, or plain checked. + static constexpr policy_flag result_policy = rep != none ? rep : checked; + }; +} + //--------------------------------------------------------------------------- // addition — `add(L, R, policy, action) -> bound`, G = Grid + Grid. // The grid arithmetic is sound by construction (the result interval contains @@ -5672,23 +5787,9 @@ namespace bnd::detail "addition: result grid's notch/interval exceeds the representable rational " "range — coarsen the operand grids"); static constexpr grid result_grid = (Grid + Grid).value(); - // Propagate fp storage only when the result grid stays exactly representable - // in the chosen width; otherwise demote (f32→f64) or drop it so storage_pick - // deduces an exact representation (the fp sum would diverge from the exact sum - // — see grid::double_exact / float_exact). Widest-wins: prefer f32 only when - // both operands are f32-only and the result fits float; an f64 operand or a - // too-fine-for-float result widens to f64; too fine for double → exact. - static constexpr bool any_f64 = - (BoundPolicy & bnd::real) == bnd::real || (BoundPolicy & bnd::real) == bnd::real; - static constexpr bool any_f32 = - (BoundPolicy & bnd::f32) == bnd::f32 || (BoundPolicy & bnd::f32) == bnd::f32; - static constexpr bool keep_f32 = any_f32 && !any_f64 && float_exact; - static constexpr bool keep_f64 = !keep_f32 && (any_f64 || any_f32) && double_exact; - // Carry both operands' representation flags (widest-wins at storage selection). - static constexpr policy_flag rep = - ((BoundPolicy | BoundPolicy) & (bnd::exact | bnd::direct | bnd::indexed)) - | (keep_f64 ? bnd::real : none) | (keep_f32 ? bnd::f32 : none); - using result = bound; + // fp / representation propagation — shared rule in detail/rep.hpp. + using rep_t = fp_rep; + using result = bound; template static constexpr bool needs_overflow_check = @@ -5792,22 +5893,11 @@ namespace bnd::detail "multiplication: result grid's notch/interval exceeds the representable " "rational range — coarsen the operand grids"); static constexpr grid result_grid = (Grid * Grid).value(); - // fp storage propagation (see addition.hpp): the product grid (notch = N_L·N_R) - // is finer, so demote f32→f64 / drop f64 when it outgrows the width — the fp - // product would round below the result notch. Widest-wins: f32 only if both - // operands f32-only and product fits float; else widen to f64; else exact. - static constexpr bool any_f64 = - (BoundPolicy & bnd::real) == bnd::real || (BoundPolicy & bnd::real) == bnd::real; - static constexpr bool any_f32 = - (BoundPolicy & bnd::f32) == bnd::f32 || (BoundPolicy & bnd::f32) == bnd::f32; - static constexpr bool keep_f32 = any_f32 && !any_f64 && float_exact; - static constexpr bool keep_f64 = !keep_f32 && (any_f64 || any_f32) && double_exact; - static constexpr bool dropped_fp = (any_f64 || any_f32) && !keep_f64 && !keep_f32; - // Carry both operands' representation flags (widest-wins at storage selection). - static constexpr policy_flag rep = - ((BoundPolicy | BoundPolicy) & (bnd::exact | bnd::direct | bnd::indexed)) - | (keep_f64 ? bnd::real : none) | (keep_f32 ? bnd::f32 : none); - using result = bound; + // fp / representation propagation — shared rule in detail/rep.hpp. The product + // grid (notch = N_L·N_R) is finer, so demotion/dropping is the common case. + using rep_t = fp_rep; + static constexpr bool dropped_fp = rep_t::dropped_fp; + using result = bound; // The dropped-fp case lands on a rational result when the product grid outgrows // uint index space; its product numerator can exceed `umax`, so check it (the @@ -6080,22 +6170,10 @@ namespace bnd::detail ? grid{interval{rational{0}, (Upper / Notch).value()}, Notch} : *(Grid / Grid); - static constexpr bool any_f64 = - (BoundPolicy & bnd::real) == bnd::real || (BoundPolicy & bnd::real) == bnd::real; - static constexpr bool any_f32 = - (BoundPolicy & bnd::f32) == bnd::f32 || (BoundPolicy & bnd::f32) == bnd::f32; - // Keep fp for a continuous result (Notch 0: the raw stores the quotient - // verbatim) or an fp-exact dyadic result; otherwise demote f32→f64 / drop f64 - // (the fp quotient would not land on the result grid). Widest-wins as in mul. - static constexpr bool keep_f32 = - any_f32 && !any_f64 && (result_grid.Notch == 0 || float_exact); - static constexpr bool keep_f64 = - !keep_f32 && (any_f64 || any_f32) && (result_grid.Notch == 0 || double_exact); - // Carry both operands' representation flags (widest-wins at storage selection). - static constexpr policy_flag rep = - ((BoundPolicy | BoundPolicy) & (bnd::exact | bnd::direct | bnd::indexed)) - | (keep_f64 ? bnd::real : none) | (keep_f32 ? bnd::f32 : none); - using result = bound; + // fp / representation propagation — shared rule in detail/rep.hpp. + // AllowContinuous: a continuous quotient (Notch 0) keeps fp verbatim. + using rep_t = fp_rep; + using result = bound; template static constexpr bool needs_overflow_check = @@ -7426,6 +7504,36 @@ namespace bnd boundable auto t) { return a + (b - a) * t; } + //--------------------------------------------------------------------------- + // common_bound — the "hull" type able to hold every value of L and R exactly: + // interval hull + notch gcd (grid `hull`), representation propagated by the + // same widest-wins rule as arithmetic results (detail::fp_rep). Backs the + // std::common_type specialisation (numeric_limits.hpp) and mixed-grid + // min/max below. The primary has no `type` when the hull grid is + // unrepresentable, so common_type_t SFINAEs away instead of erroring. + //--------------------------------------------------------------------------- + namespace detail + { + template + struct common_bound {}; + + // Same type stays itself (policy included) — mirrors std::common_type. + template + struct common_bound { using type = Same; }; + + template + requires (!std::same_as) && (hull(Grid, Grid).has_value()) + struct common_bound + { + static constexpr grid hull_grid = *hull(Grid, Grid); + using type = bound::result_policy>; + }; + } + + template + using common_bound_t = typename detail::common_bound::type; + //--------------------------------------------------------------------------- // std-vocabulary helpers — ADL-found `min` / `max` / `midpoint` for generic // code. min/max mirror std; midpoint returns the *exact* average on a refined @@ -7438,6 +7546,16 @@ namespace bnd template [[nodiscard]] constexpr T max(T a, T b) { return (a < b) ? b : a; } + // Mixed-grid forms return the common hull type (both operands convert + // losslessly — the hull is assignable from each by construction). + template requires (!std::same_as) + [[nodiscard]] constexpr auto min(Lhs a, Rhs b) -> common_bound_t + { common_bound_t ca{a}, cb{b}; return (cb < ca) ? cb : ca; } + + template requires (!std::same_as) + [[nodiscard]] constexpr auto max(Lhs a, Rhs b) -> common_bound_t + { common_bound_t ca{a}, cb{b}; return (ca < cb) ? cb : ca; } + template [[nodiscard]] constexpr auto midpoint(T a, T b) { return (a + b) * just>; } @@ -11247,11 +11365,18 @@ struct std::formatter //--------------------------------------------------------------------------- -// numeric_limits / hash — std:: specialisations for bound. numeric_limits -// reports the *grid* bounds (Lower/Upper), not the raw type's limits. std::hash -// hashes the Raw member (rational raw: Numerator+Denominator, boost-style combine). +// numeric_limits / hash / common_type — std:: specialisations for bound. +// numeric_limits reports the *grid* bounds (Lower/Upper), not the raw type's +// limits. std::hash hashes the Raw member (rational raw: Numerator+Denominator, +// boost-style combine). std::common_type is the grid hull (see bnd::common_bound +// in arithmetic.hpp) so mixed-grid bounds interoperate with generic code. //--------------------------------------------------------------------------- +template +struct std::common_type, bnd::bound> + : bnd::detail::common_bound, bnd::bound> {}; + + template struct std::numeric_limits> { diff --git a/tests/accuracy.cpp b/tests/accuracy.cpp new file mode 100644 index 0000000..75f9e2c --- /dev/null +++ b/tests/accuracy.cpp @@ -0,0 +1,156 @@ +//--------------------------------------------------------------------------- +// accuracy — characterizes each math engine (cordic / dbl / flt) against a +// long-double reference, fpm-style. Reports max/mean error in units +// of the OUTPUT grid notch (an error < 0.5 notch means the engine returned +// the correctly rounded grid point). The long-double reference and +// live only in this test harness — the library itself stays float-free on +// the cordic path. +// +// ./bound_accuracy # markdown table to stdout +// ./bound_accuracy # write the table to +// +// The `accuracy_report` CMake target regenerates docs/accuracy.md with it. +//--------------------------------------------------------------------------- +#include "bound/bound.hpp" +#include "bound/cmath.hpp" + +#include +#include + +using namespace bnd; + +namespace +{ + // Exact value of a bound (or grid constant) as long double, via rational. + long double to_long_double(detail::rational r) + { + return static_cast(r.Numerator) + / static_cast(r.Denominator); + } + + struct stats + { + long double max_err = 0.0L; // in output notches + long double mean_err = 0.0L; + long samples = 0; + }; + + // Sweeps `engine_fn` over the whole input grid (strided to ~`target_samples` + // points) and accumulates the error against `reference_fn` in output-notch + // units. + template + stats sweep(EngineFn engine_fn, ReferenceFn reference_fn, + long target_samples = 20000) + { + const long double lo = to_long_double(Lower); + const long double in_notch = to_long_double(Notch); + const long double span = to_long_double(Upper) - lo; + const long count = static_cast(span / in_notch); + const long stride = count > target_samples ? count / target_samples : 1; + + stats s; + long double sum = 0.0L; + for (long i = 0; i <= count; i += stride) + { + const long double value = lo + static_cast(i) * in_notch; + const In x{static_cast(value)}; + + // Some ops return expected (runtime-failable domains); + // unwrap and skip the rare out-of-domain sample. + const auto outcome = engine_fn(x); + constexpr bool failable = + requires { outcome.has_value(); *outcome; }; + if constexpr (failable) + if (!outcome.has_value()) + continue; + const auto result = [&]{ + if constexpr (failable) return *outcome; else return outcome; }(); + const long double out_notch = + to_long_double(Notch>); + const long double err = + std::fabs(to_long_double(detail::rational{result}) - reference_fn(value)) + / (out_notch > 0.0L ? out_notch : 1.0L); + + if (err > s.max_err) s.max_err = err; + sum += err; + ++s.samples; + } + s.mean_err = s.samples ? sum / static_cast(s.samples) : 0.0L; + return s; + } + + void row(std::FILE* out, const char* function, const char* engine, + const char* domain, stats s) + { + std::fprintf(out, "| %-7s | %-6s | %-14s | %10.4Lf | %10.6Lf | %ld |\n", + function, engine, domain, s.max_err, s.mean_err, s.samples); + } + + // 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.) + 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>; + using sqrt_grid = bound<{{0, 16}, notch<1, 16384>}, round_nearest | real>; + using log_grid = bound<{{1, 1000}, notch<1, 16384>}, round_nearest | real>; +} + +int main(int argc, char** argv) +{ + std::FILE* out = stdout; + if (argc > 1) + { + out = std::fopen(argv[1], "w"); + if (!out) { std::perror(argv[1]); return 1; } + } + + std::fprintf(out, + "# Engine accuracy\n\n" + "Generated by `bound_accuracy` (tests/accuracy.cpp) — do not edit by hand.\n" + "Regenerate with `cmake --build --target accuracy_report`.\n\n" + "Error is measured against a long-double `` reference, in units of\n" + "the **output grid notch**: max < 0.5 means every sampled result is the\n" + "correctly rounded grid point. Engines: `cordic` (integer/CORDIC, FPU-free,\n" + "bit-exact everywhere), `dbl` (reproducible binary64), `flt` (binary32).\n" + "See [math.md](math.md) for the engine model.\n\n" + "| fn | engine | domain | max (notch) | mean (notch) | samples |\n" + "|---------|--------|----------------|-------------|--------------|---------|\n"); + + // One macro row per (function, grid, reference): expands to the three + // engine namespaces so all engines run the identical sweep. +#define BND_ACCURACY(fn, In, domain, ref) \ + row(out, #fn, "cordic", domain, \ + sweep([](In x){ return math::cordic::fn(x); }, \ + [](long double v){ return ref; })); \ + row(out, #fn, "dbl", domain, \ + sweep([](In x){ return math::dbl::fn(x); }, \ + [](long double v){ return ref; })); \ + row(out, #fn, "flt", domain, \ + sweep([](In x){ return math::flt::fn(x); }, \ + [](long double v){ return ref; })); + + BND_ACCURACY(sin, angle_grid, "[-8, 8]", sinl(v)) + BND_ACCURACY(cos, angle_grid, "[-8, 8]", cosl(v)) + BND_ACCURACY(tan, tan_grid, "[-1.5, 1.5]", tanl(v)) + BND_ACCURACY(asin, unit_grid, "[-1, 1]", asinl(v)) + BND_ACCURACY(acos, unit_grid, "[-1, 1]", acosl(v)) + BND_ACCURACY(atan, angle_grid, "[-8, 8]", atanl(v)) + BND_ACCURACY(sinh, angle_grid, "[-8, 8]", sinhl(v)) + BND_ACCURACY(cosh, angle_grid, "[-8, 8]", coshl(v)) + BND_ACCURACY(tanh, angle_grid, "[-8, 8]", tanhl(v)) + BND_ACCURACY(exp, angle_grid, "[-8, 8]", expl(v)) + BND_ACCURACY(exp2, angle_grid, "[-8, 8]", exp2l(v)) + BND_ACCURACY(log, log_grid, "[1, 1000]", logl(v)) + BND_ACCURACY(log2, log_grid, "[1, 1000]", log2l(v)) + BND_ACCURACY(log10, log_grid, "[1, 1000]", log10l(v)) + BND_ACCURACY(sqrt, sqrt_grid, "[0, 16]", sqrtl(v)) + BND_ACCURACY(cbrt, angle_grid, "[-8, 8]", cbrtl(v)) +#undef BND_ACCURACY + + if (out != stdout) std::fclose(out); + return 0; +} diff --git a/tests/fail/CMakeLists.txt b/tests/fail/CMakeLists.txt new file mode 100644 index 0000000..9efc591 --- /dev/null +++ b/tests/fail/CMakeLists.txt @@ -0,0 +1,44 @@ +#--------------------------------------------------------------------------- +# Compile-fail diagnostics suite — each *.cpp here is deliberately ill-formed +# and must FAIL to compile with the message named on its first line +# (`// EXPECT: `; keep it free of regex metacharacters). The test +# invokes the build for the object target; PASS_REGULAR_EXPRESSION passes the +# test only when the build output contains the expected diagnostic — a TU that +# wrongly compiles produces no diagnostic and fails the test. RUN_SERIAL keeps +# concurrent `cmake --build` invocations out of each other's way. +# control_compiles.cpp is the harness control: valid code that must build. +#--------------------------------------------------------------------------- +file(GLOB BOUND_FAIL_SOURCES CONFIGURE_DEPENDS *.cpp) + +foreach(src ${BOUND_FAIL_SOURCES}) + get_filename_component(case_name ${src} NAME_WE) + if(case_name STREQUAL "control_compiles") + continue() + endif() + + file(STRINGS ${src} first_line LIMIT_COUNT 1) + string(REGEX REPLACE "^// EXPECT: *" "" expected_message "${first_line}") + if(NOT first_line MATCHES "^// EXPECT: ") + message(FATAL_ERROR "compile-fail case ${src} lacks a '// EXPECT: ' first line") + endif() + + add_library(fail_${case_name} OBJECT EXCLUDE_FROM_ALL ${src}) + target_link_libraries(fail_${case_name} PRIVATE bound::bound) + + add_test(NAME fail_${case_name} + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} + --target fail_${case_name} --config $) + set_tests_properties(fail_${case_name} PROPERTIES + PASS_REGULAR_EXPRESSION "${expected_message}" + RUN_SERIAL TRUE + LABELS "compile-fail") +endforeach() + +add_library(fail_control OBJECT EXCLUDE_FROM_ALL control_compiles.cpp) +target_link_libraries(fail_control PRIVATE bound::bound) +add_test(NAME fail_control_compiles + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} + --target fail_control --config $) +set_tests_properties(fail_control_compiles PROPERTIES + RUN_SERIAL TRUE + LABELS "compile-fail") diff --git a/tests/fail/control_compiles.cpp b/tests/fail/control_compiles.cpp new file mode 100644 index 0000000..d6ff9fb --- /dev/null +++ b/tests/fail/control_compiles.cpp @@ -0,0 +1,12 @@ +// Control case: valid code that MUST compile through the same harness — a +// broken harness (misconfigured compiler, always-failing builds) fails here. +#include "bound/bound.hpp" + +using namespace bnd; + +int main() +{ + bnd::bound<{0, 100}> percent{50}; + auto sum = percent + 1_b; + (void)sum; +} diff --git a/tests/fail/disjoint_interval_assign.cpp b/tests/fail/disjoint_interval_assign.cpp new file mode 100644 index 0000000..18400ac --- /dev/null +++ b/tests/fail/disjoint_interval_assign.cpp @@ -0,0 +1,11 @@ +// EXPECT: entirely outside lhs interval +// Assigning from a bound whose interval cannot overlap the target is rejected +// at compile time (bound_assignable_why names the failing clause). +#include "bound/bound.hpp" + +int main() +{ + bnd::bound<{0, 10}> small{5}; + bnd::bound<{20, 30}> big{25}; + small = big; // ill-formed: [20,30] never fits in [0,10] +} diff --git a/tests/fail/exclusive_policies.cpp b/tests/fail/exclusive_policies.cpp new file mode 100644 index 0000000..3e1f6ac --- /dev/null +++ b/tests/fail/exclusive_policies.cpp @@ -0,0 +1,10 @@ +// EXPECT: clamp and wrap are mutually exclusive +// Contradictory out-of-range policies on one bound are rejected in the class +// body's mutual-exclusion static_asserts. +#include "bound/bound.hpp" + +int main() +{ + bnd::bound<{0, 100}, bnd::clamp | bnd::wrap> contradictory{}; + (void)contradictory; +} diff --git a/tests/fail/math_without_snap.cpp b/tests/fail/math_without_snap.cpp new file mode 100644 index 0000000..3052260 --- /dev/null +++ b/tests/fail/math_without_snap.cpp @@ -0,0 +1,12 @@ +// EXPECT: must permit rounding +// Transcendentals round their result onto the grid, so the operand must carry +// `snap` (via round_nearest / a round_* mode / real) — require_snap fires. +#include "bound/bound.hpp" +#include "bound/cmath.hpp" + +int main() +{ + bnd::bound<{0, 3}> plain{1}; // no snap permission + auto s = bnd::math::sin(plain); + (void)s; +} diff --git a/tests/fail/modulo_non_integer.cpp b/tests/fail/modulo_non_integer.cpp new file mode 100644 index 0000000..d22c01b --- /dev/null +++ b/tests/fail/modulo_non_integer.cpp @@ -0,0 +1,14 @@ +// EXPECT: modulo requires integer-valued grids and snap +// `%` is integer-only by design — non-integer remainders are not representable +// on a fractional grid without a rounding story. +#include "bound/bound.hpp" + +using namespace bnd; + +int main() +{ + bound<{{0, 8}, notch<1, 2>}> halves{2.5_b}; + bound<{{1, 4}, notch<1, 2>}> divisor{1.5_b}; + auto rem = halves % divisor; // ill-formed: fractional grids + (void)rem; +} diff --git a/tests/fail/notch_mismatch_assign.cpp b/tests/fail/notch_mismatch_assign.cpp new file mode 100644 index 0000000..a71454f --- /dev/null +++ b/tests/fail/notch_mismatch_assign.cpp @@ -0,0 +1,14 @@ +// EXPECT: incompatible notches +// A source lattice that does not land on the target lattice needs explicit +// rounding permission (`with_snap()` / `policy()`); without it the +// assignment is ill-formed. +#include "bound/bound.hpp" + +using namespace bnd; + +int main() +{ + bound<{{0, 1}, notch<1, 2>}> halves{0.5_b}; + bound<{{0, 1}, notch<1, 3>}> thirds{}; + thirds = halves; // ill-formed: 1/2 grid points miss the 1/3 lattice +} diff --git a/tests/fail/scalar_operand.cpp b/tests/fail/scalar_operand.cpp new file mode 100644 index 0000000..23ee5f0 --- /dev/null +++ b/tests/fail/scalar_operand.cpp @@ -0,0 +1,11 @@ +// EXPECT: can only be added to another bound +// Grid-less scalars are rejected with a fix-it (write `a + 1_b`), not a wall +// of overload-resolution noise — the guidance overload's static_assert fires. +#include "bound/bound.hpp" + +int main() +{ + bnd::bound<{0, 100}> percent{50}; + auto sum = percent + 1; // ill-formed: scalar has no grid + (void)sum; +} diff --git a/tests/test_common_type.cpp b/tests/test_common_type.cpp new file mode 100644 index 0000000..5b7b07b --- /dev/null +++ b/tests/test_common_type.cpp @@ -0,0 +1,122 @@ +#include "bound/bound.hpp" +#include "bound/numeric_limits.hpp" + +#include + +#include + +using namespace bnd; + +// Dependent context: a non-dependent invalid requirement inside `requires{}` +// would be a hard error, so probe through a variable template. +template +inline constexpr bool has_common_type = requires { typename std::common_type_t; }; +template +inline constexpr bool has_common_bound = requires { typename common_bound_t; }; + +//--------------------------------------------------------------------------- +// grid hull +//--------------------------------------------------------------------------- +TEST_CASE("hull is the interval hull with the notch gcd", "[grid][hull]") +{ + constexpr grid percent{0, 100}; // notch 1 + constexpr grid halves{interval{-50, 50}, notch<1, 2>}; + + constexpr auto hulled = hull(percent, halves); + STATIC_REQUIRE(hulled.has_value()); + STATIC_REQUIRE(hulled->Interval.Lower == detail::rational{-50}); + STATIC_REQUIRE(hulled->Interval.Upper == detail::rational{100}); + STATIC_REQUIRE(hulled->Notch == (notch<1, 2>)); + + // Coprime denominators combine to the lcm. + constexpr grid half_steps{interval{0, 1}, notch<1, 2>}; + constexpr grid third_steps{interval{0, 1}, notch<1, 3>}; + STATIC_REQUIRE(hull(half_steps, third_steps)->Notch == (notch<1, 6>)); + + // hull(G, G) == G. + STATIC_REQUIRE(hull(percent, percent)->Interval.Lower == percent.Interval.Lower); + STATIC_REQUIRE(hull(percent, percent)->Interval.Upper == percent.Interval.Upper); + STATIC_REQUIRE(hull(percent, percent)->Notch == percent.Notch); + + // A continuous operand makes the hull continuous. + constexpr grid continuous{interval{-1, 1}, detail::rational{0}}; + STATIC_REQUIRE(hull(percent, continuous)->Notch == detail::rational{0}); + STATIC_REQUIRE(hull(percent, continuous)->Interval.Upper == detail::rational{100}); +} + +TEST_CASE("hull result is a valid grid by construction", "[grid][hull]") +{ + // Anchored lattices: both Lowers are notch multiples, so the hull is too. + constexpr grid quarter_steps{interval{-8, 8}, notch<1, 4>}; + constexpr grid unit_steps{interval{-3, 7}, detail::rational{1}}; + constexpr auto hulled = hull(quarter_steps, unit_steps); + STATIC_REQUIRE(hulled.has_value()); + STATIC_REQUIRE(grid::try_make(hulled->Interval, hulled->Notch).has_value()); +} + +//--------------------------------------------------------------------------- +// std::common_type / common_bound_t +//--------------------------------------------------------------------------- +TEST_CASE("common_type of a bound with itself is the bound, policy included", + "[bound][common_type]") +{ + using percent = bound<{0, 100}>; + STATIC_REQUIRE(std::same_as, percent>); + + using clamped_percent = bound<{0, 100}, clamp>; + STATIC_REQUIRE( + std::same_as, clamped_percent>); +} + +TEST_CASE("common_type of mixed grids is the hull type", "[bound][common_type]") +{ + using percent = bound<{0, 100}>; // notch 1 + using halves = bound<{{-50, 50}, notch<1, 2>}>; + using common = std::common_type_t; + + STATIC_REQUIRE(Lower == detail::rational{-50}); + STATIC_REQUIRE(Upper == detail::rational{100}); + STATIC_REQUIRE(Notch == (notch<1, 2>)); + STATIC_REQUIRE(std::same_as>); + STATIC_REQUIRE(std::same_as>); + + // Both operand types convert into the hull losslessly. + constexpr common from_percent{percent{42}}; + constexpr common from_halves{halves{-0.5_b}}; + STATIC_REQUIRE(from_percent == common{42}); + STATIC_REQUIRE(from_halves == common{-0.5_b}); +} + +TEST_CASE("common_type SFINAEs away when the hull notch is unrepresentable", + "[bound][common_type]") +{ + // Coprime ~2^32 denominators: the notch gcd's lcm denominator exceeds imax. + using prime_notch_a = bound<{{0, 1}, notch<1, 4294967291>}>; + using prime_notch_b = bound<{{0, 1}, notch<1, 4294967311>}>; + STATIC_REQUIRE(!has_common_type); + STATIC_REQUIRE(!has_common_bound); +} + +//--------------------------------------------------------------------------- +// mixed-grid min / max +//--------------------------------------------------------------------------- +TEST_CASE("mixed-grid min/max return the hull type", "[bound][min][max]") +{ + using percent = bound<{0, 100}>; + using halves = bound<{{-50, 50}, notch<1, 2>}>; + using common = common_bound_t; + + constexpr percent three{3}; + constexpr halves minus_half{-0.5_b}; + + constexpr auto lo = bnd::min(three, minus_half); + constexpr auto hi = bnd::max(three, minus_half); + STATIC_REQUIRE(std::same_as); + STATIC_REQUIRE(lo == common{-0.5_b}); + STATIC_REQUIRE(hi == common{3}); + + // Same-type overload still returns the operand type. + constexpr auto same = bnd::min(percent{7}, percent{5}); + STATIC_REQUIRE(std::same_as); + STATIC_REQUIRE(same == percent{5}); +} diff --git a/tests/test_rational.cpp b/tests/test_rational.cpp index 6e04f4c..e469d65 100644 --- a/tests/test_rational.cpp +++ b/tests/test_rational.cpp @@ -463,3 +463,25 @@ TEST_CASE("rational trunc / floor / round", "[rational][reduce]") REQUIRE(round((0_r)) == 0); } } + +//--------------------------------------------------------------------------- +// 2026-07: mixed-sign addition rescued in 128-bit. `1024 − m/2^54` forms the +// cross-product 1024·2^54 == 2^64 (one past umax) before the subtraction +// brings the numerator back into range — the dbl-engine store path hit +// exactly this via `rhs - Lower` and terminated through noexcept. +//--------------------------------------------------------------------------- +TEST_CASE("mixed-sign add rescues a cross-product overflow when the difference fits", + "[rational][overflow][regression]") +{ + // -0.49996929771979287 as an exact double fraction: den 2^54. + const rational fp_value{umax{9006646171630191}, imax{-18014398509481984}}; + + const auto offset = fp_value + rational{1024}; + REQUIRE(offset.has_value()); + REQUIRE(offset->Numerator == umax{18437737427537921425u}); + REQUIRE(offset->Denominator == imax{18014398509481984}); + + // Same-sign sums past umax stay nullopt (the result truly needs > 64 bits). + const rational positive{umax{9006646171630191}, imax{18014398509481984}}; + REQUIRE_FALSE((positive + rational{1024}).has_value()); +} diff --git a/tests/test_storage_bugs.cpp b/tests/test_storage_bugs.cpp index ff28047..890d0aa 100644 --- a/tests/test_storage_bugs.cpp +++ b/tests/test_storage_bugs.cpp @@ -162,3 +162,36 @@ TEST_CASE("Bug F: real div-by-zero is reported, not a silent inf", REQUIRE(z.error() == errc::division_by_zero); } #endif // !BND_MATH_FIXED + +//--------------------------------------------------------------------------- +// 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. +//--------------------------------------------------------------------------- +TEST_CASE("fp-derived rational store on a wide snap grid rounds or reports", + "[storage][overflow][regression]") +{ + using wide = bound<{{-1024, 1024}, notch<1, 16384>}, round_nearest>; + + SECTION("negative value: offset fits after the 128-bit rescue") + { + wide slot{}; + slot = rational{umax{9006646171630191}, imax{-18014398509481984}}; // ≈ -0.4999693 + // ·16384 = -8191.49692… → round_nearest → -8191/16384 + REQUIRE(rational{slot} == rational{umax{8191}, imax{-16384}}); + } + + SECTION("positive value: exact offset needs > 64 bits → errc::overflow, no terminate") + { + wide slot{}; + try + { + slot = rational{umax{9006646171630191}, imax{18014398509481984}}; // ≈ +0.4999693 + FAIL("expected the default handler to throw"); + } + catch (bound_error const& e) { REQUIRE(e.code == errc::overflow); } + } +}