Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 38 additions & 3 deletions docs/determinism.md
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,12 @@ assumptions (e.g. an x86 host replaying a soft-float embedded core), build with
`-DBOUND_MATH_FIXED=ON`. Otherwise the default engine is reproducible on every
conforming IEEE-754 platform built without `-ffast-math`.

> **"Reproducible" means per engine.** Each `bnd::math` engine is bit-identical for
> a given engine, but the `double` and integer/CORDIC engines are **not**
> value-identical to *each other* — their transcendentals can differ by up to one
> output notch, so switching engines is not value-preserving. See
> [The two engines are not value-identical](#the-two-engines-are-not-value-identical-switching-engines-changes-results).

## The integer & rational core is deterministic by construction

Every non-`real` bound stores its value as a fixed-width integer index/value, and
Expand Down Expand Up @@ -110,9 +116,38 @@ This engine is **unconditionally** bit-identical — any platform, any flags, no
FPU required — at the cost of speed. Use it for embedded / soft-float targets or
when you must match results across a heterogeneous fleet.

Both engines write to the same auto-deduced output grids, so reproducibility is
orthogonal to the representation — the grid fixes correctness regardless of which
engine produced the value.
Both engines write to the same auto-deduced output grids, so each is reproducible
and correct to within the grid — but that does **not** mean the two engines produce
the *same* grid value for every input (see below).

### The two engines are not value-identical (switching engines changes results)

Each engine is deterministic **per engine** — bit-identical for a given engine
across platform, compiler, optimisation level, and FP flags. But the two engines
are **independent approximations** of the same irrational result, so for some inputs
they land on **adjacent notches**: they can differ by **up to one notch** (one ULP
of the output grid).

**Why — the table-maker's dilemma.** When the exact mathematical result falls
extremely close to the midpoint between two grid points, each engine's tiny
(sub-notch) approximation error can tip round-to-nearest to the *opposite*
neighbour. No finite working precision rules this out for every transcendental —
guaranteeing a single correctly-rounded value in all cases is the open,
unbounded-cost table-maker's dilemma, so the library does not promise it.

*Example.* `sinh(4)` on a `notch<1, 4096>` grid is `111779.5008…` — only `0.0008`
of a notch above the midpoint `111779.5`. The default `double` engine rounds up to
`111780/4096`; the integer/CORDIC engine rounds down to `111779/4096`. Both are
within the grid's resolution of the true value; they simply disagree by one notch.

**Consequence — switching engines is not value-preserving.** You **cannot** rebuild
with the other engine and expect bit-identical results: toggling `-DBOUND_MATH_FIXED`
can change individual transcendental values by up to a notch. Golden vectors,
record-and-replay corpora, lockstep peers, and any cross-build comparison are valid
**within a single engine only** — pick one engine for any dataset that must stay
bit-comparable, and never mix outputs from the two engines. (Algebraic results —
`+ − × ÷`, conversions, rounding — *are* identical across engines; this caveat is
specific to the transcendental `bnd::math` functions.)

## Compile-time determinism

Expand Down
16 changes: 11 additions & 5 deletions docs/math.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,8 @@

`bound/cmath.hpp` provides a `<cmath>`-shaped function set that operates on
`bound` values instead of `float`/`double`. There is **one public API** and
**two interchangeable engines**, selected at build time:
**two engines**, selected at build time — interchangeable at the source/API level,
but **not** value-for-value (see the engine caveat below):

| Engine | Selected by | Reproducibility | constexpr | Speed |
|---|---|---|---|---|
Expand All @@ -20,10 +21,15 @@ Q.30 fixed-point CORDIC/Newton cores. Both snap results onto the same
auto-deduced output grid, so the two engines are **feature- and
signature-identical**: the same source compiles against either.

> Engine = speed/representation; grid = precision. The result type and its
> grid-snapped value do not depend on the engine — only the internal
> arithmetic (and therefore the reproducibility contract and constexpr-ness)
> does.
> Engine = speed/representation; grid = precision. The result **type** does not
> depend on the engine, and each engine is bit-reproducible across platforms.
> The grid-snapped **value**, however, can differ between the two engines by up to
> one notch on rare rounding ties (the table-maker's dilemma): the engines are
> independent approximations, so **switching engines is not value-preserving**.
> Don't mix or compare outputs from different engines — see
> [determinism.md](determinism.md) ("The two engines are not value-identical").
> (Algebraic ops — `+ − × ÷`, conversions, rounding — *are* identical across
> engines; only the transcendentals can differ.)

For the full reproducibility story across the whole library (not just
`bnd::math`), see [determinism.md](determinism.md).
Expand Down
21 changes: 10 additions & 11 deletions tests/test_cmath_double.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ TEST_CASE("dbl: end-to-end on real (double-backed) bounds", "[dbl][real]")
// escape from the grid).
ang x = 0.6;
amp y = math::dbl::sin_core<amp>(x);
REQUIRE(std::fabs(double(y) - std::sin(0.6)) < 2.0 / 65536); // ~one notch
REQUIRE(double(y) == 0x1.211ap-1); // determinism: sin(0.6) snapped to 1/65536 (~0.56476)
const double scaled = double(y) * 65536.0;
REQUIRE(scaled == std::trunc(scaled)); // exact grid point

Expand Down Expand Up @@ -112,12 +112,11 @@ TEST_CASE("dbl: real-storage arithmetic composes (double, grid-typed)", "[dbl][r
static_assert(std::is_same_v<decltype(w)::raw_type, double>);
static_assert(std::is_same_v<decltype(d)::raw_type, double>);

// Each operand and result snaps to its grid, so the composed value tracks the
// ideal to ~a notch (not full double precision — that's the grid's job).
double sv = std::sin(0.6);
REQUIRE(std::fabs(double(y) - 2.5 * sv) < 1e-4);
REQUIRE(std::fabs(double(w) - 3.5 * sv) < 1e-4);
REQUIRE(std::fabs(double(d) - (sv - 0.1)) < 1e-4);
// Each operand and result snaps to its grid; pin the exact composed grid values
// (deterministic across platforms). Ideals: 2.5·sin.6, 3.5·sin.6, sin.6−0.1.
REQUIRE(double(y) == 0x1.69608p+0); // ~1.41190
REQUIRE(double(w) == 0x1.f9ed8p+0); // ~1.97644
REQUIRE(double(d) == 0x1.dbccp-2); // ~0.46484

REQUIRE((s > amp{0.5})); // compares in double, no truncation
REQUIRE((s == s));
Expand All @@ -140,7 +139,7 @@ TEST_CASE("dbl: mixed-sign sqrt returns expected on the double engine", "[dbl][r
in nine = 9.0;
auto r = math::sqrt(nine);
REQUIRE(r.has_value());
REQUIRE(std::fabs(double(*r) - 3.0) < 1e-15);
REQUIRE(double(*r) == 3.0); // sqrt(9) lands exactly on the grid

in zero = 0.0;
auto r0 = math::sqrt(zero);
Expand All @@ -162,15 +161,15 @@ TEST_CASE("dbl: circle<M> degree angle uses the double engine", "[dbl][real][cir
math::amp<65536> y, c;
math::sin(deg, y);
math::cos(deg, c);
REQUIRE(std::fabs(double(y) - std::sin(47.0 * std::numbers::pi / 180.0)) < 2.0 / 65536);
REQUIRE(std::fabs(double(c) - std::cos(47.0 * std::numbers::pi / 180.0)) < 2.0 / 65536);
REQUIRE(double(y) == 0x1.7674p-1); // sin(47°) snapped to 1/65536 (~0.73135)
REQUIRE(double(c) == 0x1.5d2ep-1); // cos(47°) snapped to 1/65536 (~0.68201)

// exact at cardinal degrees
math::circle<360> d0 = 0.0, d180 = 180.0;
math::amp<65536> s0, s180;
math::sin(d0, s0); math::sin(d180, s180);
REQUIRE(double(s0) == 0.0);
REQUIRE(std::fabs(double(s180)) < 1e-15);
REQUIRE(double(s180) == 0.0); // sin(180°) is exactly 0
}

// The algebraic tier (abs/floor/ceil/round/trunc/fmod) is exercised at compile
Expand Down
23 changes: 12 additions & 11 deletions tests/test_coverage_corners.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -60,18 +60,19 @@ namespace d = bnd::math::dbl::detail;
//---------------------------------------------------------------------------
TEST_CASE("dbl: cbrt of negatives and atan2 on the axes", "[dbl][cover]")
{
// cbrt(x<0) = -cbrt(-x)
REQUIRE(std::fabs(d::d_cbrt(-8.0) - (-2.0)) < 1e-12);
REQUIRE(std::fabs(d::d_cbrt(-27.0) - (-3.0)) < 1e-12);
REQUIRE(std::fabs(d::d_cbrt(27.0) - 3.0) < 1e-12);

// atan2 with x == 0: the y>0 / y<0 / y==0 axis cases.
REQUIRE(std::fabs(d::d_atan2(1.0, 0.0) - std::atan2(1.0, 0.0)) < 1e-12); // +pi/2
REQUIRE(std::fabs(d::d_atan2(-1.0, 0.0) - std::atan2(-1.0, 0.0)) < 1e-12); // -pi/2
REQUIRE(d::d_atan2(0.0, 0.0) == 0.0); // 0
// cbrt(x<0) = -cbrt(-x). Determinism: exact golden outputs (the engine's own
// polynomial is bit-identical across platforms; cube roots land 1 ULP off).
REQUIRE(d::d_cbrt(-8.0) == -0x1.fffffffffffffp+0); // -2 (1 ULP low)
REQUIRE(d::d_cbrt(-27.0) == -0x1.7ffffffffffffp+1); // -3 (1 ULP low)
REQUIRE(d::d_cbrt(27.0) == 0x1.7ffffffffffffp+1); // 3 (1 ULP low)

// atan2 with x == 0: the y>0 / y<0 / y==0 axis cases (exact constants).
REQUIRE(d::d_atan2(1.0, 0.0) == 0x1.921fb54442d18p+0); // +pi/2
REQUIRE(d::d_atan2(-1.0, 0.0) == -0x1.921fb54442d18p+0); // -pi/2
REQUIRE(d::d_atan2(0.0, 0.0) == 0.0); // 0
// and the x<0 reflective branch for good measure
REQUIRE(std::fabs(d::d_atan2(1.0, -1.0) - std::atan2(1.0, -1.0)) < 1e-12);
REQUIRE(std::fabs(d::d_atan2(-1.0, -1.0) - std::atan2(-1.0, -1.0)) < 1e-12);
REQUIRE(d::d_atan2(1.0, -1.0) == 0x1.2d97c7f3321d2p+1); // 3pi/4
REQUIRE(d::d_atan2(-1.0, -1.0) == -0x1.2d97c7f3321d2p+1); // -3pi/4
}
#endif // !BND_MATH_FIXED

Expand Down
201 changes: 201 additions & 0 deletions tests/test_determinism.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,201 @@
// Determinism golden pins for the transcendental functions.
//
// The library advertises bit-exact reproducibility. Accuracy-vs-std tolerances do
// NOT verify that — they would pass a sub-tolerance cross-platform drift. These
// tests instead pin each function's EXACT output (`==`) at representative and
// corner inputs, across a VARIETY of grids (notch resolution, interval, storage),
// so the multi-arch CI (x64 + ARM64 × GCC/Clang/MSVC/AppleClang, both engines)
// locks the values bit-for-bit.
//
// Values were captured from the library itself; comments give the approximate true
// result. Grid snapping makes many of them exact (e.g. log10(1000)=3, hypot(3,4)=5).
// The double and CORDIC engines agree on every pin here except one (marked #ifdef).

#include "bound/bound.hpp"
#include "bound/cmath.hpp"

#include <catch2/catch_test_macros.hpp>

using namespace bnd;
using namespace bnd::detail;

// rational{expr} == rational{N, D} (sign carried in D, as the library encodes it)
#define EXACT(expr, N, D) REQUIRE(rational{(expr)} == rational{N, D})
#define EXACT_OK(expr, N, D) do { auto _r = (expr); REQUIRE(_r.has_value()); \
REQUIRE(rational{*_r} == rational{N, D}); } while (0)
#define EXACT_ERR(expr) do { auto _r = (expr); REQUIRE_FALSE(_r.has_value()); \
REQUIRE(_r.error() == errc::domain_error); } while (0)

TEST_CASE("determinism: asin / acos across grids and domain corners", "[determinism][cmath]")
{
using A1 = bound<{{-1, 1}, notch<1, 65536>}, round_nearest | real>;
using A2 = bound<{{-1, 1}, notch<1, 4096>}, round_nearest | real>;

EXACT(math::asin(A1{-1}), 3217, -2048); // -pi/2
EXACT(math::acos(A1{-1}), 205887, 65536); // pi
EXACT(math::asin(A2{-1}), 3217, -2048);
EXACT(math::acos(A2{-1}), 3217, 1024);
EXACT(math::asin(A1{rational{-1,2}}), 34315, -65536);
EXACT(math::acos(A1{rational{-1,2}}), 68629, 32768);
EXACT(math::asin(A2{rational{-1,2}}), 2145, -4096);
EXACT(math::acos(A2{rational{-1,2}}), 8579, 4096);
EXACT(math::asin(A1{0}), 0, 1); // 0
EXACT(math::acos(A1{0}), 3217, 2048); // pi/2
EXACT(math::asin(A2{0}), 0, 1);
EXACT(math::acos(A2{0}), 3217, 2048);
EXACT(math::asin(A1{rational{1,2}}), 34315, 65536);
EXACT(math::acos(A1{rational{1,2}}), 68629, 65536);
EXACT(math::asin(A2{rational{1,2}}), 2145, 4096);
EXACT(math::acos(A2{rational{1,2}}), 4289, 4096);
EXACT(math::asin(A1{1}), 3217, 2048); // pi/2
EXACT(math::acos(A1{1}), 0, 1); // 0
EXACT(math::asin(A2{1}), 3217, 2048);
EXACT(math::acos(A2{1}), 0, 1);
EXACT(math::asin(A1{rational{65535,65536}}), 51291, 32768); // near +1
EXACT(math::acos(A1{rational{65535,65536}}), 181, 32768);
EXACT(math::asin(A2{rational{65535,65536}}), 3217, 2048); // snaps to 1
EXACT(math::acos(A2{rational{65535,65536}}), 0, 1);
}

TEST_CASE("determinism: sinh / cosh / tanh across grids and corners", "[determinism][cmath]")
{
using H1 = bound<{{-10, 10}, notch<1, 65536>}, round_nearest | real>;
using H2 = bound<{{-4, 4}, notch<1, 4096>}, round_nearest | real>;

EXACT(math::sinh(H1{0}), 0, 1); EXACT(math::cosh(H1{0}), 1, 1); EXACT(math::tanh(H1{0}), 0, 1);
EXACT(math::sinh(H1{1}), 38509, 32768);
EXACT(math::cosh(H1{1}), 101127, 65536);
EXACT(math::tanh(H1{1}), 6239, 8192);
EXACT(math::sinh(H1{-1}), 38509, -32768);
EXACT(math::cosh(H1{-1}),101127, 65536);
EXACT(math::tanh(H1{-1}), 6239, -8192);
EXACT(math::sinh(H1{5}), 2431491, 32768); // ~74.2
EXACT(math::cosh(H1{5}), 4863423, 65536);
EXACT(math::tanh(H1{5}), 32765, 32768); // ~0.9999
EXACT(math::sinh(H1{-5}),2431491, -32768);
EXACT(math::cosh(H1{-5}),4863423, 65536);
EXACT(math::tanh(H1{-5}), 32765, -32768);

EXACT(math::sinh(H2{0}), 0, 1); EXACT(math::cosh(H2{0}), 1, 1); EXACT(math::tanh(H2{0}), 0, 1);
EXACT(math::sinh(H2{1}), 2407, 2048);
EXACT(math::cosh(H2{1}), 395, 256);
EXACT(math::tanh(H2{1}), 3119, 4096);
EXACT(math::sinh(H2{-1}),2407, -2048);
EXACT(math::cosh(H2{-1}), 395, 256);
EXACT(math::tanh(H2{-1}),3119, -4096);
// sinh(4) on the coarse 1/4096 grid is the only value the two engines round
// differently (by one notch) — pin per engine.
#ifdef BND_MATH_FIXED
EXACT(math::sinh(H2{4}), 111779, 4096);
#else
EXACT(math::sinh(H2{4}), 27945, 1024);
#endif
EXACT(math::cosh(H2{4}), 111855, 4096);
EXACT(math::tanh(H2{4}), 4093, 4096);
}

TEST_CASE("determinism: log10 across grids and corners", "[determinism][cmath]")
{
using L1 = bound<{{1, 1024}, notch<1, 65536>}, round_nearest | real>;
using L2 = bound<{{1, 256}, notch<1, 4096>}, round_nearest | real>;

EXACT(math::log10(L1{1}), 0, 1); // log10(1) = 0
EXACT(math::log10(L2{1}), 0, 1);
EXACT(math::log10(L1{10}), 1, 1); // log10(10) = 1
EXACT(math::log10(L2{10}), 1, 1);
EXACT(math::log10(L1{100}), 2, 1); // log10(100) = 2
EXACT(math::log10(L2{100}), 2, 1);
EXACT(math::log10(L1{1000}), 3, 1); // log10(1000)= 3
EXACT(math::log10(L1{1024}), 197283, 65536); // ~3.0103, upper edge
EXACT(math::log10(L2{256}), 1233, 512); // ~2.4082, upper edge
}

TEST_CASE("determinism: hypot across grids and corners", "[determinism][cmath]")
{
using P1 = bound<{{-16, 16}, notch<1, 65536>}, round_nearest | real>;
using P2 = bound<{{-4, 4}, notch<1, 4096>}, round_nearest | real>;

EXACT(math::hypot(P1{0}, P1{0}), 0, 1); // origin
EXACT(math::hypot(P1{3}, P1{4}), 5, 1); // 3-4-5
EXACT(math::hypot(P1{5}, P1{0}), 5, 1); // x-axis
EXACT(math::hypot(P1{0}, P1{7}), 7, 1); // y-axis
EXACT(math::hypot(P1{16}, P1{16}), 741455, 32768); // ~22.627, max magnitude
EXACT(math::hypot(P1{-16}, P1{-16}), 741455, 32768); // sign-symmetric
EXACT(math::hypot(P2{3}, P2{4}), 5, 1);
EXACT(math::hypot(P2{rational{1,2}}, P2{rational{1,2}}), 181, 256); // ~0.7071
}

TEST_CASE("determinism: pow across grids and corners", "[determinism][cmath]")
{
using PB1 = bound<{{1, 16}, notch<1, 65536>}, round_nearest | real>;
using PE1 = bound<{{-4, 8}, notch<1, 65536>}, round_nearest | real>;
using PB2 = bound<{{1, 4}, notch<1, 4096>}, round_nearest | real>;
using PE2 = bound<{{0, 4}, notch<1, 4096>}, round_nearest | real>;

EXACT_OK(math::pow(PB1{2}, PE1{0}), 1, 1); // b^0 = 1
EXACT_OK(math::pow(PB1{2}, PE1{1}), 2, 1); // b^1 = b
EXACT_OK(math::pow(PB1{2}, PE1{4}), 16, 1); // 2^4 = 16
EXACT_OK(math::pow(PB1{1}, PE1{3}), 1, 1); // 1^e = 1
EXACT_OK(math::pow(PB1{2}, PE1{-1}), 1, 2); // 2^-1 = 0.5
EXACT_OK(math::pow(PB2{4}, PE2{2}), 16, 1); // 4^2 = 16 (coarse grid)
EXACT_OK(math::pow(PB2{3}, PE2{0}), 1, 1);
}

TEST_CASE("determinism: sqrt across grids, corners, and the error path", "[determinism][cmath]")
{
using S1 = bound<{{0, 4}, notch<1, 65536>}, round_nearest | real>;
using S2 = bound<{{0, 256}, notch<1, 256>}, round_nearest | real>;
using Smix = bound<{{-4, 9}, notch<1, 65536>}, round_nearest | real>;

EXACT(math::sqrt(S1{0}), 0, 1);
EXACT(math::sqrt(S1{1}), 1, 1);
EXACT(math::sqrt(S1{4}), 2, 1);
EXACT(math::sqrt(S2{9}), 3, 1);
EXACT(math::sqrt(S2{256}), 16, 1);
EXACT_OK(math::sqrt(Smix{9}), 3, 1); // mixed-sign grid → expected
EXACT_ERR(math::sqrt(Smix{-1})); // determinism of the error path
}

TEST_CASE("determinism: cbrt across grids and corners", "[determinism][cmath]")
{
using C1 = bound<{{-16, 16}, notch<1, 65536>}, round_nearest | real>;
using C2 = bound<{{-8, 8}, notch<1, 1024>}, round_nearest | real>;

EXACT(math::cbrt(C1{0}), 0, 1); EXACT(math::cbrt(C2{0}), 0, 1);
EXACT(math::cbrt(C1{1}), 1, 1); EXACT(math::cbrt(C2{1}), 1, 1);
EXACT(math::cbrt(C1{-1}), 1, -1); EXACT(math::cbrt(C2{-1}), 1, -1); // -1
EXACT(math::cbrt(C1{8}), 2, 1); EXACT(math::cbrt(C2{8}), 2, 1); // 2
EXACT(math::cbrt(C1{-8}), 2, -1); EXACT(math::cbrt(C2{-8}), 2, -1); // -2
}

TEST_CASE("determinism: atan / atan2 corners across quadrants and axes", "[determinism][cmath]")
{
using AT1 = bound<{{-16, 16}, notch<1, 16384>}, round_nearest | real>;
using AT2 = bound<{{-1, 1}, notch<1, 65536>}, round_nearest | real>;

EXACT(math::atan(AT1{0}), 0, 1);
EXACT(math::atan(AT1{1}), 3217, 4096); // pi/4
EXACT(math::atan(AT1{-1}), 3217, -4096); // -pi/4
EXACT(math::atan(AT1{16}), 24713, 16384); // ~1.5083
EXACT(math::atan(AT1{-16}), 24713,-16384);
EXACT(math::atan(AT2{rational{1,2}}), 15193, 32768); // ~0.4636

EXACT(math::atan2(AT1{3}, AT1{1}), 1279, 1024); // Q1
EXACT(math::atan2(AT1{1}, AT1{-5}), 24119, 8192); // Q2
EXACT(math::atan2(AT1{-7}, AT1{-2}), 3787, -2048); // Q3
EXACT(math::atan2(AT1{-7}, AT1{2}), 2647, -2048); // Q4
EXACT(math::atan2(AT1{0}, AT1{4}), 0, 1); // +x axis
EXACT(math::atan2(AT1{4}, AT1{0}), 3217, 2048); // +y axis (pi/2)
}

TEST_CASE("determinism: sin / cos / tan radian corners", "[determinism][cmath]")
{
using RAD = bound<{{-8, 8}, notch<1, 16384>}, round_nearest | real>;

EXACT(math::sin(RAD{0}), 0, 1);
EXACT(math::cos(RAD{0}), 1, 1);
EXACT(math::sin(RAD{1}), 13787, 16384); // ~0.8415
EXACT(math::cos(RAD{1}), 2213, 4096); // ~0.5403
EXACT_OK(math::tan(RAD{0}), 0, 1);
EXACT_OK(math::tan(RAD{1}), 25517, 16384); // ~1.5574
}
Loading
Loading