From 7ae915c3abefedb2e0f20190c7230eeaa35aa2a5 Mon Sep 17 00:00:00 2001 From: Christian Parpart Date: Thu, 8 Oct 2026 10:00:35 +0200 Subject: [PATCH 1/4] fix(rounding): round at up to 38 decimal places either way Every rounding path refused a place outside -18 to 18, the largest power of ten in a 64-bit integer. Rational has held Int128 since 0.4.0, and from_decimal and _r span 10^-38 to 10^38, so checked_round with SignificantDigits failed for every value below 10^(n - 19): 15 significant digits of anything below 1e-4. MaximumDecimalPlaces, 38, is now the one limit checked_round (both forms), round_wide_ratio, and with it rounded_in_unit and decide_rounding, and rounded_square_root refuse past. An intermediate that does not fit is still Overflow, never a wrong or clamped value: 2/3 at 38 places, an irrational square root past 19 places. checked_decimal_text keeps refusing fewer than -18 places on its own, which it used to inherit from checked_round. rational_from_double's documented limits are recomputed and pinned on both sides: a double below 2^53 rounds at every place up to 22, a fraction over at least 2^N up to 31. Refs #34 Signed-off-by: Christian Parpart --- include/formula-cpp/detail/wide_rounding.hpp | 20 +++-- include/formula-cpp/number_text.hpp | 12 +-- include/formula-cpp/rounded_root.hpp | 10 ++- include/formula-cpp/rounding.hpp | 37 ++++++--- include/formula-cpp/unit.hpp | 7 +- test/measured_tests.cpp | 8 +- test/number_text_tests.cpp | 6 +- test/rounded_root_tests.cpp | 30 ++++++- test/rounding_tests.cpp | 85 +++++++++++++++++--- test/unit_tests.cpp | 9 ++- test/wide_rounding_tests.cpp | 38 +++++++-- 11 files changed, 201 insertions(+), 61 deletions(-) diff --git a/include/formula-cpp/detail/wide_rounding.hpp b/include/formula-cpp/detail/wide_rounding.hpp index 9ed8432..4e01dc4 100644 --- a/include/formula-cpp/detail/wide_rounding.hpp +++ b/include/formula-cpp/detail/wide_rounding.hpp @@ -8,11 +8,13 @@ /// /// `round_wide_ratio` is built to agree with `checked_round`: for a value /// `checked_round` accepts, it gives the same result in all seven modes at every -/// place from -18 to 18 (`test/wide_rounding_tests.cpp` checks 43729 such -/// cases). It also answers some values `checked_round` refuses -- where a -/// numerator times 10^places leaves 128 bits but the rounded result fits -- and -/// refuses, with `Overflow`, a result that does not fit `Rational`, as -/// `checked_round` does. +/// place from -38 to 38 (`MaximumDecimalPlaces`), in a width that holds any +/// `Rational::Int` times 10^38 -- 256 bits, `WideRatio<8>` +/// (`test/wide_rounding_tests.cpp` checks 95498 such cases). A narrower width +/// refuses with `Overflow` where the scaled value leaves it. It also answers +/// some values `checked_round` refuses -- where a numerator times 10^places +/// leaves 128 bits but the rounded result fits -- and refuses, with +/// `Overflow`, a result that does not fit `Rational`, as `checked_round` does. /// /// Every fraction is reduced before it is scaled, so an operation may hand over /// one far from lowest terms: the width only has to hold the value, not the way @@ -80,8 +82,10 @@ template /// The decimal @p unrounded rounds to at @p places under @p roundingMode, as /// `Rational::from_decimal(kept, -places)`. `Overflow` for places outside -/// -18 to 18, or a kept integer that does not fit `Rational::Int`, as -/// `checked_round` reports them; `DivisionByZero` for a zero denominator. +/// -`MaximumDecimalPlaces` to `MaximumDecimalPlaces`, or a kept integer that +/// does not fit `Rational::Int`, as `checked_round` reports them, and where +/// the value scaled by 10^places leaves `L` limbs; `DivisionByZero` for a +/// zero denominator. /// /// The magnitude is rounded and the sign applied afterwards: the directed /// modes read the sign (Floor rounds a negative magnitude up, Ceiling a @@ -92,7 +96,7 @@ template DecimalPlaces places, RoundingMode roundingMode) noexcept { - if (places.value > 18 || places.value < -18) + if (places.value > MaximumDecimalPlaces || places.value < -MaximumDecimalPlaces) return std::unexpected { ArithmeticError::Overflow }; if (unrounded.denominator.is_zero()) return std::unexpected { ArithmeticError::DivisionByZero }; diff --git a/include/formula-cpp/number_text.hpp b/include/formula-cpp/number_text.hpp index cd12adc..6c96cc7 100644 --- a/include/formula-cpp/number_text.hpp +++ b/include/formula-cpp/number_text.hpp @@ -427,15 +427,15 @@ namespace detail /// `checked_round`'s arithmetic, so it can overflow where the rounding does, /// and it is the only one here that reports to the overflow census. /// -/// @return `Overflow` for more than 18 places -- `checked_round`'s own -/// refusal, since 10^-19 is not a representable step -- or wherever -/// `checked_round` overflows at a negative @p places; `DomainError` -/// for a tie under a rounding mode that is none of the seven. -/// Nothing saturates. +/// @return `Overflow` for places outside -18 to 18 (`ExactDecimalPlaces`) +/// -- the text layer's own limit, though `checked_round` reaches 38 +/// -- or wherever `checked_round` overflows at a negative @p places; +/// `DomainError` for a tie under a rounding mode that is none of the +/// seven. Nothing saturates. [[nodiscard]] constexpr std::expected checked_decimal_text( Rational unrounded, DecimalPlaces places, RoundingMode roundingMode, DecimalPadding decimalPadding) noexcept { - if (places.value > detail::ExactDecimalPlaces) + if (places.value > detail::ExactDecimalPlaces || places.value < -detail::ExactDecimalPlaces) return std::unexpected { ArithmeticError::Overflow }; if (places.value < 0) diff --git a/include/formula-cpp/rounded_root.hpp b/include/formula-cpp/rounded_root.hpp index 62d98b1..3ccdeb4 100644 --- a/include/formula-cpp/rounded_root.hpp +++ b/include/formula-cpp/rounded_root.hpp @@ -136,8 +136,12 @@ namespace detail /// does for `checked_round`: then S is 1/10^(2|p|), and it is B that grows /// instead of q. /// - /// Every intermediate is a 128-bit unsigned integer (`detail::UInt128`). - /// So the headroom is `floor(v) * 10^(2p) < 2^128` (and + /// Every intermediate is a 128-bit unsigned integer (`detail::UInt128`), + /// and so is S = 10^(2|p|) itself: an irrational root rounds at no more + /// than 19 places either way, whatever the radicand -- 10^38 fits 128 bits + /// and 10^40 does not -- though `MaximumDecimalPlaces` allows 38. A + /// rational root is `checked_round`'s, and takes every place it does. + /// Within those 19 the headroom is `floor(v) * 10^(2p) < 2^128` (and /// `b * 10^(2p) < 2^128` for the remainder): an integer radicand of about /// 10^6 fits at 16 places and overflows at 17. **The bound is on the /// denominator b too**, whatever the value: at p places a denominator above @@ -155,7 +159,7 @@ namespace detail if (radicandInUnitSquared.sign() < 0) return std::unexpected { ArithmeticError::DomainError }; // As `checked_round` reports them, whichever path the root takes. - if (places.value > 18 || places.value < -18) + if (places.value > MaximumDecimalPlaces || places.value < -MaximumDecimalPlaces) return std::unexpected { ArithmeticError::Overflow }; std::expected const exactRoot = checked_exact_nth_root(radicandInUnitSquared, 2); diff --git a/include/formula-cpp/rounding.hpp b/include/formula-cpp/rounding.hpp index 76f201e..406fb2d 100644 --- a/include/formula-cpp/rounding.hpp +++ b/include/formula-cpp/rounding.hpp @@ -88,6 +88,16 @@ struct DecimalPlaces [[nodiscard]] constexpr bool operator==(DecimalPlaces const&) const noexcept = default; }; +/// The most decimal places rounding takes either way: a place outside +/// `-MaximumDecimalPlaces` to `MaximumDecimalPlaces` is `ArithmeticError::Overflow`. +/// +/// The scale `Rational::from_decimal` spans, 10^-38 to 10^38: the largest power +/// of ten `Rational::Int` holds. A place inside the range can still overflow +/// where an intermediate does not fit, and that is reported the same way, never +/// a wrong or clamped value. Text output stops at 18 places of its own +/// (`number_text.hpp`). +inline constexpr std::int32_t MaximumDecimalPlaces = 38; + /// A count of significant digits. Must be at least 1. struct SignificantDigits { @@ -258,12 +268,16 @@ namespace detail } /// Rounds `unrounded` to `places` decimal places under `roundingMode`. +/// +/// `Overflow` for places outside -`MaximumDecimalPlaces` to +/// `MaximumDecimalPlaces`, or where `unrounded` scaled by 10^places does not +/// fit `Rational::Int`. [[nodiscard]] constexpr std::expected checked_round(Rational unrounded, DecimalPlaces places, RoundingMode roundingMode) noexcept { // The step is 10^-places, which must itself be representable. - if (places.value > 18 || places.value < -18) + if (places.value > MaximumDecimalPlaces || places.value < -MaximumDecimalPlaces) return std::unexpected { ArithmeticError::Overflow }; std::expected const increment = Rational::from_decimal(1, -places.value); @@ -295,7 +309,7 @@ namespace detail // Keeping `significant` digits of a value whose leading digit sits at 10^e // means rounding at the 10^(e - significant + 1) place. long long const places = static_cast(significant.value) - 1 - static_cast(*exponent); - if (places > 18 || places < -18) + if (places > MaximumDecimalPlaces || places < -MaximumDecimalPlaces) return std::unexpected { ArithmeticError::Overflow }; return checked_round(unrounded, DecimalPlaces { static_cast(places) }, roundingMode); @@ -325,15 +339,20 @@ namespace detail /// `Rational::Int` is `|numerator| * (10^N / gcd(10^N, denominator))` -- for /// a binary denominator, `|numerator| * 5^N`. There the limit is set by the /// **numerator's** magnitude, not the denominator's and not the value's -/// size: `1 / 2^121` rounds at all 18 places, while a 100-bit numerator over -/// the same denominator does not. A `double` below 2^53 in magnitude has a -/// numerator of at most 53 bits, and rounding it forms at most -/// 2^53 * 5^18 * 2^18, below 2^113, so it rounds at every place from 0 to -/// 18. A whole `double` past that has a numerator of its own magnitude and -/// nothing to cancel: 1e21 is refused at 18 places, and 1e38 at 1. +/// size: `1 / 2^121` rounds at all 38 places, while a 100-bit numerator over +/// the same denominator rounds at 11 and is refused at 12. A `double` below +/// 2^53 in magnitude has a numerator of at most 53 bits, and rounding it +/// forms at most 2^53 * 10^N, below 2^127 up to N = 22, so it rounds at +/// every place from 0 to 22; 2^53 - 1 is refused at 23. A fraction whose +/// denominator is at least 2^N cancels all of the 2^N and forms at most +/// 2^53 * 5^N, which fits up to N = 31: 0.45, `8106479329266893 / 2^54`, +/// rounds at 31 places and is refused at 32. A whole `double` past 2^53 has +/// a numerator of its own magnitude and nothing to cancel: 1e21 is refused +/// at 18 places, and 1e38 at 1. /// - Rounding to a NEGATIVE number of places -- to whole tens or hundreds -- /// uses an integer step, which multiplies the **denominator** instead. There -/// the denominator is the constraint: `2^-100` is refused at -18 places. +/// the denominator is the constraint: `2^-100` rounds at -8 places and is +/// refused at -9, since 2^100 * 10^9 is past 2^127. /// /// Prefer `from_decimal` for an exact decimal; use this function only for a /// genuinely measured `double`. diff --git a/include/formula-cpp/unit.hpp b/include/formula-cpp/unit.hpp index 859c034..e218ef8 100644 --- a/include/formula-cpp/unit.hpp +++ b/include/formula-cpp/unit.hpp @@ -903,9 +903,10 @@ inline constexpr bool formats_by_describe = true; } /// A rounding to the decimal places @p roundedIn declares, under @p roundingMode. -/// A unit whose declared decimals lie outside -18 to 18 is not refused here: -/// `checked_round` returns `ArithmeticError::Overflow` for that many places, so -/// evaluating the rounding does. +/// A unit whose declared decimals lie outside -38 to 38 +/// (`MaximumDecimalPlaces`) is not refused here: `checked_round` returns +/// `ArithmeticError::Overflow` for that many places, so evaluating the +/// rounding does. [[nodiscard]] constexpr DecimalRounding declared_rounding(Unit roundedIn, RoundingMode roundingMode) noexcept { return DecimalRounding { roundedIn, declared_decimals(roundedIn), roundingMode }; diff --git a/test/measured_tests.cpp b/test/measured_tests.cpp index 3cb91aa..9a400c6 100644 --- a/test/measured_tests.cpp +++ b/test/measured_tests.cpp @@ -343,12 +343,12 @@ TEST_CASE("a conversion overflow is reported as an error, not re-read as absence namespace { -// decimals = 19: the same unrepresentable precision unit_tests.cpp uses to -// exercise checked_round_to_declared's own error path, carried through a -// quantity this time. +// decimals = 39, one past `MaximumDecimalPlaces`: the same unrepresentable +// precision unit_tests.cpp uses to exercise checked_round_to_declared's own +// error path, carried through a quantity this time. inline constexpr formula::Unit UnrepresentablePrecision { .dimension = formula::dim::Scalar, .symbolText = formula::symbol("bad"), - .decimals = 19 }; + .decimals = 39 }; struct UnroundableReading: formula::Quantityis_exact()); } -TEST_CASE("more places than DecimalPlaces spans is refused as checked_round refuses it", "[number_text]") +TEST_CASE("more places than text is written at, 18 either way, is refused", "[number_text]") { + // Text's own limit, not checked_round's, which rounds 1/3 at 19 places and at -19. + REQUIRE(formula::checked_round(Rational { 1, 3 }, DecimalPlaces { 19 }, RoundingMode::HalfEven).has_value()); + REQUIRE(formula::checked_round(Rational { 1, 3 }, DecimalPlaces { -19 }, RoundingMode::HalfEven).has_value()); + auto const tooFine = formula::checked_decimal_text( Rational { 1, 3 }, DecimalPlaces { 19 }, RoundingMode::HalfEven, DecimalPadding::Padded); REQUIRE_FALSE(tooFine.has_value()); diff --git a/test/rounded_root_tests.cpp b/test/rounded_root_tests.cpp index 9854823..b43c9dd 100644 --- a/test/rounded_root_tests.cpp +++ b/test/rounded_root_tests.cpp @@ -14,6 +14,7 @@ using formula::DecimalPlaces; using formula::Rational; using formula::RoundingMode; using formula::var; +using namespace formula::literals; /// A gram squared, the unit a variance of masses in grams is stated in. /// Invented here rather than shipped: the library has no squared mass unit. @@ -182,15 +183,36 @@ TEST_CASE("rounded_sqrt reports overflow rather than a wrapped result", "[rounde STATIC_REQUIRE(outcome.error() == formula::ArithmeticError::Overflow); } -TEST_CASE("rounded_sqrt takes places from -18 to 18, as checked_round does", "[rounded_root]") +TEST_CASE("rounded_sqrt takes places from -38 to 38, as checked_round does", "[rounded_root]") { using formula::detail::rounded_square_root; - // The root of 2 is irrational, so this is the digit-by-digit path, not checked_round's. + // The root of 2 is irrational, so this is the digit-by-digit path, not checked_round's: 18 places, and + // 19, which 2 * 10^38 still holds. At 20 places 10^40 leaves 128 bits, for every irrational root. STATIC_REQUIRE(rounded_square_root(Rational { 2 }, DecimalPlaces { 18 }, RoundingMode::Floor).value() == *Rational::from_decimal(1'414'213'562'373'095'048, -18)); - STATIC_REQUIRE(rounded_square_root(Rational { 2 }, DecimalPlaces { 19 }, RoundingMode::Floor).error() + STATIC_REQUIRE(rounded_square_root(Rational { 2 }, DecimalPlaces { 19 }, RoundingMode::Floor).value() + == 1.4142135623730950488_r); + STATIC_REQUIRE(rounded_square_root(Rational { 2 }, DecimalPlaces { 19 }, RoundingMode::Ceiling).value() + == 1.4142135623730950489_r); + STATIC_REQUIRE(rounded_square_root(Rational { 2 }, DecimalPlaces { 20 }, RoundingMode::Floor).error() == formula::ArithmeticError::Overflow); - STATIC_REQUIRE(rounded_square_root(Rational { 2 }, DecimalPlaces { -19 }, RoundingMode::Floor).error() + // At -19 places the divisor is 10^38, which fits; at -20 it does not. + STATIC_REQUIRE(rounded_square_root(Rational { 2 }, DecimalPlaces { -19 }, RoundingMode::Ceiling).value() == 1e19_r); + STATIC_REQUIRE(rounded_square_root(Rational { 2 }, DecimalPlaces { -20 }, RoundingMode::Ceiling).error() + == formula::ArithmeticError::Overflow); + // A rational root is checked_round's, and takes every place it does: the root of 1/9 at 38 places, and + // of 4 at -38. + STATIC_REQUIRE(rounded_square_root(Rational { 1, 9 }, DecimalPlaces { 38 }, RoundingMode::Floor).value() + == 0.33333333333333333333333333333333333333_r); + STATIC_REQUIRE(rounded_square_root(Rational { 4 }, DecimalPlaces { -38 }, RoundingMode::Ceiling).value() == 1e38_r); + // One past either way is refused, whichever path the root takes. + STATIC_REQUIRE(rounded_square_root(Rational { 2 }, DecimalPlaces { 39 }, RoundingMode::Floor).error() + == formula::ArithmeticError::Overflow); + STATIC_REQUIRE(rounded_square_root(Rational { 2 }, DecimalPlaces { -39 }, RoundingMode::Floor).error() + == formula::ArithmeticError::Overflow); + STATIC_REQUIRE(rounded_square_root(Rational { 1, 9 }, DecimalPlaces { 39 }, RoundingMode::Floor).error() + == formula::ArithmeticError::Overflow); + STATIC_REQUIRE(rounded_square_root(Rational { 1, 9 }, DecimalPlaces { -39 }, RoundingMode::Floor).error() == formula::ArithmeticError::Overflow); } diff --git a/test/rounding_tests.cpp b/test/rounding_tests.cpp index 8e1dd40..bf03c27 100644 --- a/test/rounding_tests.cpp +++ b/test/rounding_tests.cpp @@ -10,6 +10,7 @@ using formula::ArithmeticError; using formula::ArithmeticException; using formula::Rational; using formula::RoundingMode; +using namespace formula::literals; namespace { @@ -208,7 +209,45 @@ static_assert(formula::checked_round(Rational { 5 }, formula::SignificantDigits == ArithmeticError::DomainError); static_assert(formula::checked_round(Rational { 5 }, formula::SignificantDigits { -1 }, RoundingMode::Floor).error() == ArithmeticError::DomainError); -static_assert(formula::checked_round(Rational { 5 }, formula::DecimalPlaces { 19 }, RoundingMode::Floor).error() + +// ---- the places rounding reaches: MaximumDecimalPlaces either way, the scale from_decimal spans ---- + +static_assert(formula::MaximumDecimalPlaces == 38); +// 19 places either way, one past the limit while Rational was 64-bit. +static_assert(formula::checked_round(Rational { 1, 3 }, formula::DecimalPlaces { 19 }, RoundingMode::Floor) + == 0.3333333333333333333_r); +static_assert(formula::checked_round(5e19_r + 1, formula::DecimalPlaces { -19 }, RoundingMode::HalfEven) == 5e19_r); +// 38 places, the last: 1/3 to 10^-38, and 1/6 rounded up in that place. +static_assert(formula::checked_round(Rational { 1, 3 }, formula::DecimalPlaces { 38 }, RoundingMode::Floor) + == 0.33333333333333333333333333333333333333_r); +static_assert(formula::checked_round(Rational { 1, 6 }, formula::DecimalPlaces { 38 }, RoundingMode::HalfEven) + == 0.16666666666666666666666666666666666667_r); +// A place inside the range still overflows where the scaled value does not fit: 2/3 * 10^38 is past 2^127, +// though the rounding, 0.66666666666666666666666666666666666667, would fit. +static_assert(formula::checked_round(Rational { 2, 3 }, formula::DecimalPlaces { 38 }, RoundingMode::HalfEven).error() + == ArithmeticError::Overflow); +// -38 places, to whole 10^38s: 1.4 * 10^38 rounds down to 10^38, and 7 to 0 or up to 10^38. +static_assert(formula::checked_round(14e37_r, formula::DecimalPlaces { -38 }, RoundingMode::HalfEven) == 1e38_r); +static_assert(formula::checked_round(Rational { 7 }, formula::DecimalPlaces { -38 }, RoundingMode::Ceiling) == 1e38_r); +static_assert(formula::checked_round(Rational { 7 }, formula::DecimalPlaces { -38 }, RoundingMode::Floor) == Rational {}); +// One past, either way, on zero, which rounds to itself at any place that is accepted. +static_assert(formula::checked_round(Rational {}, formula::DecimalPlaces { 39 }, RoundingMode::Floor).error() + == ArithmeticError::Overflow); +static_assert(formula::checked_round(Rational {}, formula::DecimalPlaces { -39 }, RoundingMode::Floor).error() + == ArithmeticError::Overflow); +static_assert(formula::checked_round(Rational { 5 }, formula::DecimalPlaces { 39 }, RoundingMode::Floor).error() + == ArithmeticError::Overflow); + +// Significant digits reach the same places. 15 digits of -2.12345678901234e-5 are -212345678901234 / 10^19, +// at 19 places, and of -1/30000 -0.0000333333333333333. +static_assert(formula::checked_round(-2.12345678901234e-5_r, formula::SignificantDigits { 15 }, RoundingMode::HalfEven) + == -0.0000212345678901234_r); +static_assert(formula::checked_round(Rational { -1, 30'000 }, formula::SignificantDigits { 15 }, RoundingMode::HalfEven) + == -0.0000333333333333333_r); +// One digit of 3 * 10^-38 is at place 38; one digit of 1 / (1.5 * 10^38), 6.7 * 10^-39, would be at place 39. +static_assert(formula::checked_round(3e-38_r, formula::SignificantDigits { 1 }, RoundingMode::HalfEven) == 3e-38_r); +static_assert(formula::checked_round(Rational { 1 } / 15e37_r, formula::SignificantDigits { 1 }, RoundingMode::Floor) + .error() == ArithmeticError::Overflow); TEST_CASE("intermediate and final rounding compose, and the order matters", "[rounding]") @@ -295,30 +334,46 @@ TEST_CASE("rounding precision is limited by the numerator, not the magnitude", " // |numerator| * (10^N / gcd(10^N, denominator)), // which for a binary denominator is |numerator| * 5^N. The NUMERATOR // decides. A double below 2^53 in magnitude has one of at most 53 bits, - // and rounding it at up to 18 places forms at most 2^53 * 5^18 * 2^18, - // below 2^113, so every place from 0 to 18 fits. + // and rounding it at N places forms at most 2^53 * 10^N, below 2^127 up + // to N = 22, so every place from 0 to 22 fits. A fraction whose + // denominator is at least 2^N forms at most 2^53 * 5^N, which fits up to + // N = 31. // 0.45 as a double is exactly 8106479329266893 / 2^54 -- a 53-bit - // numerator. Six places, which overflowed 64 bits, and eighteen fit. + // numerator. Six places, which overflowed 64 bits, eighteen, and 31 fit; + // 32 does not: 8106479329266893 * 5^32 is past 2^127. CHECK(formula::rational_from_double(0.45, formula::DecimalPlaces { 2 }, RoundingMode::HalfAwayFromZero).has_value()); CHECK(formula::rational_from_double(0.45, formula::DecimalPlaces { 6 }, RoundingMode::HalfAwayFromZero) == Rational { 9, 20 }); CHECK(formula::rational_from_double(0.45, formula::DecimalPlaces { 18 }, RoundingMode::HalfAwayFromZero) == Rational { 450'000'000'000'000'011, 1'000'000'000'000'000'000 }); + CHECK(formula::rational_from_double(0.45, formula::DecimalPlaces { 31 }, RoundingMode::HalfAwayFromZero) + == 0.4500000000000000111022302462516_r); + CHECK(error_of(formula::rational_from_double(0.45, formula::DecimalPlaces { 32 }, RoundingMode::HalfAwayFromZero)) + == ArithmeticError::Overflow); // Just below 2^52, half a unit off a whole number, it is exact at 18 places. CHECK(formula::rational_from_double(4503599627370495.5, formula::DecimalPlaces { 18 }, RoundingMode::Floor) == Rational { 9007199254740991, 2 }); - // And 2^53 - 1, the largest whole number below 2^53, is exact at 18. + // And 2^53 - 1, the largest whole number below 2^53, is exact at 18 and + // at 22, where it cancels nothing; at 23, 2^53 * 10^23 is past 2^127. CHECK(formula::rational_from_double(9007199254740991.0, formula::DecimalPlaces { 18 }, RoundingMode::Floor) == Rational { 9007199254740991 }); + CHECK(formula::rational_from_double(9007199254740991.0, formula::DecimalPlaces { 22 }, RoundingMode::Floor) + == Rational { 9007199254740991 }); + CHECK(error_of(formula::rational_from_double(9007199254740991.0, formula::DecimalPlaces { 23 }, RoundingMode::Floor)) + == ArithmeticError::Overflow); // Past 2^53 the claim stops: a whole double's numerator is its own - // magnitude, with nothing to cancel, so 1e21 overflows at 18 places and - // 1e38 at 1; and at -18 places 2^-100's denominator times 10^18 does. + // magnitude, with nothing to cancel, so 1e21 rounds at 17 places and + // overflows at 18, and 1e38 overflows at 1. + CHECK(formula::rational_from_double(1e21, formula::DecimalPlaces { 17 }, RoundingMode::Floor) == 1e21_r); CHECK(error_of(formula::rational_from_double(1e21, formula::DecimalPlaces { 18 }, RoundingMode::Floor)) == ArithmeticError::Overflow); CHECK(error_of(formula::rational_from_double(1e38, formula::DecimalPlaces { 1 }, RoundingMode::Floor)) == ArithmeticError::Overflow); - CHECK(error_of(formula::rational_from_double(0x1p-100, formula::DecimalPlaces { -18 }, RoundingMode::Floor)) + // At a negative number of places the denominator is scaled: 2^-100 rounds + // up to 10^8 at -8 places, and at -9 2^100 * 10^9 is past 2^127. + CHECK(formula::rational_from_double(0x1p-100, formula::DecimalPlaces { -8 }, RoundingMode::Ceiling) == 1e8_r); + CHECK(error_of(formula::rational_from_double(0x1p-100, formula::DecimalPlaces { -9 }, RoundingMode::Ceiling)) == ArithmeticError::Overflow); // What does fail is from_double_exact, before any rounding happens, when @@ -339,12 +394,20 @@ TEST_CASE("rounding precision is limited by the numerator, not the magnitude", " .has_value()); // The denominator is not what limits it: numerator 1 over 2^121 rounds at - // every supported place, while a 100-bit numerator over the SAME - // denominator does not. This pair is what distinguishes the two - // explanations. + // every supported place, up to 38, while a 100-bit numerator over the SAME + // denominator rounds at 11 places and not at 12. This pair is what + // distinguishes the two explanations. constexpr Rational::Int twoTo121 = Rational::Int { 1 } << 121; CHECK(formula::checked_round(Rational { 1, twoTo121 }, formula::DecimalPlaces { 18 }, RoundingMode::Floor) .has_value()); + CHECK(formula::checked_round(Rational { 1, twoTo121 }, formula::DecimalPlaces { 38 }, RoundingMode::Ceiling) + == 38e-38_r); // 1/2^121 is 37.6 units of 10^-38 + CHECK(formula::checked_round( + Rational { (Rational::Int { 1 } << 100) - 1, twoTo121 }, formula::DecimalPlaces { 11 }, RoundingMode::Floor) + == 0.00000047683_r); + CHECK(error_of(formula::checked_round( + Rational { (Rational::Int { 1 } << 100) - 1, twoTo121 }, formula::DecimalPlaces { 12 }, RoundingMode::Floor)) + == ArithmeticError::Overflow); CHECK(error_of(formula::checked_round( Rational { (Rational::Int { 1 } << 100) - 1, twoTo121 }, formula::DecimalPlaces { 18 }, RoundingMode::Floor)) == ArithmeticError::Overflow); diff --git a/test/unit_tests.cpp b/test/unit_tests.cpp index ecee87f..d78b472 100644 --- a/test/unit_tests.cpp +++ b/test/unit_tests.cpp @@ -968,10 +968,11 @@ TEST_CASE("rounding to a unit's declared precision uses that unit's decimals", " TEST_CASE("round_to_declared throws exactly where checked_round_to_declared reports an error", "[unit]") { // A unit whose declared decimals fall outside what a rounding step can - // represent -- the same guard checked_round's own tests exercise directly. - // Pinned the way convert()/checked_convert() are: the checked form reports - // the error, and the throwing form throws it, on the very same input. - Unit const unrepresentableDecimals { .dimension = dim::Scalar, .symbolText = formula::symbol("bad"), .decimals = 19 }; + // represent, 39, one past `MaximumDecimalPlaces` -- the same guard + // checked_round's own tests exercise directly. Pinned the way + // convert()/checked_convert() are: the checked form reports the error, + // and the throwing form throws it, on the very same input. + Unit const unrepresentableDecimals { .dimension = dim::Scalar, .symbolText = formula::symbol("bad"), .decimals = 39 }; auto const checked = formula::checked_round_to_declared( *Rational::make(1, 1), unrepresentableDecimals, formula::RoundingMode::HalfAwayFromZero); diff --git a/test/wide_rounding_tests.cpp b/test/wide_rounding_tests.cpp index 123ed2f..c89185e 100644 --- a/test/wide_rounding_tests.cpp +++ b/test/wide_rounding_tests.cpp @@ -26,6 +26,7 @@ using formula::detail::narrow_wide_ratio; using formula::detail::round_wide_ratio; using formula::detail::rounded_in_unit; using formula::detail::wide_from_rational; +using namespace formula::literals; constexpr std::array everyMode { RoundingMode::HalfAwayFromZero, RoundingMode::HalfTowardZero, RoundingMode::HalfEven, RoundingMode::Ceiling, RoundingMode::Floor, RoundingMode::TowardZero, @@ -55,13 +56,14 @@ constexpr W4 pastIntMin = W4::from_u128(UInt128 { std::uint64_t { 1 } << 63, 1 } TEST_CASE("wide rounding: a wide fraction rounds as checked_round rounds the same Rational", "[wide-rounding]") { // Wherever checked_round answers, round_wide_ratio gives the same, in - // every mode and at every place it accepts. Counted: 48833 cases, those - // of the grid and the four extremes, at every place from -18 to 18 and in - // every mode, where checked_round answers. + // every mode and at every place it accepts. Counted: 95498 cases, those + // of the grid and the four extremes, at every place from -38 to 38 and in + // every mode, where checked_round answers. In 256 bits, which hold any + // `Rational::Int` times 10^38. int checkedCases = 0; int agreedCases = 0; auto const judge = [&](Rational exact) { - for (int places = -18; places <= 18; ++places) + for (int places = -formula::MaximumDecimalPlaces; places <= formula::MaximumDecimalPlaces; ++places) for (RoundingMode const each: everyMode) { std::expected const expected = @@ -70,7 +72,7 @@ TEST_CASE("wide rounding: a wide fraction rounds as checked_round rounds the sam continue; ++checkedCases; std::expected const wide = - round_wide_ratio(from(exact), DecimalPlaces { places }, each); + round_wide_ratio(wide_from_rational<8>(exact), DecimalPlaces { places }, each); if (wide.has_value() && *wide == *expected) ++agreedCases; } @@ -81,10 +83,25 @@ TEST_CASE("wide rounding: a wide fraction rounds as checked_round rounds the sam for (Rational const extreme: { Rational { IntMin }, Rational { IntMax }, Rational { IntMin, 3 }, Rational { IntMax, 2 } }) judge(extreme); - CHECK(checkedCases == 48833); + CHECK(checkedCases == 95498); CHECK(agreedCases == checkedCases); } +TEST_CASE("wide rounding: a width narrower than 256 bits refuses what leaves it", "[wide-rounding]") +{ + // 7/2^100 at 38 places: checked_round cancels 2^38 and forms 7 * 5^38 / 2^62, which fits; the wide + // fraction scales 7 by 10^38 first, past 2^128, so 128 bits refuse it and 256 bits answer as + // checked_round does. + constexpr Rational sevenOverTwoTo100 { 7, Rational::Int { 1 } << 100 }; + constexpr std::expected rounded = + formula::checked_round(sevenOverTwoTo100, DecimalPlaces { 38 }, RoundingMode::HalfEven); + STATIC_REQUIRE(rounded.has_value()); + STATIC_REQUIRE(round_wide_ratio(from(sevenOverTwoTo100), DecimalPlaces { 38 }, RoundingMode::HalfEven).error() + == formula::ArithmeticError::Overflow); + STATIC_REQUIRE(round_wide_ratio(wide_from_rational<8>(sevenOverTwoTo100), DecimalPlaces { 38 }, RoundingMode::HalfEven) + == rounded); +} + TEST_CASE("wide rounding: a value checked_round refuses is answered when its rounding fits", "[wide-rounding]") { // (10^21 + 1)/9 at 18 dp: checked_round forms (10^21 + 1) * 10^18, past @@ -203,9 +220,14 @@ TEST_CASE("wide rounding: the kept integer must fit Rational and the places must DecimalPlaces { -18 }, RoundingMode::HalfEven) == Rational { Rational::Int { 1181 } * 1'000'000'000'000'000'000 }); - STATIC_REQUIRE(round_wide_ratio(from(Rational { 1, 3 }), DecimalPlaces { 19 }, RoundingMode::HalfEven).error() + // Places from -38 to 38, as checked_round takes them: 1/3 at 38 places and 7 at -38, and one past + // either way refused -- on zero in 256 bits, where 10^39 fits and zero would otherwise round to itself. + STATIC_REQUIRE(*round_wide_ratio(from(Rational { 1, 3 }), DecimalPlaces { 38 }, RoundingMode::HalfEven) + == 0.33333333333333333333333333333333333333_r); + STATIC_REQUIRE(*round_wide_ratio(from(Rational { 7 }), DecimalPlaces { -38 }, RoundingMode::Ceiling) == 1e38_r); + STATIC_REQUIRE(round_wide_ratio(wide_from_rational<8>(Rational {}), DecimalPlaces { 39 }, RoundingMode::HalfEven).error() == formula::ArithmeticError::Overflow); - STATIC_REQUIRE(round_wide_ratio(from(Rational { 1, 3 }), DecimalPlaces { -19 }, RoundingMode::HalfEven).error() + STATIC_REQUIRE(round_wide_ratio(wide_from_rational<8>(Rational {}), DecimalPlaces { -39 }, RoundingMode::HalfEven).error() == formula::ArithmeticError::Overflow); STATIC_REQUIRE( round_wide_ratio(R4 { false, W4::from_u64(1), W4 {} }, DecimalPlaces { 0 }, RoundingMode::HalfEven).error() From 061436b655eb16de92b2e0d85e91a79ed5a780b1 Mon Sep 17 00:00:00 2001 From: Christian Parpart Date: Thu, 8 Oct 2026 10:01:57 +0200 Subject: [PATCH 2/4] fix(rounding): rounded logarithms and exponential reach 38 places rounded_ln, rounded_log10 and rounded_exp now take places from -38 to 38, as checked_round does. decide_rounding scales an end by up to 10^38, so the transcendental kernel widens from 384 to 448 bits: e^-80 to whole 10^38s scales the denominator 2^(192+116) past 2^384. The exponential's shortcut below -43 assumed e^x was below a quarter of the last kept unit at any place accepted, which holds only up to 18 places: at 22 places e^-50 is 2e-22, not 0. The shortcut now starts below -89, the argument whose exponential is below a quarter of 10^-38, and the kernel's reduction takes arguments down to -89 (m up to 129); the error bound and the width bounds in the kernel's comment are restated for it. A logarithm's enclosure is at least 1.6e-37 wide, so a logarithm rounded at 37 or 38 places is Overflow, never a guess; up to 30 places the reference table decides every row. The documented constant-evaluation costs are re-measured: 10 to 15 percent more per enclosure. Refs #34 Signed-off-by: Christian Parpart --- include/formula-cpp/detail/transcendental.hpp | 58 +++++----- .../formula-cpp/rounded_transcendental.hpp | 35 ++++-- test/rounded_transcendental_tests.cpp | 101 +++++++++++++++--- test/transcendental_tests.cpp | 29 +++-- 4 files changed, 159 insertions(+), 64 deletions(-) diff --git a/include/formula-cpp/detail/transcendental.hpp b/include/formula-cpp/detail/transcendental.hpp index bd088c3..18bb41b 100644 --- a/include/formula-cpp/detail/transcendental.hpp +++ b/include/formula-cpp/detail/transcendental.hpp @@ -3,7 +3,7 @@ /// @file /// An integer kernel that encloses the natural logarithm, the decimal logarithm and the exponential of a -/// rational: two ends between which the value certainly lies, computed in 384-bit fixed point with only +/// rational: two ends between which the value certainly lies, computed in 448-bit fixed point with only /// integer operations the language defines exactly, so that no floating-point mode enters and the same /// inputs are meant to give the same bits, at compile time and at run time. It is what the rounded forms /// (`rounded_transcendental.hpp`) round: when both ends round to the same decimal, that is the rounding of @@ -14,12 +14,12 @@ /// - **The argument** is a/b in lowest terms, as a `Rational` holds it: |a| <= 2^127 (the magnitude of /// `Int128`'s minimum) and 0 < b < 2^127. The kernel reads the two magnitudes as `UInt128` and never narrows /// them. -/// - **Two fixed points.** A value v is an integer V in `WideUnsigned<12>` (384 bits), V = floor(v 2^F): F = +/// - **Two fixed points.** A value v is an integer V in `WideUnsigned<14>` (448 bits), V = floor(v 2^F): F = /// 128 fraction bits for a logarithm, whose value is below 89, and F = 192 for an exponential, whose value /// can be as wide as a `Rational` -- up to 2^127 -- and must still be enclosed more narrowly than its last /// kept unit. Every operation truncates a non-negative value, so a lower bound stays one; each upper bound /// is the lower bound plus a slack derived here. No ``, no floating point, no intrinsics, **and no -/// call of the general `divmod`** (too costly in a constant evaluation over 384 bits): the scaled quotients +/// call of the general `divmod`** (too costly in a constant evaluation over 448 bits): the scaled quotients /// are long divisions over `UInt128` words, whose whole part is `u128_divmod` -- the compiler's own 128-bit /// integer where it has one, portable code elsewhere, the same quotient either way -- k is a binary search, /// and the series' divisors are below 2^32 (`divmod_small`). @@ -38,33 +38,34 @@ /// upper = floor(upper_ln (M + 1) / 2^128) + 1, under 254 · 0.44 + 89 + 2 < 203 units apart, since /// |ln(a/b)| <= ln(2^127) < 89. upper_ln is below 89 · 2^128 + 254 < 2^135 and M + 1 below 2^127, so /// log10's widest product, upper_ln (M + 1), is below 2^262. -/// - **exp(x)**, x = a/b != 0, -43 <= x <= 887/10 (the rounded forms answer outside it), in F = 192 bits. +/// - **exp(x)**, x = a/b != 0, -89 <= x <= 887/10 (the rounded forms answer outside it), in F = 192 bits. /// X = floor(|x| 2^192), exact or one below. L' = floor(ln 2 2^192). For x > 0, k is the largest integer /// in [0, 127] with k (L' + 1) <= X (128 ln 2 > 887/10 bounds it), and R = X - k (L' + 1) <= r 2^192 for -/// r = x - k ln 2; for x < 0, m is the smallest in [1, 63] with m L' >= X' (X' = X + 1 when X is inexact; -/// 63 ln 2 > 43 bounds it), R = m L' - X' <= r 2^192 for r = m ln 2 - |x|, and exp(x) = 2^-m exp(r). -/// Either way 0 <= R < L' + 1, so r < ln 2 + 2^-184, and r 2^192 - R < 128: one unit for X, and one per -/// multiple of ln 2's unit, k <= 127 or m <= 63. E = sum T_j, T_0 = 2^192, +/// r = x - k ln 2; for x < 0, m is the smallest in [1, 129] with m L' >= X' (X' = X + 1 when X is inexact; +/// 129 ln 2 > 89 bounds it), R = m L' - X' <= r 2^192 for r = m ln 2 - |x|, and exp(x) = 2^-m exp(r). +/// Either way 0 <= R < L' + 1, so r < ln 2 + 2^-184, and r 2^192 - R < 130: one unit for X, and one per +/// multiple of ln 2's unit, k <= 127 or m <= 129. E = sum T_j, T_0 = 2^192, /// T_j = floor(floor(T_{j-1} R / 2^192) / j), until T_j = 0 (by j = 43 for r < 0.7; `TaylorTermLimit` = 50 /// refuses a longer one), so E <= exp(R 2^-192) 2^192 <= exp(r) 2^192. T_{j-1} <= 2^192 and R <= L' < -/// 0.7 2^192, so the product T_{j-1} R is below 2^384: the kernel's widest, and the reason 192 fraction -/// bits are the most 384 bits allow. Each T_j is short by e_j < e_{j-1} r / j + 1 < 2, and the tail after -/// the last term is below 3, so exp(R 2^-192) 2^192 - E < 2 · 50 + 3; the 128 units of r add less than -/// 2 · 1.0001 · 128 < 257. So exp(r) 2^192 < E + 360 <= E + `ExponentialSlack` = 512: lower = E 2^k / 2^192, +/// 0.7 2^192, so the product T_{j-1} R is below 2^384. Each T_j is short by e_j < e_{j-1} r / j + 1 < 2, +/// and the tail after the last term is below 3, so exp(R 2^-192) 2^192 - E < 2 · 50 + 3; the 130 units of +/// r add less than 2 · 1.0001 · 130 < 261. So exp(r) 2^192 < E + 364 <= E + `ExponentialSlack` = 512: +/// lower = E 2^k / 2^192, /// upper = (E + 512) 2^k / 2^192 (for x < 0, denominator 2^(192+m)). Since E >= 2^192, the ends are at /// most 2^-183 of the value apart: at the largest result a `Rational` holds, 2^127 last kept units, that is -/// 2^-56 of one unit. E < 2^193, so the widest numerator, (E + 512) 2^127, is below 2^321; -/// `decide_rounding`'s scaling by up to 10^18 keeps it below 2^381. +/// 2^-56 of one unit. E < 2^193, so the widest numerator, (E + 512) 2^127, is below 2^321, and the widest +/// denominator, 2^(192+129), is 2^321 exactly; `decide_rounding`'s scaling by up to 10^38 (below 2^126.3) +/// keeps each below 2^448. /// - **Undecided.** When the two ends round differently, `decide_rounding` answers `Overflow`: the rounding /// needs more bits than the kernel holds. The rounded decimal exists, so `Inexact` would be wrong, and a new /// error enumerator would break `describe()` and consumers' switches. /// /// ## What it costs /// -/// Measured on cl 19.51.36257, whose default constant-evaluation budget measured about 1 049 000 steps: one -/// enclosure costs between 23 100 and 43 700 steps (ln 3: 23 100; ln (2^127 - 1) / 2^126: 43 700; log10 7: 28 600; -/// exp 1: 32 900; exp -43: 38 200; exp 88: 38 400), and a whole rounding of log10 2 to 3 places, with -/// `decide_rounding`, about 242 300. +/// Measured on cl 19.51.36260, whose default constant-evaluation budget measured about 1 049 000 steps: one +/// enclosure costs between 26 400 and 48 100 steps (ln 3: 26 400; ln (2^127 - 1) / 2^126: 48 100; log10 7: 32 600; +/// exp 1: 37 600; exp -43: 43 900; exp -89: 40 100; exp 88: 44 000), a whole rounding of log10 2 to 3 places, with +/// `decide_rounding`, about 279 200, of ln 2 to 30 places about 353 900, and of exp 1 to 30 places about 505 800. #include #include @@ -77,11 +78,12 @@ namespace formula::detail { -/// The kernel's width, 384 bits: room for its widest product, the exponential series' T_{j-1} R (under -/// 2^384), for the logarithm's widest product, upper_ln (M + 1) (under 2^262), for the exponential's widest -/// numerator (under 2^321), and for `decide_rounding`'s scaling by up to 10^18 of every end it is handed -/// (under 2^381). -inline constexpr std::size_t KernelLimbs = 12; +/// The kernel's width, 448 bits: room for `decide_rounding`'s scaling by up to 10^38 -- the most decimal +/// places rounding takes, `MaximumDecimalPlaces` -- of every end it is handed (under 2^448), the widest of +/// them the exponential's numerator or denominator (at most 2^321): e^-80 to whole 10^38s scales the +/// denominator 2^(192+116) by 10^38, past 2^384; for the exponential series' T_{j-1} R (under 2^384); and +/// for the logarithm's widest product, upper_ln (M + 1) (under 2^262). +inline constexpr std::size_t KernelLimbs = 14; /// A value in the kernel's fixed point. using KernelWord = WideUnsigned; /// How many of a logarithm's fixed-point bits are fraction. @@ -109,14 +111,14 @@ struct Enclosure WideRatio upper; }; -/// @p highHalf * 2^64 + @p lowHalf as a kernel word. 128 bits in 384: neither step can overflow. +/// @p highHalf * 2^64 + @p lowHalf as a kernel word. 128 bits in 448: neither step can overflow. [[nodiscard]] constexpr KernelWord kernel_word(std::uint64_t highHalf, std::uint64_t lowHalf) noexcept { return *add_checked_or_none(*shift_left_checked_or_none(KernelWord::from_u64(highHalf), 64), KernelWord::from_u64(lowHalf)); } -/// @p topWord * 2^128 + @p middleWord * 2^64 + @p bottomWord as a kernel word. 192 bits in 384: no step can +/// @p topWord * 2^128 + @p middleWord * 2^64 + @p bottomWord as a kernel word. 192 bits in 448: no step can /// overflow. [[nodiscard]] constexpr KernelWord kernel_word(std::uint64_t topWord, std::uint64_t middleWord, @@ -331,7 +333,7 @@ struct LogarithmMagnitude return signed_enclosure(shift_right(*lowerProduct, KernelFractionBits), *farther, natural->sign); } -/// exp(@p argument), enclosed -- see the file comment. @pre @p argument != 0 and -43 <= @p argument <= 887/10. +/// exp(@p argument), enclosed -- see the file comment. @pre @p argument != 0 and -89 <= @p argument <= 887/10. [[nodiscard]] constexpr std::optional exponential_enclosure(Rational argument) noexcept { bool const negative = argument.sign() < 0; @@ -362,14 +364,14 @@ struct LogarithmMagnitude } else { - // The smallest m in [1, 63] with m L' >= X' (X rounded up): 63 ln 2 > 43 >= |x|. + // The smallest m in [1, 129] with m L' >= X' (X rounded up): 129 ln 2 > 89 >= |x|. std::optional const roundedUp = fixedMagnitude.exact ? std::optional { fixedMagnitude.below } : add_small_checked_or_none(fixedMagnitude.below, 1U); if (!roundedUp) return std::nullopt; std::uint32_t below = 0; - std::uint32_t above = 63; + std::uint32_t above = 129; while (below + 1 < above) { std::uint32_t const middle = (below + above) / 2; diff --git a/include/formula-cpp/rounded_transcendental.hpp b/include/formula-cpp/rounded_transcendental.hpp index b1ec3f1..c5ada42 100644 --- a/include/formula-cpp/rounded_transcendental.hpp +++ b/include/formula-cpp/rounded_transcendental.hpp @@ -45,15 +45,26 @@ namespace detail /// largest `Rational` at every places. Above it the answer is `Overflow` without the kernel. inline constexpr Rational ExponentialArgumentCap { 887, 10 }; + /// The smallest argument the exponential's kernel takes, -89: e^-89 = 2.2 * 10^-39 is below a quarter of + /// 10^-38, the finest last kept unit (`MaximumDecimalPlaces`), so below it e^x is below a quarter of the + /// last kept unit at every places accepted and its rounding is decided without the kernel. 129 ln 2 > + /// 89 bounds the kernel's reduction (`detail/transcendental.hpp`). + inline constexpr Rational ExponentialArgumentFloor { -89 }; + /// @p F of @p argument, rounded to @p places decimal places under @p roundingMode -- the correctly /// rounded decimal of the true value, rational or not. The decision, in order: a logarithm of zero or - /// below is `DomainError`; places outside -18...18 are `Overflow`, as for `checked_round`; a special - /// point -- the only values that can tie -- goes to `checked_round`; the exponential of more than 887/10 - /// is `Overflow` (`ExponentialArgumentCap`: no such value fits a `Rational`), and of less than -43 is - /// below a quarter of the last kept unit at any places accepted, so 0, or one unit under `Ceiling` and - /// `AwayFromZero`; everything else is the kernel's enclosure, rounded by `decide_rounding`, which - /// answers `Overflow` when the kept integer does not fit -- e^88.5, at every places -- and when the two - /// ends round differently. The kernel takes every argument a `Rational` holds. The special points are + /// below is `DomainError`; places outside -`MaximumDecimalPlaces`...`MaximumDecimalPlaces` are + /// `Overflow`, as for `checked_round`; a special point -- the only values that can tie -- goes to + /// `checked_round`; the exponential of more than 887/10 is `Overflow` (`ExponentialArgumentCap`: no such + /// value fits a `Rational`), and of less than -89 (`ExponentialArgumentFloor`) is below a quarter of the + /// last kept unit at any places accepted, so 0, or one unit under `Ceiling` and `AwayFromZero`; + /// everything else is the kernel's enclosure, rounded by `decide_rounding`, which answers `Overflow` + /// when the kept integer does not fit -- e^88.5, at every places -- and when the two ends round + /// differently. A logarithm's two ends are at least 1.6 * 10^-37 apart, more than one unit at 37 places, + /// so a logarithm the kernel computes is always `Overflow` at 37 or 38 places; its exact points -- ln 1, + /// log10 10^k -- are answered without the kernel, wherever their value fits: ln 1 at every places, log10 10 + /// and log10 1/10 at 38, log10 10^k for |k| <= 17 at 37. The kernel takes every argument a + /// `Rational` holds. The special points are /// `RepFunctions`'s, through `transcendental_of`: their value, and `Inexact` elsewhere. template [[nodiscard]] constexpr std::expected rounded_transcendental( @@ -64,7 +75,7 @@ namespace detail if (argument.sign() <= 0) return std::unexpected { ArithmeticError::DomainError }; } - if (places.value > 18 || places.value < -18) + if (places.value > MaximumDecimalPlaces || places.value < -MaximumDecimalPlaces) return std::unexpected { ArithmeticError::Overflow }; std::expected const exact = transcendental_of(argument); @@ -82,7 +93,7 @@ namespace detail { if (argument > ExponentialArgumentCap) return std::unexpected { ArithmeticError::Overflow }; - if (argument < Rational { -43 }) + if (argument < ExponentialArgumentFloor) { if (roundingMode == RoundingMode::Ceiling || roundingMode == RoundingMode::AwayFromZero) return Rational::from_decimal(1, -places.value); @@ -127,7 +138,9 @@ struct RoundedTranscendentalNode: NodeBase /// places: `rounded_ln(var / var)`. /// /// The integer kernel (`detail/transcendental.hpp`) takes every argument a `Rational` holds: ln (2^127 - 1) -/// is 88.029691931113054295 at 18 places, under `Floor`. A rounding the kernel cannot decide is `Overflow`. +/// is 88.029691931113054295 at 18 places, under `Floor`. A rounding the kernel cannot decide is `Overflow`: +/// so is every logarithm it computes at 37 or 38 places, where its enclosure is wider than the last kept +/// unit. ln 1 is answered exactly before the kernel is asked, at any places. template [[nodiscard]] constexpr auto rounded_ln(Operand operand) noexcept { @@ -165,7 +178,7 @@ template /// The exponential of `operand`, rounded exactly to `Places` decimal places. /// /// Checked in this order: an argument above 887/10 is `Overflow`, since e^x is then past the largest -/// `Rational`; one below -43 is 0, or one last kept unit under `Ceiling` and `AwayFromZero`, whatever its +/// `Rational`; one below -89 is 0, or one last kept unit under `Ceiling` and `AwayFromZero`, whatever its /// width, so exp(-2^70) is 0; any other goes to the integer kernel (`detail/transcendental.hpp`), which /// answers wherever the result fits the declared places -- e^45 to 18 places, e^88 to whole units -- and is /// `Overflow` where it does not, as e^88.5 is at every places. A rounding the kernel cannot decide is diff --git a/test/rounded_transcendental_tests.cpp b/test/rounded_transcendental_tests.cpp index e2124d4..43795d6 100644 --- a/test/rounded_transcendental_tests.cpp +++ b/test/rounded_transcendental_tests.cpp @@ -21,6 +21,7 @@ using formula::DecimalPlaces; using formula::Rational; using formula::RoundingMode; using formula::var; +using namespace formula::literals; struct Ratio: formula::Quantity { @@ -173,17 +174,85 @@ TEST_CASE("rounded_transcendental: an exponential answers wherever it fits a Rat TEST_CASE("rounded_transcendental: a tiny exponential is zero or one unit by mode", "[rounded_transcendental]") { - // exp -50 = 1.9 * 10^-22, below a quarter of any last kept unit: 0, except one unit upward. - STATIC_REQUIRE(expAt(Rational { -50 }) == Rational {}); - STATIC_REQUIRE(expAt(Rational { -50 }) == Rational {}); - STATIC_REQUIRE(expAt(Rational { -50 }) == Rational {}); - STATIC_REQUIRE(expAt(Rational { -50 }) == Rational { 1, 1000 }); - STATIC_REQUIRE(expAt(Rational { -50 }) == Rational { 1, 1000 }); - STATIC_REQUIRE(expAt(Rational { -50 }) == Rational { 100 }); - STATIC_REQUIRE(expAt(Rational { -50 }) == Rational {}); - // exp -43 = 2.1 * 10^-19, through the kernel at its lowest power of two, 2^-63. + // exp -100 = 3.7 * 10^-44, below a quarter of any last kept unit: 0, except one unit upward. + STATIC_REQUIRE(expAt(Rational { -100 }) == Rational {}); + STATIC_REQUIRE(expAt(Rational { -100 }) == Rational {}); + STATIC_REQUIRE(expAt(Rational { -100 }) == Rational {}); + STATIC_REQUIRE(expAt(Rational { -100 }) == Rational { 1, 1000 }); + STATIC_REQUIRE(expAt(Rational { -100 }) == Rational { 1, 1000 }); + STATIC_REQUIRE(expAt(Rational { -100 }) == Rational { 100 }); + STATIC_REQUIRE(expAt(Rational { -100 }) == Rational {}); + // And at 38 places, the last: still below a quarter of 10^-38. + STATIC_REQUIRE(expAt(Rational { -100 }) == Rational {}); + STATIC_REQUIRE(expAt(Rational { -100 }) == 1e-38_r); + // exp -43 = 2.1 * 10^-19, through the kernel. CHECK(expAt(Rational { -43 }) == Rational {}); CHECK(expAt(Rational { -43 }) == Rational { 1, 1'000'000'000'000'000'000 }); + // exp -50 = 1.9287498479639177830 * 10^-22 is no longer below a quarter of the last kept unit past 21 + // places: through the kernel, and its digits. + CHECK(expAt(Rational { -50 }) == 2e-22_r); + CHECK(expAt(Rational { -50 }) + == 0.00000000000000000000019287498479639178_r); + // exp -89 = 2.2 * 10^-39, the smallest argument the kernel takes, at its lowest power of two, 2^-129: + // below a quarter of 10^-38, 2.5 * 10^-39, so 0 to nearest at 38 places and one unit up under Ceiling. + CHECK(expAt(Rational { -89 }) == Rational {}); + CHECK(expAt(Rational { -89 }) == 1e-38_r); + CHECK(expAt(Rational { -89 }) == Rational {}); +} + +TEST_CASE("rounded_transcendental: the exponential's shortcut starts where the kernel stops and no earlier", + "[rounded_transcendental]") +{ + // The shortcut below -89 (ExponentialArgumentFloor) is right only where e^x is below a quarter of + // 10^-38, and the kernel reaches only to -129 ln 2 = -89.416. Pinned from both sides. + // exp -88.19 = 5.0069 * 10^-39, just past half a unit at 38 places: through the kernel, 10^-38 to + // nearest. A floor above -88.19 would answer 0. + CHECK(expAt(Rational { -8819, 100 }) == 1e-38_r); + // exp -70 = 3.9754497359086468 * 10^-31, its digits at 38 places; a floor above -70 would answer 0. + CHECK(expAt(Rational { -70 }) + == 0.00000000000000000000000000000039754497_r); + // exp -89.5 = 1.35 * 10^-39, past the kernel's reach: the shortcut's 0, and one unit under Ceiling. A + // floor below -89.416 would send it to the kernel, which answers Overflow. + STATIC_REQUIRE(expAt(Rational { -895, 10 }) == Rational {}); + STATIC_REQUIRE(expAt(Rational { -895, 10 }) == 1e-38_r); +} + +TEST_CASE("rounded_transcendental: places reach 38 either way, as checked_round's do", "[rounded_transcendental]") +{ + // Through the kernel, so at run time. ln 2 = 0.693147180559945309417232121458|1765..., to 30 places. + CHECK(lnAt(Rational { 2 }) == 0.693147180559945309417232121458_r); + CHECK(lnAt(Rational { 2 }) == 0.693147180559945309417232121459_r); + // log10 7 = 0.845098040014256830712216258592|6363..., to 30 places. + CHECK(log10At(Rational { 7 }) == 0.845098040014256830712216258593_r); + // exp 1 = 2.7182818284590452353602874713526624977|572..., to 37 places; at 38 the kept integer, + // 2.7 * 10^38, is past 2^127. + CHECK(expAt(Rational { 1 }) + == 2.7182818284590452353602874713526624978_r); + CHECK(expAt(Rational { 1 }) == 2.7182818284590452353602874713526624977_r); + CHECK(expAt(Rational { 1 }) == overflow); + // A logarithm's two ends are at least 1.6 * 10^-37 apart, wider than a unit at 37 places: there and at + // 38 it cannot place ln 2 between two units, and the answer is Overflow, never a guess. + CHECK(lnAt(Rational { 2 }) == overflow); + CHECK(lnAt(Rational { 2 }) == overflow); + CHECK(log10At(Rational { 7 }) == overflow); + // The exact points never reach the kernel, and answer wherever their value fits: ln 1 at every places, + // log10 10 at 38 and log10 1/10 at 37. log10 100 at 38 places is 2 * 10^38 units, and log10 10^18 at 37 + // is 1.8 * 10^38: past 2^127 both, so Overflow. + STATIC_REQUIRE(lnAt(Rational { 1 }) == Rational {}); + STATIC_REQUIRE(log10At(Rational { 10 }) == Rational { 1 }); + STATIC_REQUIRE(log10At(Rational { 1, 10 }) == Rational { -1 }); + STATIC_REQUIRE(log10At(Rational { 100 }) == overflow); + STATIC_REQUIRE(log10At(Rational { 1'000'000'000'000'000'000 }) + == overflow); + // To whole 10^38s at -38 places: exp 88 = 1.65 * 10^38 is 10^38 under Floor, and 2 * 10^38, past + // 2^127, to nearest; ln 2 is 10^38 under Ceiling. + CHECK(expAt(Rational { 88 }) == 1e38_r); + CHECK(expAt(Rational { 88 }) == overflow); + CHECK(lnAt(Rational { 2 }) == 1e38_r); + // exp -80 = 1.8 * 10^-35, whose ends lie over 2^(192+116): to whole 10^38s that denominator times + // 10^38 is past 2^434, which the kernel's 448 bits hold. 0, and 10^38 under Ceiling. + CHECK(expAt(Rational { -80 }) == Rational {}); + CHECK(expAt(Rational { -80 }) == 1e38_r); } TEST_CASE("rounded_transcendental: absence and failures come first and in order", "[rounded_transcendental]") @@ -206,19 +275,19 @@ TEST_CASE("rounded_transcendental: absence and failures come first and in order" // A logarithm of zero or below; before places out of range. STATIC_REQUIRE(lnAt(Rational {}) == domainError); STATIC_REQUIRE(log10At(Rational { -1 }) == domainError); - STATIC_REQUIRE(lnAt(Rational { -1 }) == domainError); - STATIC_REQUIRE(lnAt(Rational { 2 }) == overflow); - STATIC_REQUIRE(lnAt(Rational { 1 }) == overflow); + STATIC_REQUIRE(lnAt(Rational { -1 }) == domainError); + STATIC_REQUIRE(lnAt(Rational { 2 }) == overflow); + STATIC_REQUIRE(lnAt(Rational { 1 }) == overflow); // Places out of range before the tiny rule: without that check the tiny rule would answer 0 in the modes // that round down or to nearest (only the two upward modes overflow on their own). - STATIC_REQUIRE(expAt(Rational { -50 }) == overflow); - STATIC_REQUIRE(expAt(Rational { -50 }) == overflow); - STATIC_REQUIRE(expAt(Rational { -50 }) == overflow); + STATIC_REQUIRE(expAt(Rational { -100 }) == overflow); + STATIC_REQUIRE(expAt(Rational { -100 }) == overflow); + STATIC_REQUIRE(expAt(Rational { -100 }) == overflow); // A wide power of ten is a special point, answered before the kernel is asked: log10 10^30 is 30. constexpr Rational::Int tenToFifteen = 1'000'000'000'000'000; STATIC_REQUIRE(log10At(Rational { tenToFifteen * tenToFifteen }) == Rational { 30 }); - // The rule below -43 comes first: exp -2^70 is 0 without the kernel. + // The rule below -89 comes first: exp -2^70 is 0 without the kernel. STATIC_REQUIRE(expAt(Rational { -(Rational::Int { 1 } << 70) }) == Rational {}); } diff --git a/test/transcendental_tests.cpp b/test/transcendental_tests.cpp index fc568a4..b7e2caa 100644 --- a/test/transcendental_tests.cpp +++ b/test/transcendental_tests.cpp @@ -148,7 +148,7 @@ constexpr Rational::Int largestInt = std::numeric_limits::max(); constexpr Rational::Int smallestInt = std::numeric_limits::min(); // -2^127 // clang-format off -constexpr std::array references { { +constexpr std::array references { { { Transcendental::NaturalLogarithm, 2, 1, false, "6931471805599453094172321214581765680755", 40 }, { Transcendental::NaturalLogarithm, 1, 2, true, "6931471805599453094172321214581765680755", 40 }, { Transcendental::NaturalLogarithm, 3, 1, false, "1098612288668109691395245236922525704647", 39 }, @@ -181,6 +181,9 @@ constexpr std::array references { { { Transcendental::Exponential, -1, 1, false, "3678794411714423215955237701614608674458", 40 }, { Transcendental::Exponential, 43, 1, false, "4727839468229346561474457562744280370819", 21 }, { Transcendental::Exponential, -43, 1, false, "2115131037591080486631401007022651470196", 58 }, + { Transcendental::Exponential, -50, 1, false, "1928749847963917783017342816527012574752", 61 }, + { Transcendental::Exponential, -88, 1, false, "6054601895401185884531860533810599287768", 78 }, + { Transcendental::Exponential, -89, 1, false, "2227363561795743739222909281640826188721", 78 }, { Transcendental::Exponential, 1, 2, false, "1648721270700128146848650787814163571653", 39 }, { Transcendental::Exponential, -1, 2, false, "6065306597126334236037995349911804534419", 40 }, { Transcendental::Exponential, 10, 1, false, "2202646579480671651695790064528424436635", 35 }, @@ -264,9 +267,10 @@ TEST_CASE("transcendental kernel: the kernel re-derives its stored constants fro TEST_CASE("transcendental kernel: every reference value is enclosed and rounds as the reference does in every mode", "[transcendental]") { - constexpr std::array placesTried { -2, -1, 0, 1, 2, 4, 9, 17, 18 }; + constexpr std::array placesTried { -38, -2, -1, 0, 1, 2, 4, 9, 17, 18, 30, 37 }; std::size_t compared = 0; std::size_t undecided = 0; + std::size_t undecidedLogarithms = 0; for (std::size_t rowAt = 0; rowAt < references.size(); ++rowAt) { Reference const& row = references[rowAt]; @@ -291,18 +295,25 @@ TEST_CASE("transcendental kernel: every reference value is enclosed and rounds a if (!decided.has_value() && referenceDecided.has_value()) { CHECK(decided.error() == formula::ArithmeticError::Overflow); - ++undecided; + if (row.function != Transcendental::Exponential && places == 37) + ++undecidedLogarithms; + else + ++undecided; } else CHECK(decided == referenceDecided); ++compared; } } - // 51 rows, 9 places, 7 modes: a loop over nothing fails here. - REQUIRE(compared == 3213); - // None: an exponential's 192 fraction bits and a logarithm's 128 decide every row the reference's 40 - // digits decide. At 128 fraction bits the exponentials of 43 to 44 at 17 and 18 places were not. + // 54 rows, 12 places, 7 modes: a loop over nothing fails here. + REQUIRE(compared == 4536); + // None but the logarithms at 37 places: an exponential's 192 fraction bits decide every row the + // reference's 40 digits decide, at every places tried, and a logarithm's 128 at every places up to 30. + // At 128 fraction bits the exponentials of 43 to 44 at 17 and 18 places were not. CHECK(undecided == 0); + // A logarithm's ends are at least 1.6 * 10^-37 apart, wider than a unit at 37 places: every logarithm + // the reference decides there, 147 rows and modes, is Overflow. + CHECK(undecidedLogarithms == 147); // Three of them written out, so that a reader sees the digits. CHECK(kernel_rounding(Transcendental::NaturalLogarithm, Rational { 2 }, 18, RoundingMode::Floor) == Rational::from_decimal(693'147'180'559'945'309, -18)); @@ -362,7 +373,7 @@ TEST_CASE("transcendental kernel: a scaled quotient of 128-bit operands carries "[transcendental]") { // The divisor 2^128 - 1 leaves a remainder above 2^127, whose doubling passes 2^128: the shifted-out bit - // must still count. Checked against the general long division of the 384-bit word, a different route. + // must still count. Checked against the general long division of the 448-bit word, a different route. detail::UInt128 const divisor { ~std::uint64_t { 0 }, ~std::uint64_t { 0 } }; detail::UInt128 const dividend { std::uint64_t { 1 } << 63, 5 }; for (std::size_t const fractionBits: { std::size_t { 128 }, std::size_t { 192 } }) @@ -388,7 +399,7 @@ TEST_CASE("transcendental kernel: the kernel answers at compile time", "[transce { // The one deliberate compile-time check of the kernel. A whole rounding -- the enclosure and // decide_rounding -- in one constant evaluation, measured at about - // 242 300 steps on cl 19.51.36257, against a default budget of about 1 049 000. Every other check that runs the kernel + // 279 200 steps on cl 19.51.36260, against a default budget of about 1 049 000. Every other check that runs the kernel // runs at run time. STATIC_REQUIRE(kernel_rounding(Transcendental::DecimalLogarithm, Rational { 2 }, 3, RoundingMode::HalfEven) == Rational { 301, 1000 }); From 43751144f9100cb7a562179df960f2e7bc263eb3 Mon Sep 17 00:00:00 2001 From: Christian Parpart Date: Thu, 8 Oct 2026 10:02:13 +0200 Subject: [PATCH 3/4] docs(text): state text output's 18 places as its own limit Text output keeps refusing more than 18 decimal places either way: fixed 64-bit digit buffers, long division and the compile-time format-spec checks. Its comments said that limit was checked_round's, or the -18 to 18 DecimalPlaces spans; rounding now reaches 38, so they state it as the text layer's own. format.hpp's two bare 18s name ExactDecimalPlaces. No behaviour changes. Refs #34 Signed-off-by: Christian Parpart --- docs/display.md | 4 ++-- include/formula-cpp/format.hpp | 12 ++++++------ include/formula-cpp/number_text.hpp | 18 ++++++++++-------- include/formula-cpp/render.hpp | 2 +- include/formula-cpp/trace_render.hpp | 3 ++- test/calculation_trace_tests.cpp | 6 +++--- test/format_tests.cpp | 4 ++-- test/number_text_tests.cpp | 4 ++-- test/trace_render_tests.cpp | 4 ++-- 9 files changed, 30 insertions(+), 27 deletions(-) diff --git a/docs/display.md b/docs/display.md index b956f69..aca774b 100644 --- a/docs/display.md +++ b/docs/display.md @@ -653,8 +653,8 @@ before the call stack of the evaluation: ``` test\negative\format_places_without_mode.cpp(17): error C7595: 'std::basic_format_string::basic_format_string': call to immediate function is not a constant expression -include\formula-cpp/format.hpp(347): note: failure was caused by call of undefined function or one not declared 'constexpr' -include\formula-cpp/format.hpp(347): note: see usage of 'formula::detail::formula_number_format_needs_a_rounding_mode' +include\formula-cpp/format.hpp(348): note: failure was caused by call of undefined function or one not declared 'constexpr' +include\formula-cpp/format.hpp(348): note: see usage of 'formula::detail::formula_number_format_needs_a_rounding_mode' ``` clang and g++ name the same function, in their own words. diff --git a/include/formula-cpp/format.hpp b/include/formula-cpp/format.hpp index e98479f..cd5a815 100644 --- a/include/formula-cpp/format.hpp +++ b/include/formula-cpp/format.hpp @@ -87,7 +87,8 @@ namespace formula::detail /// Refuses a number format that rounds to more than 18 decimal places, or /// that rounds a `Measured` at its unit's decimals when those lie outside -18 -/// to 18. Not `constexpr`, for the reason +/// to 18: text stops at 18 places either way (`detail::ExactDecimalPlaces`), +/// though rounding itself reaches 38. Not `constexpr`, for the reason /// `formula_number_format_needs_a_rounding_mode` gives. /// @throws std::format_error always. [[noreturn]] inline void formula_number_format_places_out_of_range() @@ -334,7 +335,7 @@ inline constexpr RoundingModeName RoundingModeNames[] { { if (!beyondRange) placesAsked = placesAsked * 10 + (specText[at] - '0'); - beyondRange = beyondRange || placesAsked > 18; + beyondRange = beyondRange || placesAsked > ExactDecimalPlaces; ++at; } if (beyondRange) @@ -442,7 +443,7 @@ template /// Reads a `Measured` or `Outcome` replacement field's spec, as /// `parse_number_format_field` does, and refuses `~Mode` without `.N` when /// the unit @p Q's value is shown in (`shown_unit_of`) declares decimals -/// outside the -18 to 18 that `DecimalPlaces` spans. +/// outside the -18 to 18 text is written at (`ExactDecimalPlaces`). template [[nodiscard]] constexpr std::format_parse_context::iterator parse_measured_format_field( std::format_parse_context& parseContext, NumberFormatSpec& parsed) @@ -450,7 +451,7 @@ template auto const specEnd = parse_number_format_field(parseContext, parsed); constexpr int declaredPlaces = shown_unit_of(Describe::unit, Describe::dimension).decimals; if (parsed.body == NumberFormatBody::Approximated && !parsed.places.has_value() - && (declaredPlaces > 18 || declaredPlaces < -18)) + && (declaredPlaces > ExactDecimalPlaces || declaredPlaces < -ExactDecimalPlaces)) formula_number_format_places_out_of_range(); return specEnd; } @@ -749,8 +750,7 @@ struct formatter, char> { /// Reads the spec up to its `}`, as `formatter` does. `~Mode` /// without `.N` rounds at `Q`'s declared decimals, so a unit whose declared - /// decimals lie outside the -18 to 18 that `DecimalPlaces` spans is - /// refused here. + /// decimals lie outside the -18 to 18 text is written at is refused here. constexpr auto parse(std::format_parse_context& parseContext) { return formula::detail::parse_measured_format_field(parseContext, _spec); diff --git a/include/formula-cpp/number_text.hpp b/include/formula-cpp/number_text.hpp index 6c96cc7..68640dc 100644 --- a/include/formula-cpp/number_text.hpp +++ b/include/formula-cpp/number_text.hpp @@ -270,13 +270,15 @@ namespace detail static_assert(LongestNumberText <= NumberTextCapacity, "formula: NumberTextCapacity is too small for the longest number this library spells"); - /// 10^18, the finest scale this header writes a decimal on: rounding's - /// places stop at 18, as `DecimalPlaces` does, though `Rational::Int` - /// holds powers of ten up to 10^38. + /// 10^18, the finest scale this header writes a decimal on: text stops at + /// 18 places, the most a 64-bit fraction part holds, though rounding + /// reaches 38 (`MaximumDecimalPlaces`) and `Rational::Int` holds powers of + /// ten up to 10^38. inline constexpr std::uint64_t ExactDecimalScale = 1'000'000'000'000'000'000ULL; - /// The places `ExactDecimalScale` spans -- also the largest - /// `DecimalPlaces` `checked_round` accepts. + /// The places `ExactDecimalScale` spans: the most decimal places text is + /// written at, either way -- the text layer's own limit, refused here and + /// in `format.hpp`, not `checked_round`'s. inline constexpr int ExactDecimalPlaces = 18; /// @p shownValue as an exact decimal with at least @p minimumPlaces @@ -554,9 +556,9 @@ namespace detail /// is `0.600`, `123/1000` in a unit of 1 decimal is `0.123`. /// /// @return `Overflow` when the style reads the unit's decimals -- it pads, or -/// it approximates -- and they lie outside the -18 to 18 that -/// `DecimalPlaces` spans, whatever the value; otherwise any error of -/// `checked_decimal_text`. +/// it approximates -- and they lie outside the -18 to 18 text is +/// written at (`checked_decimal_text`), whatever the value; otherwise +/// any error of `checked_decimal_text`. [[nodiscard]] constexpr std::expected checked_number_text(Rational shownValue, NumberStyle shownStyle, Unit const& shownIn) noexcept diff --git a/include/formula-cpp/render.hpp b/include/formula-cpp/render.hpp index c6460c8..bc5475f 100644 --- a/include/formula-cpp/render.hpp +++ b/include/formula-cpp/render.hpp @@ -349,7 +349,7 @@ namespace detail /// For a helper with no clause to say a number is not shown: a table's /// bound, an envelope's limit. The fraction style never fails. Any other /// can -- a padded or approximated number in a unit whose declared - /// decimals lie outside the -18 to 18 `DecimalPlaces` spans, for one -- + /// decimals lie outside the -18 to 18 text is written at, for one -- /// and the fraction is then the one text still exact. [[nodiscard]] inline std::string styled_number_text(Rational shownNumber, NumberStyle numberStyle, Unit const& shownIn) { diff --git a/include/formula-cpp/trace_render.hpp b/include/formula-cpp/trace_render.hpp index 7c11d23..0c12c9a 100644 --- a/include/formula-cpp/trace_render.hpp +++ b/include/formula-cpp/trace_render.hpp @@ -123,7 +123,8 @@ struct TraceRenderOptions /// value, a decimal, so it reads without `≈` in every style, its mode in /// brackets after it. A value the style cannot spell in its /// unit -- one padded or rounded in a unit whose declared decimals lie - /// outside -18 to 18, say -- reads `(not shown: ...)`, as a value its unit + /// outside the -18 to 18 text is written at, say -- reads + /// `(not shown: ...)`, as a value its unit /// cannot show does; a bound or a limit the author typed falls back to its /// exact fraction instead. /// diff --git a/test/calculation_trace_tests.cpp b/test/calculation_trace_tests.cpp index 49496de..74fb0af 100644 --- a/test/calculation_trace_tests.cpp +++ b/test/calculation_trace_tests.cpp @@ -212,7 +212,7 @@ inline constexpr auto sharesAdded = /// A vocabulary writing the typed share as the computed one is written, `s`. inline constexpr auto alikeShares = formula::vocabulary(formula::renames("s")); -/// A unit declaring more decimals than `DecimalPlaces` spans: 19. Invented. +/// A unit declaring more decimals than text is written at: 19. Invented. inline constexpr formula::Unit OverPrecise { .dimension = formula::dim::Length, .symbolText = formula::symbol("u"), .decimals = 19 }; @@ -1208,8 +1208,8 @@ TEST_CASE("a derived quantity standing for a typed value is rounded, as the valu TEST_CASE("a derivation's header says a value is not shown where its style cannot spell it", "[calculation][worksheet][trace][decimals]") { - // A length in a unit declaring 19 decimals, more than a rounding or a - // padding can take: under a style that pads or rounds, its header and + // A length in a unit declaring 19 decimals, more than text is written + // at: under a style that pads or rounds, its header and // the line reading it say it is not shown, while its root, a width in // millimetres scaled by a pure number and so in millimetres too, spells // it. In fractions every value is shown. diff --git a/test/format_tests.cpp b/test/format_tests.cpp index c25b3f6..6f1a2b3 100644 --- a/test/format_tests.cpp +++ b/test/format_tests.cpp @@ -37,7 +37,7 @@ struct Share: formula::Quantity { }; -/// A unit declaring more decimals than `DecimalPlaces` spans: 19. Invented. +/// A unit declaring more decimals than text is written at: 19. Invented. inline constexpr formula::Unit OverPrecise { .dimension = formula::dim::Length, .symbolText = formula::symbol("u"), .decimals = 19 }; @@ -45,7 +45,7 @@ struct OverPreciseLength: formula::Quantity Date: Thu, 8 Oct 2026 10:05:38 +0200 Subject: [PATCH 4/4] docs: rounding reaches 38 decimal places, text output 18 The numbers guide states the new rounding limit and recomputes the limits rational_from_double meets; the expressions guide gives the rounded logarithms' and exponential's places, the logarithms' 36-place ceiling and the exponential's threshold at -89, and the exact log10 range of 10^-38 to 10^38 the code has. The display guide and the numeric-headroom page keep 18 places as text output's own limit. The changelog records both, and MaximumDecimalPlaces. Refs #34 Signed-off-by: Christian Parpart --- CHANGELOG.md | 23 +++++++++++++++++++++ docs/dimensions.md | 5 +++-- docs/display.md | 10 +++++---- docs/expressions.md | 24 ++++++++++++++-------- docs/numbers.md | 44 ++++++++++++++++++++++++---------------- docs/numeric-headroom.md | 8 ++++---- 6 files changed, 78 insertions(+), 36 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e949d3d..ca383aa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,29 @@ change is recorded here. ## [Unreleased] +### Added + +- **`MaximumDecimalPlaces`** (`rounding.hpp`), 38: the most decimal places rounding takes either way, the scale + `Rational::from_decimal` and `_r` span. A place outside -38 to 38 is `ArithmeticError::Overflow`. + +### Changed + +- **Rounding reaches 38 decimal places in each direction**, up from 18, the largest power of ten in a 64-bit integer: + `checked_round` and `round` with `DecimalPlaces` or `SignificantDigits`, `checked_round_to_declared`, + `rational_from_double`, the rounding nodes, `rounded_output`, `rounded_sqrt`, `rounded_ln`, `rounded_log10` and + `rounded_exp`. 15 significant digits of a value below 1e-4 -- `-2.12345678901234e-5_r` is + `-212345678901234 / 10^19` -- now round instead of reporting `Overflow`. A place inside the range still reports + `Overflow` wherever an intermediate does not fit, never a wrong or clamped value: an irrational square root rounds + at no more than 19 places, and a logarithm the kernel computes at no more than 36 (ln 1 and log10 of a power of + ten are answered without the kernel, wherever their value fits: ln 1 at every places, log10 10 and log10 1/10 at + 38). The integer kernel behind the rounded logarithms and exponential works in 448 bits, up + from 384, and `rounded_exp` computes every argument from -89 up, where it used to answer below -43 without the + kernel. Each rounded logarithm or exponential the kernel computes costs about 10 to 15 percent more + constant-evaluation steps than before, so a constant expression close to the compiler's budget may need room. +- **Text output still stops at 18 decimal places** either way, as its own limit: `decimal_text`, `number_text`, a + trace's and `render()`'s numbers and `std::format` refuse more, as before, though the value was rounded further. + `checked_decimal_text` at fewer than -18 places is `Overflow` on its own, as it was through `checked_round`. + ## [0.5.0] - 2026-10-07 The fifth release. A tutorial on the documentation site builds a compressive-strength test step by step, and the README diff --git a/docs/dimensions.md b/docs/dimensions.md index 5d23aa9..15f25b7 100644 --- a/docs/dimensions.md +++ b/docs/dimensions.md @@ -629,8 +629,9 @@ value these fields state converts to one exactly. Conversion is built on `formula::Rational` and the `checked_` arithmetic functions, so it inherits their overflow behaviour and rounding limits -exactly -- including the ±18-decimal-place ceiling on `DecimalPlaces` and the -numerator-magnitude-dependent limits on rounding described in +exactly -- including the ±38-decimal-place limit on rounding, +`MaximumDecimalPlaces`, and the numerator-magnitude-dependent limits on +rounding described in [`docs/numbers.md`](numbers.md#limits). Nothing in `dimension.hpp` or `unit.hpp` widens or narrows those limits; a `Unit`'s magnitude and offset are exactly the integer pairs `Rational` already knows how to handle exactly. diff --git a/docs/display.md b/docs/display.md index aca774b..150a4f9 100644 --- a/docs/display.md +++ b/docs/display.md @@ -485,10 +485,12 @@ std::println("not measured: {}", absentText.view()); std::println("two places: {}\n", twoPlaces.view()); ``` -**When a number cannot be spelled.** `decimal_text` throws -`ArithmeticException` for more than 18 places, and where rounding to whole -tens or thousands overflows. `number_text` throws it where its rounding -overflows so, and for a padded or approximating style in a unit whose +**When a number cannot be spelled.** Text is written at no more than 18 +places either way -- its own limit, though `round` itself reaches 38 +([Numbers](numbers.md#limits)). `decimal_text` throws +`ArithmeticException` for more than 18 places, or fewer than -18, and where +rounding to whole tens or thousands overflows. `number_text` throws it +where its rounding overflows so, and for a padded or approximating style in a unit whose declared decimals lie outside -18 to 18: a `Rational` in any unit you pass, or a `Measured` whose unit has a symbol or is dimensionless. A `Measured` in a dimensioned unit with no symbol is shown in the coherent unit, at that diff --git a/docs/expressions.md b/docs/expressions.md index 81381ed..05ebd94 100644 --- a/docs/expressions.md +++ b/docs/expressions.md @@ -480,7 +480,7 @@ STATIC_REQUIRE(**formula::checked_evaluate_si(formula::log10( | Function | Exact at | Elsewhere, in `Rational` | |---|---|---| | `ln(x)` | x = 1: 0 | `Inexact` | -| `log10(x)` | x = 10^k, k from -18 to 18: k -- `1000` and `1/1000` alike | `Inexact` | +| `log10(x)` | x = 10^k, k from -38 to 38: k -- `1000` and `1/1000` alike | `Inexact` | | `exp(x)` | x = 0: 1 | `Inexact`, however large | `Inexact` is the exact layer refusing to approximate, as it refuses @@ -510,16 +510,22 @@ CHECK(lnAt(Rational { 2 }) == Ratio `"rounded_transcendental: ln 2 to 4 dp in every mode and of 1/2 with the directions paired the other way"`, whose `lnAt` helper evaluates `rounded_ln` at the given ratio.) -The places are the method's own, and at most 18; the result must fit a -`Rational` there, which any logarithm does: the `log10` of 10^18 - 1 is -reported to all 18 places. The integer kernel takes every argument a +The places are the method's own, from -38 to 38 as for `round`; the result +must fit a `Rational` there, which any logarithm does up to 36 places: the +`log10` of 10^18 - 1 is reported to 18 places, and `ln` 2 to 30, as +0.693147180559945309417232121458. The integer kernel takes every argument a `Rational` holds: `ln` of 2^127 - 1 is 88.029691931113054295 at 18 places, -under `Floor`. `exp` of more than 88.7 is `Overflow`, since e^x is then past -the largest `Rational`; below that it answers wherever the result fits the -declared places -- e^45 to 18 places, e^88 to whole units -- and is +under `Floor`. A logarithm's computation is up to 2^-120 wide, at least +1.6 · 10^-37, so at 37 or 38 places every logarithm the kernel computes is +`Overflow`. The exact points below -- `ln` 1, `log10` of a power of ten -- +are answered without the kernel, wherever their value fits: `ln` 1 at every +places, `log10` of 10 and of 1/10 at 38. +`exp` of more than 88.7 is `Overflow`, since e^x is then past the largest +`Rational`; below that it answers wherever the result fits the declared +places -- e^1 to 37 places, e^45 to 18, e^88 to whole units -- and is `Overflow` where it does not, as e^88.5 is at every places. Two kinds of argument never reach the kernel: a power of ten, 10^-38 up to 10^38, is -answered exactly, so `log10` of 10^30 is 30; and `exp` of less than -43 is +answered exactly, so `log10` of 10^30 is 30; and `exp` of less than -89 is 0, or one unit under `Ceiling` and `AwayFromZero`, whatever its width. Only ln 1, log10 10^k and exp 0 can tie, and the mode breaks the tie as `rounded<>` does: `log10` of @@ -539,7 +545,7 @@ guess. `rounded<...>(ln(x))` is not | argument failed | its error | its error | its error | | argument zero or below | `DomainError` | -- | `DomainError` (the logarithms) | | not a point where the value is rational | `Inexact` | `Inexact` | the rounded decimal | -| result too large at the declared places, or places outside -18 to 18 | -- | -- | `Overflow` | +| result too large at the declared places, or places outside -38 to 38 | -- | -- | `Overflow` | | rounding not decidable | -- | -- | `Overflow` | | `checked_evaluate_si` | `std::log`, `std::log10`; `DomainError` for zero, below and NaN | `std::exp`; too large is `+inf` | does not compile, as `rounded_sqrt` does not | diff --git a/docs/numbers.md b/docs/numbers.md index 9f7578f..ffbe72e 100644 --- a/docs/numbers.md +++ b/docs/numbers.md @@ -313,12 +313,17 @@ the answer. `Rational`'s numerator and denominator are `formula::Int128`, signed 128-bit integers: each holds up to 2^127 − 1, and a numerator down to −2^127 -- up to 39 decimal digits. `from_decimal`'s exponent, the `_r` literal and -`parse_decimal_text` reach 10^-38 to 10^38. `DecimalPlaces` and the -decimal-place form of `round` stay limited to ±18 places; an out-of-range -`DecimalPlaces` reports `Overflow`, while an out-of-range `SignificantDigits` -(fewer than 1) reports `DomainError` -- both mean "argument outside the -domain of the operation", but a caller switching on the code should expect -either one. +`parse_decimal_text` reach 10^-38 to 10^38, and so does rounding: `round` +and `checked_round` take places from -38 to 38, `MaximumDecimalPlaces`, +whether named as `DecimalPlaces` or reached through `SignificantDigits` -- 15 +significant digits of 2.12345678901234e-5 are rounded at 19 places. A place +outside that range reports `Overflow`, while an out-of-range +`SignificantDigits` (fewer than 1) reports `DomainError` -- both mean +"argument outside the domain of the operation", but a caller switching on the +code should expect either one. A place inside the range still reports +`Overflow` where the scaled value does not fit, as below. Text is another +matter: a number is written at no more than 18 places either way +([Displaying numbers](display.md)), however many it was rounded to. Rounding to `N` decimal places scales the value by `10^N`. Common factors of two cancel against the denominator first, so what must fit in `Rational::Int` @@ -330,17 +335,22 @@ is which for a power-of-two denominator is `|numerator| * 5^N`. **The limit is set by the numerator's magnitude**, not by the denominator and not by the size of -the value: `1 / 2^121` rounds correctly at all 18 places, while a 100-bit -numerator over the same denominator is refused at 18. +the value: `1 / 2^121` rounds correctly at all 38 places, while a 100-bit +numerator over the same denominator rounds at 11 places and is refused at 12. +Even a small value can be refused well inside the range: 2/3 at 38 places +forms 2 · 10^38, past 2^127. A `double` below 2^53 in magnitude has a numerator of at most 53 bits, and -rounding it at up to 18 places forms at most 2^53 · 5^18 · 2^18, below 2^113, -so such a value from `from_double_exact` or `rational_from_double` rounds at -every place from 0 to 18: 0,45 as a `double` is exactly -`8106479329266893 / 2^54`, and rounds to 18 places as 0.450000000000000011. -Past that the limit returns: a whole `double` such as 1e21 has a numerator of -its own magnitude and nothing to cancel, so it is refused at 18 places, 1e38 -even at 1, and `2^-100` at -18 places multiplies its denominator past 2^127. +rounding it at `N` places forms at most 2^53 · 10^N, below 2^127 up to +`N` = 22, so such a value from `from_double_exact` or `rational_from_double` +rounds at every place from 0 to 22; 2^53 - 1 is refused at 23. A fraction +whose denominator is at least 2^N cancels all of the 2^N and forms at most +2^53 · 5^N, which fits up to `N` = 31: 0,45 as a `double` is exactly +`8106479329266893 / 2^54`, rounds to 18 places as 0.450000000000000011 and +to 31 as 0.4500000000000000111022302462516, and is refused at 32. Past 2^53 +the limit returns: a whole `double` such as 1e21 has a numerator of its own +magnitude and nothing to cancel, so it is refused at 18 places, 1e38 even at +1, and `2^-100` at -9 places multiplies its denominator past 2^127. `from_decimal(45, -2)` is `9/20` -- the same nominal value, and what a method that writes 0,45 means. @@ -351,8 +361,8 @@ over `2^147`, is refused. For an exact decimal, prefer `from_decimal`: its numerator is whatever you passed -- usually a handful of significant digits -- so the limit above does not -bite. Its *denominator* need not be small at all: `from_decimal(1, -18)` is -`1/10^18`, and it still rounds correctly at every decimal place from 0 to 18. +bite. Its *denominator* need not be small at all: `from_decimal(1, -38)` is +`1/10^38`, and it still rounds correctly at every decimal place from 0 to 38. That is the clearest demonstration that the denominator is not what constrains rounding to decimals. diff --git a/docs/numeric-headroom.md b/docs/numeric-headroom.md index f3f2454..682d4b7 100644 --- a/docs/numeric-headroom.md +++ b/docs/numeric-headroom.md @@ -440,10 +440,10 @@ different denominator for every point is built to overflow. Only the integer `Rational` stores changed. These stay as they were: -- **Rounding's decimal places**, which stop at ±18 - ([Numbers](numbers.md#limits)). Widening them is a separate decision. - `from_decimal`'s exponents and the `_r` literal have since widened to - 10^±38 and a 128-bit mantissa. +- **Text output's decimal places**, which stop at 18 either way + ([Displaying numbers](display.md)). Rounding's own places, `from_decimal`'s + exponents and the `_r` literal have since widened to ±38 places, 10^±38 and + a 128-bit mantissa ([Numbers](numbers.md#limits)). - **The 64-bit fields** of `Unit`, `Band` and `Breakpoint`: `band` and `breakpoint` refuse a `Rational` bound or key that does not fit them.